5.13.3-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!> - (35) SO: SM: 1.0 - SPAEF of spatially distributed soil moisture
22!> - (36) SO: ET: 1.0 - SPAEF of spatially distributed evapotranspiration
23!> - (37) SO: SM + ET: 1.0 - 0.5 * [SPAEF(SM) + SPAEF(ET)]
24!> - (38) SO: SM + Q : 1.0 - [input wt.* SPAEF(SM) + SPAEF(Q)]/(input wt. + 1)
25!> - (39) SO: ET + Q : 1.0 - [input wt.* SPAEF(ET) + SPAEF(Q)]/(input wt. + 1)
26!> - (40) SO: SM + ET + Q: 1.0 - [2.5* SPAEF(SM) + 2.5* SPAEF(ET) + SPAEF(Q)]/6
27!> - (41) SO: ET: 1.0 - Pattern dissimilarity (PD) of spatially distributed soil moisture
28!> - (42) SO: SM + ET: 1.0 - 0.5 * [PD(SM) + PD(ET)]
29!> - (43) SO: SM + Q : 1.0 - [input wt.* PD(SM) + PD(Q)]/(input wt. + 1)
30!> - (44) SO: ET + Q : 1.0 - [input wt.* PD(ET) + PD(Q)]/(input wt. + 1)
31!> - (45) SO: SM + ET + Q: 1.0 - [0.5*input wt.* PD(SM) + 0.5*input wt.* PD(ET) + SPAEF(Q)]/(input wt. + 1)
32!> - (46) SO: SPF: 1.0 - CA(SPF)
33!> - (47) SO: SPF+ Q : 1.0 - [input wt.* CA(SPF) + KGE(Q)]/(input wt. + 1)
34!> \changelog
35!! - Oldrich Rakovec Oct 2015
36!! - added obj. func. 35 to 47, related to spatial checks (SPAEF, CA, PD)
37!! - Oldrich Rakovec Oct 2015
38!! - added obj. func. 15 (objective_kge_q_rmse_tws) and extract_domain_avg_tws routine, former basin_avg
39!! - Robert Schweppe Jun 2018
40!! - refactoring and reformatting
41!> \authors Juliane Mai
42!> \date Dec 2012
43!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
44!! mHM is released under the LGPLv3+ license \license_note
45!> \ingroup f_mhm
46MODULE mo_objective_function
47
48 ! This module provides objective functions for optimization of the UFZ CHS mesoscale hydrologic model mHM.
49
50 ! Written Juliane Mai, Dec 2012
51 ! Modified Stephan Thober, Oct 2015 moved all runoff only related objectives to mRM
52
53 USE mo_kind, ONLY : i4, dp
54 use mo_optimization_utils, only : eval_interface
55 use mo_message, only : message, error_message
56
57 IMPLICIT NONE
58
59 PRIVATE
60
61 PUBLIC :: objective
62#ifdef MPI
63 PUBLIC :: objective_master, objective_subprocess ! objective function wrapper for soil moisture only
64#endif
65
66 ! ------------------------------------------------------------------
67
68CONTAINS
69
70 ! ------------------------------------------------------------------
71
72 ! NAME
73 ! objective
74
75 ! PURPOSE
76 !> \brief Wrapper for objective functions.
77
78 !> \details The functions selects the objective function case defined in a namelist,
79 !> i.e. the global variable \e opti\_function.
80 !> It return the objective function value for a specific parameter set.
81
82 ! INTENT(IN)
83 !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
84 !> \param[in] "procedure(eval_interface) :: eval"
85
86 ! INTENT(IN), OPTIONAL
87 !> \param[in] "real(dp), optional :: arg1"
88
89 ! INTENT(OUT), OPTIONAL
90 !> \param[out] "real(dp), optional :: arg2"
91 !> \param[out] "real(dp), optional :: arg3"
92
93 ! RETURN
94 !> \return real(dp) :: objective — objective function value
95 !> (which will be e.g. minimized by an optimization routine like DDS)
96
97 ! HISTORY
98 !> \authors Juliane Mai
99
100 !> \date Dec 2012
101
102 ! Modifications:
103 ! Stephan Thober Oct 2015 - moved all runoff related objective functions to mRM
104 ! Robert Schweppe Jun 2018 - refactoring and reformatting
105
106 FUNCTION objective(parameterset, eval, arg1, arg2, arg3)
107
108 use mo_common_constants, only : nodata_dp
109 use mo_common_mhm_mrm_variables, only : opti_function
110
111 implicit none
112
113 REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
114
115 procedure(eval_interface), INTENT(IN), POINTER :: eval
116
117 real(dp), optional, intent(in) :: arg1
118
119 real(dp), optional, intent(out) :: arg2
120
121 real(dp), optional, intent(out) :: arg3
122
123 REAL(dp) :: objective
124
125 real(dp), dimension(6) :: multiple_objective
126
127
128 if (present(arg1) .or. present(arg2) .or. present(arg3)) then
129 call error_message("Error mo_objective_function: Received unexpected argument, check optimization settings")
130 end if
131
132 ! set these to nan so compiler does not complain
133 if (present(arg2)) then
134 arg2 = nodata_dp
135 end if
136 if (present(arg3)) then
137 arg3 = nodata_dp
138 end if
139
140 select case (opti_function)
141 case (10)
142 ! KGE of catchment average SM
143 objective = objective_sm_kge_catchment_avg(parameterset, eval)
144 case (11)
145 ! pattern dissimilarity (PD) of SM fields
146 objective = objective_sm_pd(parameterset, eval)
147 case (12)
148 ! sum of squared errors of standard_score SM
149 objective = objective_sm_sse_standard_score(parameterset, eval)
150 case (13)
151 ! soil moisture correlation - temporal
152 objective = objective_sm_corr(parameterset, eval)
153 case (15)
154 ! KGE for Q * RMSE for domain_avg TWS (standarized scored)
155 objective = objective_kge_q_rmse_tws(parameterset, eval)
156 case (17)
157 ! KGE of catchment average SM
158 objective = objective_neutrons_kge_catchment_avg(parameterset, eval)
159 case (27)
160 ! KGE of catchment average ET
161 objective = objective_et_kge_catchment_avg(parameterset, eval)
162 case (28)
163 ! KGE for Q + SSE for SM (standarized scored)
164 objective = objective_kge_q_sm_corr(parameterset, eval)
165 case (29)
166 ! KGE for Q + KGE of catchment average ET
167 objective = objective_kge_q_et(parameterset, eval)
168 case (30)
169 ! KGE for Q * RMSE for domain_avg ET (standarized scored)
170 objective = objective_kge_q_rmse_et(parameterset, eval)
171 case (33)
172 multiple_objective = objective_q_et_tws_kge_catchment_avg(parameterset, eval)
173 objective = multiple_objective(1)
174 case (34)
175 ! KGE for Q * Absolute-Error for BFI
176 objective = objective_kge_q_bfi(parameterset, eval)
177 case (35)
178 ! SPAEF for SM fields
179 objective = objective_spaef_sm(parameterset, eval)
180 case (36)
181 ! SPAEF for ET fields
182 objective = objective_spaef_et(parameterset, eval)
183 case (37)
184 ! SPAEF for SM fields * SPAEF for ET fields
185 objective = objective_spaef_sm_spaef_et(parameterset, eval)
186 case (38)
187 ! SPAEF for SM fields * KGE for Q
188 objective = objective_spaef_sm_kge_q(parameterset, eval)
189 case (39)
190 ! SPAEF for ET fields * KGE for Q
191 objective = objective_spaef_et_kge_q(parameterset, eval)
192 case (40)
193 ! SPAEF for SM fields * SPAEF for ET fields * KGE for Q
194 objective = objective_spaef_sm_spaef_et_kge_q(parameterset, eval)
195 case (41)
196 ! pattern dissimilarity (PD) of ET fields
197 objective = objective_pd_et(parameterset, eval)
198 case (42)
199 ! PD for SM fields * PD for ET fields
200 objective = objective_pd_sm_pd_et(parameterset, eval)
201 case (43)
202 ! PD for SM fields * KGE for Q
203 objective = objective_pd_sm_kge_q(parameterset, eval)
204 case (44)
205 ! PD for ET fields * KGE for Q
206 objective = objective_pd_et_kge_q(parameterset, eval)
207 case (45)
208 ! PD for SM fields * PD for ET fields * KGE for Q
209 objective = objective_pd_sm_pd_et_kge_q(parameterset, eval)
210 case (46)
211 ! Classification accuracy for SPF fields
212 objective = objective_ca_spf(parameterset, eval)
213 case (47)
214 ! Classification accuracy for SPF fields * KGE for Q
215 objective = objective_ca_spf_kge_q(parameterset, eval)
216
217 case default
218 call error_message("Error objective: opti_function not implemented yet.")
219 end select
220
221 END FUNCTION objective
222
223 ! ------------------------------------------------------------------
224
225 ! NAME
226 ! objective_master
227
228 ! PURPOSE
229 !> \brief Wrapper for objective functions.
230
231 !> \details The functions sends the parameterset to the subprocess, receives
232 !> the partial objective and calculates the final objective
233 ! INTENT(IN)
234 !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
235 !> \param[in] "procedure(eval_interface) :: eval"
236
237 ! INTENT(IN), OPTIONAL
238 !> \param[in] "real(dp), optional :: arg1"
239
240 ! INTENT(OUT), OPTIONAL
241 !> \param[out] "real(dp), optional :: arg2"
242 !> \param[out] "real(dp), optional :: arg3"
243
244 ! RETURN
245 !> \return real(dp) :: objective — objective function value
246 !> (which will be e.g. minimized by an optimization routine like DDS)
247
248 ! HISTORY
249 !> \authors Juliane Mai
250
251 !> \date Dec 2012
252
253 ! Modifications:
254 ! Stephan Thober Oct 2015 - moved all runoff related objective functions to mRM
255 ! Robert Schweppe Jun 2018 - refactoring and reformatting
256 ! Maren Kaluza Jun 2019 - parallel version
257
258#ifdef MPI
259 FUNCTION objective_master(parameterset, eval, arg1, arg2, arg3)
260
261 use mo_common_constants, only : nodata_dp
262 use mo_common_mhm_mrm_variables, only : opti_function
263 use mo_common_mpi_tools, only : distribute_parameterset
264 use mo_common_variables, only : domainmeta
265 use mo_string_utils, only : num2str
266 use mpi_f08
267
268 implicit none
269
270 REAL(dp), DIMENSION(:), INTENT(IN) :: parameterset
271
272 procedure(eval_interface), INTENT(IN), POINTER :: eval
273
274 real(dp), optional, intent(in) :: arg1
275
276 real(dp), optional, intent(out) :: arg2
277
278 real(dp), optional, intent(out) :: arg3
279
280 REAL(dp) :: objective_master
281
282 REAL(dp) :: partial_objective
283
284 real(dp), dimension(6) :: multiple_partial_objective
285
286 real(dp), dimension(6) :: multiple_master_objective
287
288 ! for sixth root
289 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
290
291 integer(i4) :: iproc, nproc
292
293 integer(i4) :: ierror
294
295 type(mpi_status) :: status
296
297
298 if (present(arg1) .or. present(arg2) .or. present(arg3)) then
299 call error_message("Error mo_objective_function: Received unexpected argument, check optimization settings")
300 end if
301
302 ! set these to nan so compiler does not complain
303 if (present(arg2)) then
304 arg2 = nodata_dp
305 end if
306 if (present(arg3)) then
307 arg3 = nodata_dp
308 end if
309
310 call distribute_parameterset(parameterset)
311 select case (opti_function)
312 case (10 : 13, 17, 27 : 29)
313 call mpi_comm_size(domainmeta%comMaster, nproc, ierror)
314 objective_master = 0.0_dp
315 do iproc = 1, nproc - 1
316 call mpi_recv(partial_objective, 1, mpi_double_precision, iproc, 0, domainmeta%comMaster, status, ierror)
317 objective_master = objective_master + partial_objective
318 end do
319 objective_master = objective_master**onesixth
320 case (15)
321 ! KGE for Q * RMSE for domain_avg TWS (standarized scored)
322 call error_message("case 15, objective_kge_q_rmse_tws not implemented in parallel yet")
323 case (30)
324 ! KGE for Q * RMSE for domain_avg ET (standarized scored)
325 ! objective_master = objective_kge_q_rmse_et(parameterset, eval)
326 call message("case 30, objective_kge_q_rmse_et not implemented in parallel yet")
327 case(33)
328 call mpi_comm_size(domainmeta%comMaster, nproc, ierror)
329 objective_master = 0.0_dp
330 multiple_master_objective(:) = 0.0_dp
331 do iproc = 1, nproc - 1
332 call mpi_recv(multiple_partial_objective, 6, mpi_double_precision, iproc, 0, domainmeta%comMaster, status, ierror)
333 multiple_master_objective = multiple_master_objective + multiple_partial_objective
334 end do
335 objective_master = objective_master + &
336 (multiple_master_objective(1)+multiple_master_objective(2)+multiple_master_objective(3))
337 objective_master = (objective_master/multiple_master_objective(4))**onesixth
338
339 case default
340 call error_message("Error objective_master: opti_function not implemented yet.")
341 end select
342
343 select case (opti_function)
344 case(10)
345 call message(' objective_sm_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
346 case(11)
347 call message(' objective_sm_pd = ', num2str(objective_master, '(F9.5)'))
348 case(12)
349 call message(' objective_sm_sse_standard_score = ', num2str(objective_master, '(E12.5)'))
350 case(13)
351 call message(' objective_sm_corr = ', num2str(objective_master, '(F9.5)'))
352 case(17)
353 call message(' objective_neutrons_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
354 case(27)
355 call message(' objective_et_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
356 case(28)
357 call message(' objective_kge_q_sm_corr = ', num2str(objective_master, '(F9.5)'))
358 case(29)
359 call message(' objective_kge_q_et = ', num2str(objective_master, '(F9.5)'))
360 case(33)
361 call message(' objective_q_et_tws_kge_catchment_avg = ', num2str(objective_master, '(F9.5)'))
362 case(35)
363 call message(' objective_spaef_sm = ', num2str(objective_master, '(F9.5)'))
364 case(36)
365 call message(' objective_spaef_et = ', num2str(objective_master, '(F9.5)'))
366 case(37)
367 call message(' objective_spaef_sm_spaef_et = ', num2str(objective_master, '(F9.5)'))
368 case(38)
369 call message(' objective_spaef_sm_kge_q = ', num2str(objective_master, '(F9.5)'))
370 case(39)
371 call message(' objective_spaef_et_kge_q = ', num2str(objective_master, '(F9.5)'))
372 case(40)
373 call message(' objective_spaef_sm_spaef_et_kge_q = ', num2str(objective_master, '(F9.5)'))
374 case(41)
375 call message(' objective_pd_et = ', num2str(objective_master, '(F9.5)'))
376 case(42)
377 call message(' objective_pd_sm_pd_et = ', num2str(objective_master, '(F9.5)'))
378 case(43)
379 call message(' objective_pd_sm_kge_q = ', num2str(objective_master, '(F9.5)'))
380 case(44)
381 call message(' objective_pd_et_kge_q = ', num2str(objective_master, '(F9.5)'))
382 case(45)
383 call message(' objective_pd_sm_pd_et_kge_q = ', num2str(objective_master, '(F9.5)'))
384 case(46)
385 call message(' objective_ca_spf = ', num2str(objective_master, '(F9.5)'))
386 case(47)
387 call message(' objective_ca_spf_kge_q = ', num2str(objective_master, '(F9.5)'))
388 case default
389 call error_message("Error objective_master: opti_function not implemented yet, this part of the code should never execute.")
390 end select
391
392 END FUNCTION objective_master
393
394 ! ------------------------------------------------------------------
395
396 ! NAME
397 ! objective_subprocess
398
399 ! PURPOSE
400 !> \brief Wrapper for objective functions.
401
402 !> \details The function receives the parameterset from the master
403 !> process, selects the objective function case defined in a namelist,
404 !> i.e. the global variable \e opti\_function.
405 !> It returns the partial objective function value for a specific parameter set.
406
407 ! INTENT(IN)
408 !> \param[in] "REAL(dp), DIMENSION(:) :: parameterset"
409 !> \param[in] "procedure(eval_interface) :: eval"
410
411 ! INTENT(IN), OPTIONAL
412 !> \param[in] "real(dp), optional :: arg1"
413
414 ! INTENT(OUT), OPTIONAL
415 !> \param[out] "real(dp), optional :: arg2"
416 !> \param[out] "real(dp), optional :: arg3"
417
418 ! RETURN
419 !> \return real(dp) :: objective — objective function value
420 !> (which will be e.g. minimized by an optimization routine like DDS)
421
422 ! HISTORY
423 !> \authors Juliane Mai
424
425 !> \date Dec 2012
426
427 ! Modifications:
428 ! Stephan Thober Oct 2015 - moved all runoff related objective functions to mRM
429 ! Robert Schweppe Jun 2018 - refactoring and reformatting
430 ! Maren Kaluza Jun 2019 - parallel version
431
432 subroutine objective_subprocess(eval, arg1, arg2, arg3)
433
434 use mo_common_constants, only : nodata_dp
435 use mo_common_mhm_mrm_variables, only : opti_function
436 use mo_common_mpi_tools, only : get_parameterset
437 use mo_common_variables, only : domainmeta
438 use mpi_f08
439
440 implicit none
441
442 procedure(eval_interface), INTENT(IN), POINTER :: eval
443
444 real(dp), optional, intent(in) :: arg1
445
446 real(dp), optional, intent(out) :: arg2
447
448 real(dp), optional, intent(out) :: arg3
449
450 REAL(dp) :: partial_objective
451
452 real(dp), dimension(6) :: multiple_partial_objective
453
454 REAL(dp), DIMENSION(:), allocatable :: parameterset
455
456 integer(i4) :: ierror
457
458 type(mpi_status) :: status
459
460 logical :: do_obj_loop
461
462 do ! a do loop without condition runs until exit
463 call mpi_recv(do_obj_loop, 1, mpi_logical, 0, 0, domainmeta%comMaster, status, ierror)
464
465 if (.not. do_obj_loop) exit
466
467 if (present(arg1) .or. present(arg2) .or. present(arg3)) then
468 call error_message("Error mo_objective_function: Received unexpected argument, check optimization settings")
469 end if
470
471 ! set these to nan so compiler does not complain
472 if (present(arg2)) then
473 arg2 = nodata_dp
474 end if
475 if (present(arg3)) then
476 arg3 = nodata_dp
477 end if
478 call get_parameterset(parameterset)
479 select case (opti_function)
480 case (10)
481 ! KGE of catchment average SM
482 partial_objective = objective_sm_kge_catchment_avg(parameterset, eval)
483 case (11)
484 ! pattern dissimilarity (PD) of SM fields
485 partial_objective = objective_sm_pd(parameterset, eval)
486 case (12)
487 ! sum of squared errors of standard_score SM
488 partial_objective = objective_sm_sse_standard_score(parameterset, eval)
489 case (13)
490 ! soil moisture correlation - temporal
491 partial_objective = objective_sm_corr(parameterset, eval)
492 case (15)
493 ! KGE for Q * RMSE for domain_avg TWS (standarized scored)
494 ! partial_objective = objective_kge_q_rmse_tws(parameterset, eval)
495 call error_message("Error objective_subprocess: case 15 not supported with MPI.")
496 case (17)
497 ! KGE of catchment average SM
498 partial_objective = objective_neutrons_kge_catchment_avg(parameterset, eval)
499 case (27)
500 ! KGE of catchment average ET
501 partial_objective = objective_et_kge_catchment_avg(parameterset, eval)
502 case (28)
503 ! KGE for Q + SSE for SM (standarized scored)
504 partial_objective = objective_kge_q_sm_corr(parameterset, eval)
505 case (29)
506 ! KGE for Q + KGE of catchment average ET
507 partial_objective = objective_kge_q_et(parameterset, eval)
508 case (30)
509 ! KGE for Q * RMSE for domain_avg ET (standarized scored)
510 ! partial_objective = objective_kge_q_rmse_et(parameterset, eval)
511 call error_message("Error objective_subprocess: case 30 not supported with MPI.")
512 case(33)
513 multiple_partial_objective = objective_q_et_tws_kge_catchment_avg(parameterset, eval)
514 case (35)
515 ! SPAEF for SM fields
516 partial_objective = objective_spaef_sm(parameterset, eval)
517 case (36)
518 ! SPAEF for ET fields
519 partial_objective = objective_spaef_et(parameterset, eval)
520 case (37)
521 ! SPAEF for SM fields * SPAEF for ET fields
522 partial_objective = objective_spaef_sm_spaef_et(parameterset, eval)
523 case (38)
524 ! SPAEF for SM fields * KGE for Q
525 partial_objective = objective_spaef_sm_kge_q(parameterset, eval)
526 case (39)
527 ! SPAEF for ET fields * KGE for Q
528 partial_objective = objective_spaef_et_kge_q(parameterset, eval)
529 case (40)
530 ! SPAEF for SM fields * SPAEF for ET fields * KGE for Q
531 partial_objective = objective_spaef_sm_spaef_et_kge_q(parameterset, eval)
532 case (41)
533 ! pattern dissimilarity (PD) for ET fields
534 partial_objective = objective_pd_et(parameterset, eval)
535 case (42)
536 ! PD for SM fields * PD for ET fields
537 partial_objective = objective_pd_sm_pd_et(parameterset, eval)
538 case (43)
539 ! PD for SM fields * KGE for Q
540 partial_objective = objective_pd_sm_kge_q(parameterset, eval)
541 case (44)
542 ! PD for ET fields * KGE for Q
543 partial_objective = objective_pd_et_kge_q(parameterset, eval)
544 case (45)
545 ! PD for SM fields * PD for ET fields * KGE for Q
546 partial_objective = objective_pd_sm_pd_et_kge_q(parameterset, eval)
547 case (46)
548 ! Classification accuracy for SPF fields
549 partial_objective = objective_ca_spf(parameterset, eval)
550 case (47)
551 ! Classification accuracy for SPF fields * KGE for Q
552 partial_objective = objective_ca_spf_kge_q(parameterset, eval)
553 case default
554 call error_message("Error objective_subprocess: opti_function not implemented yet.")
555 end select
556
557 select case (opti_function)
558 case (10 : 13, 17, 27 : 29, 35 : 47)
559 call mpi_send(partial_objective,1, mpi_double_precision,0,0,domainmeta%comMaster,ierror)
560 case(33)
561 call mpi_send(multiple_partial_objective, 6, mpi_double_precision,0,0,domainmeta%comMaster,ierror)
562 case default
563 call error_message("Error objective_subprocess: this part should not be executed -> error in the code.")
564 end select
565
566 deallocate(parameterset)
567 end do
568
569 END subroutine objective_subprocess
570
571#endif
572 ! ------------------------------------------------------------------
573
574 ! NAME
575 ! objective_sm_kge_catchment_avg
576
577 ! PURPOSE
578 !> \brief Objective function for soil moisture.
579
580 !> \details The objective function only depends on a parameter vector.
581 !> The model will be called with that parameter vector and
582 !> the model output is subsequently compared to observed data.
583
584 !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average
585 !> soil mloisture (SM) is calculated
586 !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
587 !> where
588 !> \f$ r \f$ = Pearson product-moment correlation coefficient
589 !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM
590 !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
591 !> is calculated and the objective function for a given domain \f$ i \f$ is
592 !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
593 !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
594 !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
595
596 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
597 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
598 !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f]
599 !> The observed data L1_sm, L1_sm_mask are global in this module.
600
601 ! INTENT(IN)
602 !> \param[in] "real(dp), dimension(:) :: parameterset"
603 !> \param[in] "procedure(eval_interface) :: eval"
604
605 ! RETURN
606 !> \return real(dp) :: objective_sm_kge_catchment_avg — objective function value
607 !> (which will be e.g. minimized by an optimization routine like DDS)
608
609 ! HISTORY
610 !> \authors Matthias Zink
611
612 !> \date May 2015
613
614 ! Modifications:
615 ! Robert Schweppe Jun 2018 - refactoring and reformatting
616
617 FUNCTION objective_sm_kge_catchment_avg(parameterset, eval)
618
619 use mo_optimization_types, only : optidata_sim
620 use mo_common_constants, only : nodata_dp
621 use mo_common_variables, only : level1, domainmeta
622 use mo_errormeasures, only : kge
623 use mo_global_variables, only : l1_smobs
624 use mo_moment, only : average
625 use mo_string_utils, only : num2str
626
627 implicit none
628
629 real(dp), dimension(:), intent(in) :: parameterset
630
631 procedure(eval_interface), INTENT(IN), POINTER :: eval
632
633 real(dp) :: objective_sm_kge_catchment_avg
634
635 ! domain loop counter
636 integer(i4) :: idomain
637
638 ! time loop counter
639 integer(i4) :: itime
640
641 ! number of time steps in simulated SM
642 integer(i4) :: n_time_steps
643
644 ! ncells1 of level 1
645 integer(i4) :: ncells1
646
647 ! number of invalid timesteps
648 real(dp) :: invalid_times
649#ifndef MPI
650 ! for sixth root
651 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
652#endif
653
654 ! spatial average of observed soil moisture
655 real(dp), dimension(:), allocatable :: sm_catch_avg_domain
656
657 ! spatial avergae of modeled soil moisture
658 real(dp), dimension(:), allocatable :: sm_opti_catch_avg_domain
659
660 type(optidata_sim), dimension(:), allocatable :: smoptisim
661
662 ! mask for valid sm catchment avg time steps
663 logical, dimension(:), allocatable :: mask_times
664
665
666 allocate(smoptisim(domainmeta%nDomains))
667 call eval(parameterset, smoptisim = smoptisim)
668
669 ! initialize some variables
670 objective_sm_kge_catchment_avg = 0.0_dp
671
672 ! loop over domain - for applying power law later on
673 do idomain = 1, domainmeta%nDomains
674
675 ! get domain information
676 ncells1 = level1(idomain)%ncells
677
678 ! allocate
679 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
680 allocate(sm_catch_avg_domain(size(smoptisim(idomain)%dataSim, dim = 2)))
681 allocate(sm_opti_catch_avg_domain(size(smoptisim(idomain)%dataSim, dim = 2)))
682
683 ! initalize
684 mask_times = .true.
685 sm_catch_avg_domain = nodata_dp
686 sm_opti_catch_avg_domain = nodata_dp
687
688 invalid_times = 0.0_dp
689 ! calculate catchment average soil moisture
690 n_time_steps = size(smoptisim(idomain)%dataSim, dim = 2)
691 do itime = 1, n_time_steps
692
693 ! check for enough data points in timesteps for KGE calculation
694 ! more then 10 percent avaiable in current field
695 if (count(l1_smobs(idomain)%maskObs(:, itime)) .LE. (0.10_dp * real(ncells1, dp))) then
696 invalid_times = invalid_times + 1.0_dp
697 mask_times(itime) = .false.
698 cycle
699 end if
700 sm_catch_avg_domain(itime) = average(l1_smobs(idomain)%dataObs(:, itime), mask = l1_smobs(idomain)%maskObs(:, itime))
701 sm_opti_catch_avg_domain(itime) = average(smoptisim(idomain)%dataSim(:, itime), mask = l1_smobs(idomain)%maskObs(:, itime))
702 end do
703
704 ! user information about invalid times
705 if (invalid_times .GT. 0.5_dp) then
706 call message(.LT.' WARNING: objective_sm: Detected invalid timesteps ( 10 valid data points).')
707 call message(' Fraction of invalid timesteps: ', &
708 num2str(invalid_times / real(n_time_steps, dp), '(F4.2)'))
709 end if
710
711
712 ! calculate average soil moisture KGE over all domains with power law
713 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
714 objective_sm_kge_catchment_avg = objective_sm_kge_catchment_avg + &
715 ((1.0_dp - kge(sm_catch_avg_domain, sm_opti_catch_avg_domain, mask = mask_times)) / &
716 real(domainmeta%overallnumberofdomains, dp))**6
717
718 ! deallocate
719 deallocate(mask_times)
720 deallocate(sm_catch_avg_domain)
721 deallocate(sm_opti_catch_avg_domain)
722 call smoptisim(idomain)%destroy()
723 end do
724 deallocate(smoptisim)
725
726#ifndef MPI
727 objective_sm_kge_catchment_avg = objective_sm_kge_catchment_avg**onesixth
728
729 call message(' objective_sm_kge_catchment_avg = ', num2str(objective_sm_kge_catchment_avg, '(F9.5)'))
730#endif
731
732
733 END FUNCTION objective_sm_kge_catchment_avg
734
735 ! ------------------------------------------------------------------
736
737 ! NAME
738 ! objective_q_et_tws_kge_catchment_avg
739
740 ! PURPOSE
741 !> \brief Objective function for et, tws and q.
742
743 !> \details The feature of this objective function is the
744 !> separation of the eval call into four
745 !> calls, each with another index list. The subroutine eval then only
746 !> uses the indices from that index list internally instead of having loops
747 !> over all domains. The integer array domainMeta%optidata decides which
748 !> indices to use. Therefore the array is split into disjunct subsets, and,
749 !> if chosen wisely in the namelist, also covers all domains.
750 !>
751 !> With this the eval calls sum up in a way that for each domain eval was
752 !> called at most once, but for different opti_data.
753
754 ! HISTORY
755 !> \authors Maren Kaluza
756
757 !> \date July 2019
758
759 ! Modifications:
760
761 FUNCTION objective_q_et_tws_kge_catchment_avg(parameterset, eval)
762
763 use mo_optimization_types, only : optidata_sim
764 use mo_common_constants, only : nodata_dp
765 use mo_common_variables, only : domainmeta
766 use mo_global_variables, only : l1_etobs, l1_twsaobs
767 use mo_errormeasures, only : kge
768 use mo_moment, only : average
769 use mo_string_utils, only : num2str
770 use mo_mrm_objective_function_runoff, only : extract_runoff
771
772 implicit none
773
774 !> the parameterset passed to the eval subroutine
775 real(dp), dimension(:), intent(in) :: parameterset
776 !> the eval subroutine called by this objective function
777 procedure(eval_interface), INTENT(IN), POINTER :: eval
778 !> the return value of the objective function. In this case it is
779 !> an array to provide the possibility to weight the outcome accordingly
780 real(dp), dimension(6) :: objective_q_et_tws_kge_catchment_avg
781
782 !> domain loop counter
783 integer(i4) :: idomain
784
785 !> counter for short loops
786 integer(i4) :: i
787#ifndef MPI
788 !> for sixth root
789 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
790#endif
791
792 !> modelled runoff for a given parameter set
793 ! dim1=nTimeSteps, dim2=nGauges
794 real(dp), allocatable, dimension(:, :) :: runoff
795 !> number of all gauges, aquired via runoff
796 integer(i4) :: ngaugestotal
797
798 !> aggregated simulated runoff
799 real(dp), dimension(:), allocatable :: runoff_agg
800
801 !> measured runoff
802 real(dp), dimension(:), allocatable :: runoff_obs
803
804 !> mask for measured runoff
805 logical, dimension(:), allocatable :: runoff_obs_mask
806
807 !> kge_q(nGaugesTotal)
808 real(dp) :: kge_q
809
810 !> gauges counter
811 integer(i4) :: gg, icell
812
813 !> number of q domains
814 integer(i4) :: nqdomains
815
816 !> number of et domains
817 integer(i4) :: netdomains
818
819 !> number of tws domains
820 integer(i4) :: ntwsdomains
821
822 !> number of TWS and ET domains (providing both)
823 integer(i4) :: nettwsdomains
824
825 !> index array of ET domains
826 integer(i4), dimension(:), allocatable :: opti_domain_indices_et
827
828 !> index array of TWS domains
829 integer(i4), dimension(:), allocatable :: opti_domain_indices_tws
830
831 !> index array of TWS and ET domains (providing both)
832 integer(i4), dimension(:), allocatable :: opti_domain_indices_et_tws
833
834 !> index array of ET domains
835 integer(i4), dimension(:), allocatable :: opti_domain_indices_q
836
837 !> simulated et
838 type(optidata_sim), dimension(:), allocatable :: etoptisim
839
840 !> simulated tws
841 type(optidata_sim), dimension(:), allocatable :: twsoptisim
842
843 !> simulated twsa (anomaly)
844 type(optidata_sim), dimension(:), allocatable :: twsaoptisim
845
846 real(dp) :: kge_tws
847
848 real(dp) :: kge_et
849
850 integer(i4) :: numberofsummands
851
852
853 ! initialize some variables
854 objective_q_et_tws_kge_catchment_avg(:) = 0.0_dp
855 kge_tws = 0.0_dp
856 kge_et = 0.0_dp
857 kge_q = 0.0_dp
858 numberofsummands = 0
859 !--------------------------------------------
860 ! ET & TWS
861 !--------------------------------------------
862 ! eval runs to get simulated output for et and tws
863 ! before each eval call we generate an index list of the domains for which
864 ! eval should be called. Read details for further information
865 call init_indexarray_for_opti_data(domainmeta, 6, nettwsdomains, opti_domain_indices_et_tws)
866 if (nettwsdomains > 0) then
867 allocate( etoptisim(domainmeta%nDomains))
868 allocate(twsoptisim(domainmeta%nDomains))
869 allocate(twsaoptisim(domainmeta%nDomains))
870 call eval(parameterset, opti_domain_indices = opti_domain_indices_et_tws, &
871 twsoptisim = twsoptisim, etoptisim = etoptisim)
872 ! for all domains that have ET and TWS
873 do i = 1, size(opti_domain_indices_et_tws)
874 idomain = opti_domain_indices_et_tws(i)
875 call convert_tws_to_twsa(twsoptisim(idomain), l1_twsaobs(idomain), twsaoptisim(idomain))
876 do icell = 1, size(l1_etobs(idomain)%maskObs(:, :), dim = 1)
877 kge_et = kge_et + &
878 (1.0_dp - kge(l1_etobs(idomain)%dataObs(icell, :), etoptisim(idomain)%dataSim(icell, :),&
879 mask = l1_etobs(idomain)%maskObs(icell, :)))**6
880 numberofsummands = numberofsummands + 1
881 end do
882 do icell = 1, size(l1_twsaobs(idomain)%maskObs(:, :), dim = 1)
883 kge_tws = kge_tws + &
884 (1.0_dp - kge(l1_twsaobs(idomain)%dataObs(icell, :), twsaoptisim(idomain)%dataSim(icell, :),&
885 mask = l1_twsaobs(idomain)%maskObs(icell, :)))**6
886 numberofsummands = numberofsummands + 1
887 end do
888 ! deallocate
889 call etoptisim(idomain)%destroy()
890 call twsoptisim(idomain)%destroy()
891 call twsaoptisim(idomain)%destroy()
892 end do
893 deallocate(etoptisim)
894 deallocate(twsoptisim)
895 deallocate(twsaoptisim)
896 ! write(0,*) 'nEtTwsDomains, kge_tws', nEtTwsDomains, kge_tws
897 ! write(0,*) 'nEtTwsDomains, kge_et', nEtTwsDomains, kge_et
898 end if
899 !--------------------------------------------
900 ! TWS
901 !--------------------------------------------
902 ! eval runs to get simulated output for tws
903 ! before each eval call we generate an index list of the domains for which
904 ! eval should be called. Read details for further information
905 call init_indexarray_for_opti_data(domainmeta, 3, ntwsdomains, opti_domain_indices_tws)
906 if (ntwsdomains > 0) then
907 allocate(twsoptisim(domainmeta%nDomains))
908 allocate(twsaoptisim(domainmeta%nDomains))
909 call eval(parameterset, opti_domain_indices = opti_domain_indices_tws, twsoptisim = twsoptisim)
910 ! for all domains that have ET and TWS
911 do i = 1, size(opti_domain_indices_tws)
912 idomain = opti_domain_indices_tws(i)
913 call convert_tws_to_twsa(twsoptisim(idomain), l1_twsaobs(idomain), twsaoptisim(idomain))
914 do icell = 1, size(l1_twsaobs(idomain)%maskObs(:, :), dim = 1)
915 kge_tws = kge_tws + &
916 (1.0_dp - kge(l1_twsaobs(idomain)%dataObs(icell, :), twsaoptisim(idomain)%dataSim(icell, :),&
917 mask = l1_twsaobs(idomain)%maskObs(icell, :)))**6
918 numberofsummands = numberofsummands + 1
919 end do
920 call twsoptisim(idomain)%destroy()
921 end do
922 deallocate(twsoptisim)
923 ! write(0,*) 'nTwsDomains, kge_tws', nTwsDomains, kge_tws
924 end if
925 objective_q_et_tws_kge_catchment_avg(2) = kge_tws
926
927 !--------------------------------------------
928 ! ET
929 !--------------------------------------------
930 ! eval runs to get simulated output for et
931 ! before each eval call we generate an index list of the domains for which
932 ! eval should be called. Read details for further information
933 call init_indexarray_for_opti_data(domainmeta, 5, netdomains, opti_domain_indices_et)
934 if (netdomains > 0) then
935 allocate(etoptisim(domainmeta%nDomains))
936 call eval(parameterset, opti_domain_indices = opti_domain_indices_et, etoptisim = etoptisim)
937 ! for all domains that have ET and TWS
938 do i = 1, size(opti_domain_indices_et)
939 idomain = opti_domain_indices_et(i)
940 do icell = 1, size(l1_etobs(idomain)%maskObs(:, :), dim = 1)
941 kge_et = kge_et + &
942 (1.0_dp - kge(l1_etobs(idomain)%dataObs(icell, :), etoptisim(idomain)%dataSim(icell, :),&
943 mask = l1_etobs(idomain)%maskObs(icell, :)))**6
944 numberofsummands = numberofsummands + 1
945 end do
946 call etoptisim(idomain)%destroy()
947 end do
948 deallocate(etoptisim)
949 ! write(0,*) 'nEtDomains, kge_et', nEtDomains, kge_et
950 end if
951 objective_q_et_tws_kge_catchment_avg(3) = kge_et
952
953 !--------------------------------------------
954 ! RUNOFF
955 !--------------------------------------------
956 ! eval runs to get simulated output for runoff
957 ! before the eval call we generate an index list of the domains for which
958 ! eval should be called. Read details for further information
959 ! ToDo: The arrays for qTin, qTout, will be rewritten in the other calls when
960 ! Q is not called last. Change that for more flexibility
961 call init_indexarray_for_opti_data(domainmeta, 1, nqdomains, opti_domain_indices_q)
962
963 if (nqdomains > 0) then
964 call eval(parameterset, opti_domain_indices = opti_domain_indices_q, runoff = runoff)
965 ngaugestotal = size(runoff, dim = 2)
966 do gg = 1, ngaugestotal
967
968 ! extract runoff
969 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
970
971 kge_q = kge_q + &
972 (1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask))**6
973 numberofsummands = numberofsummands + 1
974 deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
975 end do
976 ! write(0,*) 'nQDomains, kge_q', nQDomains, kge_q
977 end if
978 objective_q_et_tws_kge_catchment_avg(1) = kge_q
979
980 objective_q_et_tws_kge_catchment_avg(4) = real(numberofsummands, dp)
981
982
983#ifndef MPI
984 objective_q_et_tws_kge_catchment_avg(1) = ((kge_q+kge_et+kge_tws)/real(numberofsummands, dp))**onesixth
985
986 call message(' objective_q_et_tws_kge_catchment_avg = ', &
987 num2str(objective_q_et_tws_kge_catchment_avg(1), '(F9.5)'))
988#endif
989
990 END FUNCTION objective_q_et_tws_kge_catchment_avg
991
992 ! ------------------------------------------------------------------
993
994 ! NAME
995 ! init_indexarray_for_opti_data
996
997 ! PURPOSE
998 !> \brief creates an index array of the inidices of the domains eval
999 !> should MPI process.
1000 !
1001 !> \details The data type domainMeta contains an array optidata of size
1002 !> domainMeta%nDomains, telling us, which domains should be
1003 !> optimized with which opti_data. This subroutine splits all
1004 !> domains assigned to a process and returns an index list
1005 !> corresponding to the value of domainMeta%optidata.
1006 !>
1007 !> The index array opti_domain_indices can then be passed
1008 !> as an optional argument to the eval subroutine. The
1009 !> eval then instead of using loops over all domains only
1010 !> uses the passed indices.
1011 !>
1012 !> This subroutine also returns the size of that array since it
1013 !> helps with the calculations of the optimization in the end.
1014
1015 ! HISTORY
1016 !> \authors Maren Kaluza
1017
1018 !> \date July 2019
1019 subroutine init_indexarray_for_opti_data(domainMeta, optidataOption, nOptiDomains, opti_domain_indices)
1020 use mo_common_types, only: domain_meta
1021 !> meta data for all domains assigned to that process
1022 type(domain_meta), intent(in) :: domainMeta
1023 !> which opti data should be used in the eval called after calling this subroutine
1024 integer(i4), intent(in) :: optidataOption
1025 !> number of domains that will be optimized in the following eval call
1026 integer(i4), intent(out) :: nOptiDomains
1027 !> the indices of the domains that are to be processed in the following eval call
1028 integer(i4), dimension(:), allocatable, intent(out) :: opti_domain_indices
1029
1030 !> domain loop counter
1031 integer(i4) :: iDomain, i
1032
1033 if (allocated(opti_domain_indices)) deallocate(opti_domain_indices)
1034 ! count domains on MPI process that use optidata
1035 noptidomains = 0
1036 do idomain = 1, domainmeta%nDomains
1037 if (domainmeta%optidata(idomain) == optidataoption) noptidomains = noptidomains + 1
1038 end do
1039 ! write indices of these domains into an array
1040 if (noptidomains > 0) then
1041 allocate(opti_domain_indices(noptidomains))
1042 i = 0
1043 do idomain = 1, domainmeta%nDomains
1044 if (domainmeta%optidata(idomain) == optidataoption) then
1045 i = i + 1
1046 opti_domain_indices(i) = idomain
1047 end if
1048 end do
1049 end if
1050 end subroutine init_indexarray_for_opti_data
1051 ! ------------------------------------------------------------------
1052
1053 ! NAME
1054 ! objective_sm_corr
1055
1056 ! PURPOSE
1057 !> \brief Objective function for soil moisture.
1058
1059 !> \details The objective function only depends on a parameter vector.
1060 !> The model will be called with that parameter vector and
1061 !> the model output is subsequently compared to observed data.
1062
1063 !> Therefore the Pearson correlation between observed and modeled soil
1064 !> moisture on each grid cell \f$ j \f$ is compared
1065 !> \f[ r_j = r^2(SM_{obs}^j, SM_{sim}^j) \f]
1066 !> where
1067 !> \f$ r^2\f$ = Pearson correlation coefficient,
1068 !> \f$ SM_{obs} \f$ = observed soil moisture,
1069 !> \f$ SM_{sim} \f$ = simulated soil moisture.
1070 !> The observed data \f$ SM_{obs} \f$ are global in this module.
1071
1072 !> The the correlation is spatially averaged as
1073 !> \f[ \phi_{i} = \frac{1}{K} \cdot \sum_{j=1}^K r_j \f]
1074 !> where \f$ K \f$ denotes the number of valid cells in the study domain.
1075 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1076 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1077 !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 }. \f]
1078 !> The observed data L1_sm, L1_sm_mask are global in this module.
1079
1080 ! INTENT(IN)
1081 !> \param[in] "real(dp), dimension(:) :: parameterset"
1082 !> \param[in] "procedure(eval_interface) :: eval"
1083
1084 ! RETURN
1085 !> \return real(dp) :: objective_sm_corr — objective function value
1086 !> (which will be e.g. minimized by an optimization routine like DDS)
1087
1088 ! HISTORY
1089 !> \authors Matthias Zink
1090
1091 !> \date March 2015
1092
1093 ! Modifications:
1094 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1095
1096 FUNCTION objective_sm_corr(parameterset, eval)
1097
1098 use mo_optimization_types, only : optidata_sim
1099 use mo_common_variables, only : level1, domainmeta
1100 use mo_global_variables, only : l1_smobs
1101 use mo_moment, only : correlation
1102 use mo_string_utils, only : num2str
1103
1104 implicit none
1105
1106 real(dp), dimension(:), intent(in) :: parameterset
1107
1108 procedure(eval_interface), INTENT(IN), POINTER :: eval
1109
1110 real(dp) :: objective_sm_corr
1111
1112 ! domain loop counter
1113 integer(i4) :: idomain
1114
1115 ! cell loop counter
1116 integer(i4) :: icell
1117
1118 ! ncells1 of level 1
1119 integer(i4) :: ncells1
1120
1121 ! number of invalid cells in catchment
1122 real(dp) :: invalid_cells
1123
1124 ! domains wise objectives
1125 real(dp) :: objective_sm_corr_domain
1126
1127#ifndef MPI
1128 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1129#endif
1130
1131 type(optidata_sim), dimension(:), allocatable :: smoptisim
1132
1133
1134 allocate(smoptisim(domainmeta%nDomains))
1135 call eval(parameterset, smoptisim = smoptisim)
1136
1137 ! initialize some variables
1138 objective_sm_corr = 0.0_dp
1139
1140 ! loop over domain - for applying power law later on
1141 do idomain = 1, domainmeta%nDomains
1142
1143 ! init
1144 objective_sm_corr_domain = 0.0_dp
1145 ! get domain information
1146 ncells1 = level1(idomain)%ncells
1147
1148 invalid_cells = 0.0_dp
1149 ! temporal correlation is calculated on individual gridd cells
1150
1151 do icell = 1, size(l1_smobs(idomain)%maskObs(:, :), dim = 1)
1152
1153 ! check for enough data points in time for correlation
1154 if (count(l1_smobs(idomain)%maskObs(icell, :)) .LE. 0.10_dp * real(size(l1_smobs(idomain)%dataObs(:, :), dim = 2), dp)) then
1155 invalid_cells = invalid_cells + 1.0_dp
1156 cycle
1157 end if
1158 objective_sm_corr_domain = objective_sm_corr_domain + &
1159 correlation(l1_smobs(idomain)%dataObs(icell, :), smoptisim(idomain)%dataSim(icell, :), &
1160 mask = l1_smobs(idomain)%maskObs(icell, :))
1161 end do
1162
1163 ! user information about invalid cells
1164 if (invalid_cells .GT. 0.5_dp) then
1165 call message(.LT.' WARNING: objective_sm: Detected invalid cells in study area ( 10 valid data points).')
1166 call message(' Fraction of invalid cells: ', &
1167 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
1168 end if
1169
1170
1171 ! calculate average soil moisture correlation over all domains with power law
1172 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1173 objective_sm_corr = objective_sm_corr + &
1174 ((1.0_dp - objective_sm_corr_domain / real(ncells1, dp)) / real(domainmeta%overallNumberOfDomains, dp))**6
1175 end do
1176#ifndef MPI
1177 objective_sm_corr = objective_sm_corr**onesixth
1178
1179 call message(' objective_sm_corr = ', num2str(objective_sm_corr, '(F9.5)'))
1180#endif
1181
1182 END FUNCTION objective_sm_corr
1183
1184 ! ------------------------------------------------------------------
1185
1186 ! NAME
1187 ! objective_sm_pd
1188
1189 ! PURPOSE
1190 !> \brief Objective function for soil moisture.
1191
1192 !> \details The objective function only depends on a parameter vector.
1193 !> The model will be called with that parameter vector and
1194 !> the model output is subsequently compared to observed data.
1195
1196 !> Therefore the Pattern Dissimilarity (PD) of observed and modeled soil
1197 !> moisture fields is calculated - aim: matching spatial patters
1198 !> \f[ E(t) = PD\left( SM_{obs}(t), SM_{sim}(t) \right) \f]
1199 !> where
1200 !> \f$ PD \f$ = pattern dissimilarity function,
1201 !> \f$ SM_{obs} \f$ = observed soil moisture,
1202 !> \f$ SM_{sim} \f$ = simulated soil moisture.
1203 !> \f$ E(t) \f$ = pattern dissimilarity at timestep \f$ t \f$.
1204 !> The the pattern dissimilaity (E) is spatially averaged as
1205 !> \f[ \phi_{i} = \frac{1}{T} \cdot \sum_{t=1}^T E_t \f]
1206 !> where \f$ T \f$ denotes the number of time steps.
1207 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1208 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1209 !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 } . \f]
1210 !> The observed data L1_sm, L1_sm_mask are global in this module.
1211 !> The observed data L1_sm, L1_sm_mask are global in this module.
1212
1213 ! INTENT(IN)
1214 !> \param[in] "real(dp), dimension(:) :: parameterset"
1215 !> \param[in] "procedure(eval_interface) :: eval"
1216
1217 ! RETURN
1218 !> \return real(dp) :: objecive_sm_pd — objective function value
1219 !> (which will be e.g. minimized by an optimization routine like DDS)
1220
1221 ! HISTORY
1222 !> \authors Matthias Zink
1223
1224 !> \date May 2015
1225
1226 ! Modifications:
1227 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1228
1229 FUNCTION objective_sm_pd(parameterset, eval)
1230
1231 use mo_optimization_types, only : optidata_sim
1232 use mo_common_constants, only : nodata_dp
1233 use mo_common_variables, only : level1, domainmeta
1234 use mo_global_variables, only : l1_smobs
1235 use mo_spatialsimilarity, only : pd
1236 use mo_string_utils, only : num2str
1237
1238 implicit none
1239
1240 real(dp), dimension(:), intent(in) :: parameterset
1241
1242 procedure(eval_interface), INTENT(IN), POINTER :: eval
1243
1244 ! objective function value
1245 real(dp) :: objective_sm_pd
1246
1247 ! domain loop counter
1248 integer(i4) :: idomain
1249
1250 ! time loop counter
1251 integer(i4) :: itime
1252
1253 ! level 1 number of culomns and rows
1254 integer(i4) :: nrows1, ncols1
1255
1256 ! for sixth root
1257#ifndef MPI
1258 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1259#endif
1260
1261 ! matrices of SM from vectorized arrays
1262 real(dp), dimension(:, :), allocatable :: mat1, mat2
1263
1264 ! pattern dissimilarity (pd) at every time step
1265 real(dp), dimension(:), allocatable :: pd_time_series
1266
1267 ! simulated soil moisture
1268 type(optidata_sim), dimension(:), allocatable :: smoptisim
1269
1270 ! mask of valid cells at level1
1271 logical, dimension(:, :), allocatable :: mask1
1272
1273 ! mask of valid sm cells
1274 logical, dimension(:, :), allocatable :: mask_sm
1275
1276 ! mask for valid sm catchment avg time steps
1277 logical, dimension(:), allocatable :: mask_times
1278
1279
1280 allocate(smoptisim(domainmeta%nDomains))
1281 call eval(parameterset, smoptisim = smoptisim)
1282
1283 ! initialize some variables
1284 objective_sm_pd = 0.0_dp
1285
1286 ! loop over domain - for applying power law later on
1287 do idomain = 1, domainmeta%nDomains
1288
1289 ! get domain information
1290 mask1 = level1(idomain)%mask
1291 ncols1 = level1(idomain)%ncols
1292 nrows1 = level1(idomain)%nrows
1293
1294 ! allocate
1295 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
1296 allocate(pd_time_series(size(smoptisim(idomain)%dataSim, dim = 2)))
1297 allocate(mat1(nrows1, ncols1))
1298 allocate(mat2(nrows1, ncols1))
1299 allocate(mask_sm(nrows1, ncols1))
1300
1301 ! initalize
1302 mask_times = .false.
1303 pd_time_series = 0.0_dp
1304
1305 ! calculate pattern similarity criterion
1306 do itime = 1, size(smoptisim(idomain)%dataSim, dim = 2)
1307 mat1 = unpack(l1_smobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
1308 mat2 = unpack(smoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
1309 mask_sm = unpack(l1_smobs(idomain)%maskObs(:, itime), mask1, .false.)
1310 pd_time_series = pd(mat1, mat2, mask = mask_sm, valid = mask_times(itime))
1311 end do
1312
1313 if (count(mask_times) > 0_i4) then
1314 ! calculate avergae PD over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1315 objective_sm_pd = objective_sm_pd + &
1316 ((1.0_dp - sum(pd_time_series, mask = mask_times) / real(count(mask_times), dp)) / &
1317 real(domainmeta%overallnumberofdomains, dp))**6
1318 else
1319 call error_message('***ERROR: mo_objective_funtion: objective_sm_pd: No soil moisture observations available!')
1320 end if
1321
1322 ! deallocate
1323 deallocate(mask_times)
1324 deallocate(pd_time_series)
1325 deallocate(mat1)
1326 deallocate(mat2)
1327 deallocate(mask_sm)
1328 call smoptisim(idomain)%destroy()
1329 end do
1330 deallocate(smoptisim)
1331
1332#ifndef MPI
1333 objective_sm_pd = objective_sm_pd**onesixth
1334
1335 call message(' objective_sm_pd = ', num2str(objective_sm_pd, '(F9.5)'))
1336#endif
1337
1338 END FUNCTION objective_sm_pd
1339
1340 ! ------------------------------------------------------------------
1341
1342 ! NAME
1343 ! objective_sm_sse_standard_score
1344
1345 ! PURPOSE
1346 !> \brief Objective function for soil moisture.
1347
1348 !> \details The objective function only depends on a parameter vector.
1349 !> The model will be called with that parameter vector and
1350 !> the model output is subsequently compared to observed data.
1351
1352 !> Therefore the sum of squared errors (SSE) of the standard score of observed and
1353 !> modeled soil moisture is calculated. The standard score or normalization (anomaly)
1354 !> make the objctive function bias insensitive and basically the dynamics of the soil moisture
1355 !> is tried to capture by this obejective function.
1356 !> \f[ phi_i = \sum_{j=1}^K \{ standard\_score( SM_{obs}(j) )- standard\_score(SM_{sim}(j)) \}^2 \f]
1357 !> where
1358 !> \f$ standard\_score \f$ = standard score function,
1359 !> \f$ SM_{obs} \f$ = observed soil moisture,
1360 !> \f$ SM_{sim} \f$ = simulated soil moisture.
1361 !> \f$ K \f$ = valid elements in study domain.
1362 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1363 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1364 !> \f[ OF = \sqrt[6]{\sum(\phi_{i}/N)^6 }. \f]
1365 !> The observed data L1_sm, L1_sm_mask are global in this module.
1366
1367 ! INTENT(IN)
1368 !> \param[in] "real(dp), dimension(:) :: parameterset"
1369 !> \param[in] "procedure(eval_interface) :: eval"
1370
1371 ! RETURN
1372 !> \return real(dp) :: objective_sm_sse_standard_score — objective function value
1373 !> (which will be e.g. minimized by an optimization routine like DDS)
1374
1375 ! HISTORY
1376 !> \authors Matthias Zink
1377
1378 !> \date March 2015
1379
1380 ! Modifications:
1381 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1382
1383 FUNCTION objective_sm_sse_standard_score(parameterset, eval)
1384
1385 use mo_optimization_types, only : optidata_sim
1386 use mo_common_variables, only : level1, domainmeta
1387 use mo_errormeasures, only : sse
1388 use mo_global_variables, only : l1_smobs
1389 use mo_standard_score, only : standard_score
1390 use mo_string_utils, only : num2str
1391
1392 implicit none
1393
1394 real(dp), dimension(:), intent(in) :: parameterset
1395
1396 procedure(eval_interface), INTENT(IN), POINTER :: eval
1397
1398 real(dp) :: objective_sm_sse_standard_score
1399
1400 ! domain loop counter
1401 integer(i4) :: idomain
1402
1403 ! cell loop counter
1404 integer(i4) :: icell
1405
1406 ! ncells1 of level 1
1407 integer(i4) :: ncells1
1408
1409 ! number of invalid cells in catchment
1410 real(dp) :: invalid_cells
1411
1412 ! domains wise objectives
1413 real(dp) :: objective_sm_sse_standard_score_domain
1414
1415#ifndef MPI
1416 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1417#endif
1418
1419 ! simulated soil moisture
1420 type(optidata_sim), dimension(:), allocatable :: smoptisim
1421
1422
1423 call eval(parameterset, smoptisim = smoptisim)
1424
1425 ! initialize some variables
1426 objective_sm_sse_standard_score = 0.0_dp
1427
1428 ! loop over domain - for applying power law later on
1429 do idomain = 1, domainmeta%nDomains
1430
1431 ! init
1432 objective_sm_sse_standard_score_domain = 0.0_dp
1433 ! get domain information
1434 ncells1 = level1(idomain)%nCells
1435
1436 invalid_cells = 0.0_dp
1437 ! standard_score signal is calculated on individual grid cells
1438 do icell = 1, size(l1_smobs(idomain)%maskObs(:, :), dim = 1)
1439
1440 ! check for enough data points in time for statistical calculations (e.g. mean, stddev)
1441 if (count(l1_smobs(idomain)%maskObs(icell, :)) .LE. (0.10_dp * real(size(l1_smobs(idomain)%dataObs, dim = 2), dp))) then
1442 invalid_cells = invalid_cells + 1.0_dp
1443 cycle
1444 end if
1445 objective_sm_sse_standard_score_domain = objective_sm_sse_standard_score_domain + &
1446 sse(standard_score(l1_smobs(idomain)%dataObs(icell, :), mask = l1_smobs(idomain)%maskObs(icell, :)), &
1447 standard_score(smoptisim(idomain)%dataSim(icell, :), mask = l1_smobs(idomain)%maskObs(icell, :)), &
1448 mask = l1_smobs(idomain)%maskObs(icell, :))
1449
1450 end do
1451
1452 ! user information about invalid cells
1453 if (invalid_cells .GT. 0.5_dp) then
1454 call message(.LT.' WARNING: objective_sm: Detected invalid cells in study area ( 10 valid data points).')
1455 call message(' Fraction of invalid cells: ', &
1456 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
1457 end if
1458
1459 ! calculate average soil moisture correlation over all domains with power law
1460 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1461 objective_sm_sse_standard_score = objective_sm_sse_standard_score + &
1462 (objective_sm_sse_standard_score_domain / real(domainmeta%overallNumberOfDomains, dp))**6
1463 end do
1464
1465#ifndef MPI
1466 objective_sm_sse_standard_score = objective_sm_sse_standard_score**onesixth
1467
1468 call message(' objective_sm_sse_standard_score = ', num2str(objective_sm_sse_standard_score, '(E12.5)'))
1469#endif
1470
1471 END FUNCTION objective_sm_sse_standard_score
1472
1473
1474 ! -----------------------------------------------------------------
1475
1476 ! NAME
1477 ! objective_kge_q_rmse_tws
1478
1479 ! PURPOSE
1480 !> \brief Objective function of KGE for runoff and RMSE for domain_avg TWS (standarized scores)
1481
1482 !> \details Objective function of KGE for runoff and RMSE for domain_avg TWS (standarized scores)
1483
1484 ! INTENT(IN)
1485 !> \param[in] "real(dp), dimension(:) :: parameterset"
1486 !> \param[in] "procedure(eval_interface) :: eval"
1487
1488 ! RETURN
1489 !> \return real(dp) :: objective_kge_q_rmse_tws — objective function value
1490 !> (which will be e.g. minimized by an optimization routine like DDS)
1491
1492 ! HISTORY
1493 !> \authors Oldrich Rakovec, Rohini Kumar
1494
1495 !> \date Oct. 2015
1496
1497 ! Modifications:
1498 ! Stephan Thober Oct 2015 - moved tws optimization from mo_mrm_objective_function_runoff here
1499 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1500 ! Maren Kaluza Oct 2019 - changed averaging function for tws, this will not produce the same output as before
1501
1502 FUNCTION objective_kge_q_rmse_tws(parameterset, eval)
1503
1504 use mo_optimization_types, only : optidata_sim
1505 use mo_common_constants, only : eps_dp, nodata_dp
1506 use mo_common_mhm_mrm_variables, only : evalper
1507 use mo_common_variables, only : domainmeta
1508 use mo_global_variables, only : l1_twsaobs
1509 use mo_errormeasures, only : rmse
1510 use mo_julian, only : caldat
1511 use mo_moment, only : mean
1512 use mo_standard_score, only : classified_standard_score
1513 use mo_string_utils, only : num2str
1514 use mo_temporal_aggregation, only : day2mon_average
1515 use mo_errormeasures, only : kge
1516 use mo_mrm_objective_function_runoff, only : extract_runoff
1517
1518 implicit none
1519
1520 real(dp), dimension(:), intent(in) :: parameterset
1521
1522 procedure(eval_interface), INTENT(IN), POINTER :: eval
1523
1524 real(dp) :: objective_kge_q_rmse_tws
1525
1526 ! modelled runoff for a given parameter set
1527 ! dim1=nTimeSteps, dim2=nGauges
1528 real(dp), allocatable, dimension(:, :) :: runoff
1529
1530 !> simulated tws
1531 type(optidata_sim), dimension(:), allocatable :: twsoptisim
1532
1533 ! domain counter, month counters
1534 integer(i4) :: domainid, idomain, pp, mmm
1535
1536 integer(i4) :: year, month, day
1537
1538 real(dp), dimension(domainMeta%nDomains) :: inittime
1539
1540 ! simulated tws
1541 real(dp), dimension(:), allocatable :: tws_catch_avg_domain
1542
1543 ! measured tws
1544 real(dp), dimension(:), allocatable :: tws_opti_catch_avg_domain
1545
1546 ! mask for measured tws
1547 logical, dimension(:), allocatable :: tws_obs_mask
1548
1549 ! total number of months
1550 integer(i4) :: nmonths
1551
1552 ! vector with months' classes
1553 integer(i4), dimension(:), allocatable :: month_classes
1554
1555 ! monthly values anomaly time series
1556 real(dp), DIMENSION(:), allocatable :: tws_sim_m_anom, tws_obs_m_anom
1557
1558 ! rmse_tws(domainMeta%nDomains)
1559 real(dp), dimension(:), allocatable :: rmse_tws
1560
1561 ! obj. functions
1562 real(dp) :: rmse_tws_avg, kge_q_avg
1563
1564 integer(i4) :: ngaugestotal
1565
1566 ! aggregated simulated runoff
1567 real(dp), dimension(:), allocatable :: runoff_agg
1568
1569 ! measured runoff
1570 real(dp), dimension(:), allocatable :: runoff_obs
1571
1572 ! mask for measured runoff
1573 logical, dimension(:), allocatable :: runoff_obs_mask
1574
1575 ! kge_q(nGaugesTotal)
1576 real(dp), dimension(:), allocatable :: kge_q
1577
1578 ! gauges counter
1579 integer(i4) :: gg
1580
1581 ! obtain hourly values of runoff and tws:
1582 allocate(twsoptisim(domainmeta%nDomains))
1583 call eval(parameterset, runoff = runoff, twsoptisim = twsoptisim)
1584
1585 !--------------------------------------------
1586 !! TWS
1587 !--------------------------------------------
1588
1589 ! allocate
1590 allocate(rmse_tws(domainmeta%nDomains))
1591 rmse_tws(:) = nodata_dp
1592
1593 do idomain = 1, domainmeta%nDomains
1594 if (.not. (l1_twsaobs(idomain)%timeStepInput == -2)) then
1595 call message('objective_kge_q_rmse_tws: current implementation of this subroutine only allows monthly timesteps')
1596 end if
1597 domainid = domainmeta%indices(idomain)
1598
1599 ! extract tws the same way as runoff using mrm
1600 ! Note that with the change from tws(iDomain, tt) to tws(tt, :) this
1601 ! will not work like before and also does maybe not make sense
1602 call create_domain_avg_tws(idomain, twsoptisim, tws_catch_avg_domain, tws_opti_catch_avg_domain, tws_obs_mask)
1603
1604 ! check for potentially 2 years of data
1605 if (count(tws_obs_mask) .lt. 12 * 2) then
1606 call message('objective_kge_q_rmse_tws: Length of TWS data of domain ', trim(adjustl(num2str(domainid))), &
1607 ' less than 2 years: this is not recommended')
1608 end if
1609
1610 ! get initial time of the evaluation period
1611 inittime(idomain) = real(evalper(idomain)%julStart, dp)
1612
1613 ! get calendar days, months, year
1614 call caldat(int(inittime(idomain)), yy = year, mm = month, dd = day)
1615
1616 nmonths = size(tws_obs_mask)
1617
1618 allocate (month_classes(nmonths))
1619 allocate (tws_obs_m_anom(nmonths))
1620 allocate (tws_sim_m_anom(nmonths))
1621
1622 month_classes(:) = 0
1623 tws_obs_m_anom(:) = nodata_dp
1624 tws_sim_m_anom(:) = nodata_dp
1625
1626 ! define months' classes
1627 mmm = month
1628 do pp = 1, nmonths
1629 month_classes(pp) = mmm
1630 if (mmm .LT. 12) then
1631 mmm = mmm + 1
1632 else
1633 mmm = 1
1634 end if
1635 end do
1636
1637 ! calculate standard score
1638 tws_obs_m_anom = classified_standard_score(tws_opti_catch_avg_domain, month_classes, mask = tws_obs_mask)
1639 tws_sim_m_anom = classified_standard_score(tws_catch_avg_domain, month_classes, mask = tws_obs_mask)
1640 rmse_tws(idomain) = rmse(tws_sim_m_anom, tws_obs_m_anom, mask = tws_obs_mask)
1641
1642 deallocate (month_classes)
1643 deallocate (tws_sim_m_anom)
1644 deallocate (tws_obs_m_anom)
1645 deallocate (tws_catch_avg_domain, tws_opti_catch_avg_domain, tws_obs_mask)
1646
1647 call twsoptisim(idomain)%destroy()
1648 end do
1649
1650 rmse_tws_avg = sum(rmse_tws(:), abs(rmse_tws - nodata_dp) .gt. eps_dp) / &
1651 real(count(abs(rmse_tws - nodata_dp) .gt. eps_dp), dp)
1652 deallocate(rmse_tws)
1653 deallocate(twsoptisim)
1654
1655 !--------------------------------------------
1656 !! RUNOFF
1657 !--------------------------------------------
1658 kge_q_avg = 0_dp
1659 ngaugestotal = size(runoff, dim = 2)
1660 allocate(kge_q(ngaugestotal))
1661 kge_q(:) = nodata_dp
1662
1663 do gg = 1, ngaugestotal
1664
1665 ! extract runoff
1666 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
1667
1668 ! check for potentially 2 years of data
1669 pp = count(runoff_agg .ge. 0.0_dp)
1670 if (pp .lt. 365 * 2) then
1671 call message('objective_kge_q_rmse_tws: The simulation at gauge ', trim(adjustl(num2str(gg))), &
1672 ' is not long enough. Please provide at least 730 days of data.')
1673 end if
1674 ! calculate KGE for each domain:
1675 kge_q(gg) = kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
1676 deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
1677
1678 end do
1679
1680 ! calculate average KGE value for runoff
1681 kge_q_avg = sum(kge_q(:), abs(kge_q - nodata_dp) .gt. eps_dp) / &
1682 real(count(abs(kge_q - nodata_dp) .gt. eps_dp), dp)
1683 deallocate(kge_q)
1684
1685 !
1686 objective_kge_q_rmse_tws = rmse_tws_avg * (1._dp - kge_q_avg)
1687
1688 call message(' objective_kge_q_rmse_tws = ', num2str(objective_kge_q_rmse_tws, '(F9.5)'))
1689
1690 END FUNCTION objective_kge_q_rmse_tws
1691
1692 ! ------------------------------------------------------------------
1693
1694 ! NAME
1695 ! objective_neutrons_kge_catchment_avg
1696
1697 ! PURPOSE
1698 !> \brief Objective function for neutrons.
1699
1700 !> \details The objective function only depends on a parameter vector.
1701 !> The model will be called with that parameter vector and
1702 !> the model output is subsequently compared to observed data.
1703
1704 !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average
1705 !> neutrons (N) is calculated
1706 !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
1707 !> where
1708 !> \f$ r \f$ = Pearson product-moment CORRELATION coefficient
1709 !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM
1710 !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
1711 !> is calculated and the objective function for a given domain \f$ i \f$ is
1712 !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
1713 !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
1714 !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
1715
1716 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1717 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1718 !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f]
1719 !> The observed data L1_neutronsdata, L1_neutronsdata_mask are global in this module.
1720
1721 ! INTENT(IN)
1722 !> \param[in] "real(dp), dimension(:) :: parameterset"
1723 !> \param[in] "procedure(eval_interface) :: eval"
1724
1725 ! RETURN
1726 !> \return real(dp) :: objective_neutrons_kge_catchment_avg — objective function value
1727 !> (which will be e.g. minimized by an optimization routine)
1728
1729 ! HISTORY
1730 !> \authors Martin Schroen
1731
1732 !> \date Jun 2015
1733
1734 ! Modifications:
1735 ! Maren Kaluza Mar 2018 - changed format string to '(I10)'
1736 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1737
1738 FUNCTION objective_neutrons_kge_catchment_avg(parameterset, eval)
1739
1740 use mo_optimization_types, only : optidata_sim
1741 use mo_common_constants, only : nodata_dp
1742 use mo_common_variables, only : domainmeta
1743 use mo_errormeasures, only : kge
1744 use mo_global_variables, only : l1_neutronsobs
1745 use mo_moment, only : average
1746 use mo_string_utils, only : num2str
1747
1748 implicit none
1749
1750 real(dp), dimension(:), intent(in) :: parameterset
1751
1752 procedure(eval_interface), INTENT(IN), POINTER :: eval
1753
1754 real(dp) :: objective_neutrons_kge_catchment_avg
1755
1756 ! domain loop counter
1757 integer(i4) :: idomain
1758
1759 ! time loop counter
1760 integer(i4) :: itime
1761
1762 ! for sixth root
1763#ifndef MPI
1764 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1765#endif
1766
1767 ! spatial average of observed neutrons
1768 real(dp), dimension(:), allocatable :: neutrons_catch_avg_domain
1769
1770 ! spatial avergae of modeled neutrons
1771 real(dp), dimension(:), allocatable :: neutrons_opti_catch_avg_domain
1772
1773 ! simulated neutrons
1774 type(optidata_sim), dimension(:), allocatable :: neutronsoptisim
1775
1776 ! mask for valid neutrons catchment avg time steps
1777 logical, dimension(:), allocatable :: mask_times
1778
1779
1780 allocate(neutronsoptisim(domainmeta%nDomains))
1781 call eval(parameterset, neutronsoptisim = neutronsoptisim)
1782
1783 ! initialize some variables
1784 objective_neutrons_kge_catchment_avg = 0.0_dp
1785
1786 ! loop over domain - for applying power law later on
1787 do idomain = 1, domainmeta%nDomains
1788
1789 ! allocate
1790 allocate(mask_times(size(neutronsoptisim(idomain)%dataSim, dim = 2)))
1791 allocate(neutrons_catch_avg_domain(size(neutronsoptisim(idomain)%dataSim, dim = 2)))
1792 allocate(neutrons_opti_catch_avg_domain(size(neutronsoptisim(idomain)%dataSim, dim = 2)))
1793
1794 ! initalize
1795 mask_times = .true.
1796 neutrons_catch_avg_domain = nodata_dp
1797 neutrons_opti_catch_avg_domain = nodata_dp
1798
1799 ! calculate catchment average soil moisture
1800 do itime = 1, size(neutronsoptisim(idomain)%dataSim, dim = 2)
1801
1802 ! check for enough data points in time for correlation
1803 if (all(.NOT. l1_neutronsobs(idomain)%maskObs(:, itime))) then
1804 call message('WARNING: neutrons data at time ', num2str(itime, '(I10)'), ' is empty.')
1805 !call message('WARNING: objective_neutrons_kge_catchment_avg: ignored current time step since less than')
1806 !call message(' 10 valid cells available in soil moisture observation')
1807 mask_times(itime) = .false.
1808 cycle
1809 end if
1810 neutrons_catch_avg_domain(itime) = average(l1_neutronsobs(idomain)%dataObs(:, itime), &
1811 mask = l1_neutronsobs(idomain)%maskObs(:, itime))
1812 neutrons_opti_catch_avg_domain(itime) = average(neutronsoptisim(idomain)%dataSim(:, itime), &
1813 mask = l1_neutronsobs(idomain)%maskObs(:, itime))
1814 end do
1815
1816 ! calculate average neutrons KGE over all domains with power law
1817 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1818 objective_neutrons_kge_catchment_avg = objective_neutrons_kge_catchment_avg + &
1819 ((1.0_dp - kge(neutrons_catch_avg_domain, neutrons_opti_catch_avg_domain, mask = mask_times)) / &
1820 real(domainmeta%overallnumberofdomains, dp))**6
1821
1822 ! deallocate
1823 deallocate(mask_times)
1824 deallocate(neutrons_catch_avg_domain)
1825 deallocate(neutrons_opti_catch_avg_domain)
1826
1827 call neutronsoptisim(idomain)%destroy()
1828 end do
1829 deallocate(neutronsoptisim)
1830
1831#ifndef MPI
1832 objective_neutrons_kge_catchment_avg = objective_neutrons_kge_catchment_avg**onesixth
1833
1834 call message(' objective_neutrons_kge_catchment_avg = ', num2str(objective_neutrons_kge_catchment_avg, '(F9.5)'))
1835#endif
1836
1837 END FUNCTION objective_neutrons_kge_catchment_avg
1838
1839 ! ------------------------------------------------------------------
1840
1841 ! NAME
1842 ! objective_et_kge_catchment_avg
1843
1844 ! PURPOSE
1845 !> \brief Objective function for evpotranspirstion (et).
1846
1847 !> \details The objective function only depends on a parameter vector.
1848 !> The model will be called with that parameter vector and
1849 !> the model output is subsequently compared to observed data.
1850
1851 !> Therefore, the Kling-Gupta model efficiency \f$ KGE \f$ of the catchment average
1852 !> evapotranspiration (et) is calculated
1853 !> \f[ KGE = 1.0 - \sqrt{( (1-r)^2 + (1-\alpha)^2 + (1-\beta)^2 )} \f]
1854 !> where
1855 !> \f$ r \f$ = Pearson product-moment correlation coefficient
1856 !> \f$ \alpha \f$ = ratio of simulated mean to observed mean SM
1857 !> \f$ \beta \f$ = ratio of similated standard deviation to observed standard deviation
1858 !> is calculated and the objective function for a given domain \f$ i \f$ is
1859 !> \f[ \phi_{i} = 1.0 - KGE_{i} \f]
1860 !> \f$ \phi_{i} \f$ is the objective since we always apply minimization methods.
1861 !> The minimal value of \f$ \phi_{i} \f$ is 0 for the optimal KGE of 1.0.
1862
1863 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
1864 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
1865 !> \f[ OF = \sqrt[6]{\sum((1.0 - KGE_{i})/N)^6 }. \f]
1866 !> The observed data L1_et, L1_et_mask are global in this module.
1867
1868 ! INTENT(IN)
1869 !> \param[in] "real(dp), dimension(:) :: parameterset"
1870 !> \param[in] "procedure(eval_interface) :: eval"
1871
1872 ! RETURN
1873 !> \return real(dp) :: objective_et_kge_catchment_avg — objective function value
1874 !> (which will be e.g. minimized by an optimization routine)
1875
1876 ! HISTORY
1877 !> \authors Johannes Brenner
1878
1879 !> \date Feb 2017
1880
1881 ! Modifications:
1882 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1883
1884 FUNCTION objective_et_kge_catchment_avg(parameterset, eval)
1885
1886 use mo_optimization_types, only : optidata_sim
1887 use mo_common_constants, only : nodata_dp
1888 use mo_common_variables, only : domainmeta
1889 use mo_errormeasures, only : kge
1890 use mo_moment, only : average
1891 use mo_string_utils, only : num2str
1892
1893 implicit none
1894
1895 real(dp), dimension(:), intent(in) :: parameterset
1896
1897 procedure(eval_interface), INTENT(IN), POINTER :: eval
1898
1899 real(dp) :: objective_et_kge_catchment_avg
1900
1901 ! domain loop counter
1902 integer(i4) ::idomain
1903
1904 ! for sixth root
1905#ifndef MPI
1906 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
1907#endif
1908
1909 ! spatial average of observed et
1910 real(dp), dimension(:), allocatable :: et_catch_avg_domain
1911
1912 ! spatial avergae of modeled et
1913 real(dp), dimension(:), allocatable :: et_opti_catch_avg_domain
1914
1915 !> simulated et
1916 type(optidata_sim), dimension(:), allocatable :: etoptisim
1917
1918 ! mask for valid et catchment avg time steps
1919 logical, dimension(:), allocatable :: mask_times
1920
1921
1922 call eval(parameterset, etoptisim = etoptisim)
1923
1924 ! initialize some variables
1925 objective_et_kge_catchment_avg = 0.0_dp
1926
1927 ! loop over domain - for applying power law later on
1928 allocate(etoptisim(domainmeta%nDomains))
1929 do idomain = 1, domainmeta%nDomains
1930
1931 ! create et array input
1932 call create_domain_avg_et(idomain, etoptisim, et_catch_avg_domain, &
1933 et_opti_catch_avg_domain, mask_times)
1934 ! calculate average ET KGE over all domains with power law
1935 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
1936
1937 objective_et_kge_catchment_avg = objective_et_kge_catchment_avg + &
1938 ((1.0_dp - kge(et_catch_avg_domain, et_opti_catch_avg_domain, mask = mask_times)) / &
1939 real(domainmeta%overallnumberofdomains, dp))**6
1940
1941 ! deallocate
1942 deallocate(mask_times)
1943 deallocate(et_catch_avg_domain)
1944 deallocate(et_opti_catch_avg_domain)
1945 call etoptisim(idomain)%destroy()
1946 end do
1947 deallocate(etoptisim)
1948
1949#ifndef MPI
1950 objective_et_kge_catchment_avg = objective_et_kge_catchment_avg**onesixth
1951
1952 call message(' objective_et_kge_catchment_avg = ', num2str(objective_et_kge_catchment_avg, '(F9.5)'))
1953#endif
1954
1955 END FUNCTION objective_et_kge_catchment_avg
1956
1957 ! -----------------------------------------------------------------
1958
1959 ! NAME
1960 ! objective_kge_q_sm_corr
1961
1962 ! PURPOSE
1963 !> \brief Objective function of KGE for runoff and correlation for SM
1964
1965 !> \details Objective function of KGE for runoff and SSE for soil moisture (standarized scores).
1966 !> Further details can be found in the documentation of objective functions
1967 !> '14 - objective_multiple_gauges_kge_power6' and '13 - objective_sm_corr'.
1968
1969 ! INTENT(IN)
1970 !> \param[in] "real(dp), dimension(:) :: parameterset"
1971 !> \param[in] "procedure(eval_interface) :: eval"
1972
1973 ! RETURN
1974 !> \return real(dp) :: objective_kge_q_sse_sm — objective function value
1975 !> (which will be e.g. minimized by an optimization routine like DDS)
1976
1977 ! HISTORY
1978 !> \authors Matthias Zink
1979
1980 !> \date Mar. 2017
1981
1982 ! Modifications:
1983 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1984
1985 FUNCTION objective_kge_q_sm_corr(parameterset, eval)
1986
1987 use mo_optimization_types, only : optidata_sim
1988 use mo_common_variables, only : level1, domainmeta
1989 use mo_global_variables, only : l1_smobs
1990 use mo_moment, only : correlation
1991 use mo_string_utils, only : num2str
1992 use mo_errormeasures, only : kge
1993 use mo_mrm_objective_function_runoff, only : extract_runoff
1994
1995 implicit none
1996
1997 real(dp), dimension(:), intent(in) :: parameterset
1998
1999 procedure(eval_interface), INTENT(IN), POINTER :: eval
2000
2001 real(dp) :: objective_kge_q_sm_corr
2002
2003 real(dp) :: objective_sm
2004
2005 real(dp) :: objective_kge
2006
2007 ! number of invalid cells in catchment
2008 real(dp) :: invalid_cells
2009
2010 ! modelled runoff for a given parameter set
2011 ! dim1=nTimeSteps, dim2=nGauges
2012 real(dp), allocatable, dimension(:, :) :: runoff
2013
2014 ! domain loop counter
2015 integer(i4) :: idomain
2016
2017 ! cell loop counter
2018 integer(i4) :: icell
2019
2020 ! ncells1 of level 1
2021 integer(i4) :: ncells1
2022
2023 ! domains wise objectives
2024 real(dp) :: objective_sm_domain
2025
2026 ! simulated soil moisture
2027 type(optidata_sim), dimension(:), allocatable :: smoptisim
2028
2029 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2030
2031 ! gauges counter
2032 integer(i4) :: gg
2033
2034 integer(i4) :: ngaugestotal
2035
2036 ! aggregated simulated runoff
2037 real(dp), dimension(:), allocatable :: runoff_agg
2038
2039 ! measured runoff
2040 real(dp), dimension(:), allocatable :: runoff_obs
2041
2042 ! mask for measured runoff
2043 logical, dimension(:), allocatable :: runoff_obs_mask
2044
2045 ! run mHM
2046 allocate(smoptisim(domainmeta%nDomains))
2047 call eval(parameterset, runoff = runoff, smoptisim = smoptisim)
2048
2049 ! -----------------------------
2050 ! SOIL MOISTURE
2051 ! -----------------------------
2052
2053 ! initialize some variables
2054 objective_sm = 0.0_dp
2055
2056 ! loop over domain - for applying power law later on
2057 do idomain = 1, domainmeta%nDomains
2058
2059 ! init
2060 objective_sm_domain = 0.0_dp
2061 ! get domain information
2062 ncells1 = level1(idomain)%nCells
2063
2064 ! correlation signal is calculated on individual grid cells
2065 invalid_cells = 0.0_dp
2066 do icell = 1, size(l1_smobs(idomain)%maskObs(:, :), dim = 1)
2067
2068 ! check for enough data points in time for statistical calculations (e.g. mean, stddev)
2069 if (count(l1_smobs(idomain)%maskObs(icell, :)) .LE. (0.10_dp * real(size(l1_smobs(idomain)%dataObs, dim = 2), dp))) then
2070 invalid_cells = invalid_cells + 1.0_dp
2071 cycle
2072 end if
2073
2074 ! calculate ojective function
2075 objective_sm_domain = objective_sm_domain + &
2076 correlation(l1_smobs(idomain)%dataObs(icell, :), smoptisim(idomain)%dataSim(icell, :), &
2077 mask = l1_smobs(idomain)%maskObs(icell, :))
2078 end do
2079
2080 ! user information about invalid cells
2081 if (invalid_cells .GT. 0.5_dp) then
2082 call message(.LT.' WARNING: objective_sm: Detected invalid cells in study area ( 10 valid data points).')
2083 call message(' Fraction of invalid cells: ', &
2084 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
2085 end if
2086
2087 ! calculate average soil moisture objective over all domains with power law
2088 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
2089 objective_sm = objective_sm + &
2090 ((1.0_dp - objective_sm_domain / real(ncells1, dp)) / real(domainmeta%overallNumberOfDomains, dp))**6
2091 call smoptisim(idomain)%destroy()
2092 end do
2093 deallocate(smoptisim)
2094
2095 ! compromise solution - sixth root
2096 objective_sm = objective_sm**onesixth
2097
2098 ! -----------------------------
2099 ! RUNOFF
2100 ! -----------------------------
2101 objective_kge = 0.0_dp
2102 ngaugestotal = size(runoff, dim = 2)
2103
2104 do gg = 1, ngaugestotal
2105
2106 ! extract runoff
2107 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2108
2109 ! KGE
2110 objective_kge = objective_kge + &
2111 ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(ngaugestotal, dp))**6
2112
2113 end do
2114
2115 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2116
2117 ! compromise solution - sixth root
2118 objective_kge = objective_kge**onesixth
2119
2120 ! equal weighted compromise objective functions for discharge and soilmoisture
2121 ! ToDo: why do we use the sixth root of of objective_sm and objective_kge
2122 ! only to take the power to 6 here again, when we never need the
2123 ! intermediate results?
2124#ifdef MPI
2125 objective_kge_q_sm_corr = (objective_sm**6 + objective_kge**6)
2126#else
2127 objective_kge_q_sm_corr = (objective_sm**6 + objective_kge**6)**onesixth
2128
2129 call message(' objective_kge_q_sm_corr = ', num2str(objective_kge_q_sm_corr, '(F9.5)'))
2130#endif
2131 ! print*, "1-SM 2-Q : ", 1.0_dp-objective_sm, 1.0_dp-objective_kge ! MZMZMZMZ
2132
2133 END FUNCTION objective_kge_q_sm_corr
2134
2135
2136 ! -----------------------------------------------------------------
2137
2138 ! NAME
2139 ! objective_kge_q_et
2140
2141 ! PURPOSE
2142 !> \brief Objective function of KGE for runoff and KGE for ET
2143
2144 !> \details Objective function of KGE for runoff and KGE for ET.
2145 !> Further details can be found in the documentation of objective functions
2146 !> '14 - objective_multiple_gauges_kge_power6'.
2147
2148 ! INTENT(IN)
2149 !> \param[in] "real(dp), dimension(:) :: parameterset"
2150 !> \param[in] "procedure(eval_interface) :: eval"
2151
2152 ! RETURN
2153 !> \return real(dp) :: objective_kge_q_et — objective function value
2154 !> (which will be e.g. minimized by an optimization routine like DDS)
2155
2156 ! HISTORY
2157 !> \authors Johannes Brenner
2158
2159 !> \date July 2017
2160
2161 ! Modifications:
2162 ! Robert Schweppe Jun 2018 - refactoring and reformatting
2163
2164 FUNCTION objective_kge_q_et(parameterset, eval)
2165
2166 use mo_optimization_types, only : optidata_sim
2167 use mo_common_variables, only : level1, domainmeta
2168 use mo_errormeasures, only : kge
2169 use mo_global_variables, only : l1_etobs
2170 use mo_string_utils, only : num2str
2171 use mo_mrm_objective_function_runoff, only : extract_runoff
2172
2173 implicit none
2174
2175 real(dp), dimension(:), intent(in) :: parameterset
2176
2177 procedure(eval_interface), INTENT(IN), POINTER :: eval
2178
2179 real(dp) :: objective_kge_q_et
2180
2181 real(dp) :: objective_et
2182
2183 real(dp) :: objective_q
2184
2185 ! number of invalid cells in catchment
2186 real(dp) :: invalid_cells
2187
2188 ! modelled runoff for a given parameter set
2189 ! dim1=nTimeSteps, dim2=nGauges
2190 real(dp), allocatable, dimension(:, :) :: runoff
2191
2192 ! domain loop counter
2193 integer(i4) :: idomain
2194
2195 ! cell loop counter
2196 integer(i4) :: icell
2197
2198 ! ncells1 of level 1
2199 integer(i4) :: ncells1
2200
2201 ! domains wise objectives
2202 real(dp) :: objective_et_domain
2203
2204 !> simulated et
2205 type(optidata_sim), dimension(:), allocatable :: etoptisim
2206
2207 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2208
2209 ! gauges counter
2210 integer(i4) :: gg
2211
2212 integer(i4) :: ngaugestotal
2213
2214 ! aggregated simulated runoff
2215 real(dp), dimension(:), allocatable :: runoff_agg
2216
2217 ! measured runoff
2218 real(dp), dimension(:), allocatable :: runoff_obs
2219
2220 ! mask for measured runoff
2221 logical, dimension(:), allocatable :: runoff_obs_mask
2222
2223 ! run mHM
2224 allocate(etoptisim(domainmeta%nDomains))
2225 call eval(parameterset, runoff = runoff, etoptisim = etoptisim)
2226
2227 ! -----------------------------
2228 ! EVAPOTRANSPIRATION
2229 ! -----------------------------
2230
2231 ! initialize some variables
2232 objective_et = 0.0_dp
2233
2234 ! loop over domain - for applying power law later on
2235 do idomain = 1, domainmeta%nDomains
2236
2237 ! init
2238 objective_et_domain = 0.0_dp
2239 ! get domain information
2240 ncells1 = level1(idomain)%nCells
2241
2242 ! correlation signal is calculated on individual grid cells
2243 invalid_cells = 0.0_dp
2244 do icell = 1, size(l1_etobs(idomain)%maskObs, dim = 1)
2245
2246 ! check for enough data points in time for statistical calculations (e.g. mean, stddev)
2247 if (count(l1_etobs(idomain)%maskObs(icell, :)) .LE. &
2248 (0.10_dp * real(size(l1_etobs(idomain)%dataObs(:, :), dim = 2), dp))) then
2249 invalid_cells = invalid_cells + 1.0_dp
2250 cycle
2251 end if
2252
2253 ! calculate ojective function
2254 objective_et_domain = objective_et_domain + &
2255 kge(l1_etobs(idomain)%dataObs(icell, :), etoptisim(idomain)%dataSim(icell, :), &
2256 mask = l1_etobs(idomain)%maskObs(icell, :))
2257 end do
2258
2259 ! user information about invalid cells
2260 if (invalid_cells .GT. 0.5_dp) then
2261 call message(.LT.' WARNING: objective_et: Detected invalid cells in study area ( 10 valid data points).')
2262 call message(' Fraction of invalid cells: ', &
2263 num2str(invalid_cells / real(ncells1, dp), '(F4.2)'))
2264 end if
2265
2266 ! calculate average soil moisture objective over all domains with power law
2267 ! domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
2268 objective_et = objective_et + &
2269 ((1.0_dp - objective_et_domain / real(ncells1, dp)) / real(domainmeta%overallNumberOfDomains, dp))**6
2270 call etoptisim(idomain)%destroy()
2271 end do
2272 deallocate(etoptisim)
2273
2274 ! compromise solution - sixth root
2275 objective_et = objective_et**onesixth
2276
2277 ! -----------------------------
2278 ! RUNOFF
2279 ! -----------------------------
2280 objective_q = 0.0_dp
2281 ngaugestotal = size(runoff, dim = 2)
2282
2283 do gg = 1, ngaugestotal
2284
2285 ! extract runoff
2286 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2287
2288 ! KGE
2289 objective_q = objective_q + &
2290 ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(ngaugestotal, dp))**6
2291
2292 end do
2293
2294 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2295
2296 ! compromise solution - sixth root
2297 objective_q = objective_q**onesixth
2298
2299 ! equal weighted compromise objective functions for discharge and soilmoisture
2300 ! ToDo: why do we use the sixth root of of objective_sm and objective_kge
2301 ! only to take the power to 6 here again, when we never need the
2302 ! intermediate results?
2303#ifdef MPI
2304 objective_kge_q_et = (objective_et**6 + objective_q**6)
2305#else
2306 objective_kge_q_et = (objective_et**6 + objective_q**6)**onesixth
2307
2308 call message(' objective_kge_q_et = ', num2str(objective_kge_q_et, '(F9.5)'))
2309#endif
2310 ! print*, "1-SM 2-Q : ", 1.0_dp-objective_sm, 1.0_dp-objective_kge ! MZMZMZMZ
2311
2312 END FUNCTION objective_kge_q_et
2313
2314
2315 ! -----------------------------------------------------------------
2316
2317 ! NAME
2318 ! objective_kge_q_BFI
2319
2320 ! PURPOSE
2321 !> \brief Objective function of KGE for runoff and BFI absulute difference
2322
2323 !> \details Objective function of KGE for runoff and KGE for ET.
2324 !> Further details can be found in the documentation of objective functions
2325 !> '14 - objective_multiple_gauges_kge_power6'.
2326
2327 ! INTENT(IN)
2328 !> \param[in] "real(dp), dimension(:) :: parameterset"
2329 !> \param[in] "procedure(eval_interface) :: eval"
2330
2331 ! RETURN
2332 !> \return real(dp) :: objective_kge_q_BFI — objective function value
2333 !> (which will be e.g. minimized by an optimization routine like DDS)
2334
2335 ! HISTORY
2336 !> \authors Sebastian Müller
2337
2338 !> \date Apr 2022
2339
2340 FUNCTION objective_kge_q_bfi(parameterset, eval)
2341
2342 use mo_optimization_types, only : optidata_sim
2343 use mo_common_variables, only : level1, domainmeta
2344 use mo_errormeasures, only : kge
2345 use mo_global_variables, only : bfi_obs
2346 use mo_string_utils, only : num2str
2347 use mo_mrm_objective_function_runoff, only : extract_runoff
2348
2349 implicit none
2350
2351 real(dp), dimension(:), intent(in) :: parameterset
2352
2353 procedure(eval_interface), INTENT(IN), POINTER :: eval
2354
2355 real(dp) :: objective_kge_q_bfi
2356 real(dp) :: objective_bfi
2357 real(dp) :: objective_q
2358
2359 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2360
2361 ! number of invalid cells in catchment
2362 real(dp) :: invalid_cells
2363
2364 ! modelled runoff for a given parameter set
2365 ! dim1=nTimeSteps, dim2=nGauges
2366 real(dp), allocatable, dimension(:, :) :: runoff
2367
2368 ! domain loop counter
2369 integer(i4) :: idomain
2370
2371 !> baseflow index for each domain
2372 real(dp), dimension(:), allocatable :: bfi
2373
2374 ! counter
2375 integer(i4) :: gg, i
2376
2377 integer(i4) :: ngaugestotal
2378
2379 ! aggregated simulated runoff
2380 integer(i4), dimension(:), allocatable :: domain_ids, domain_ids_pack
2381
2382 ! aggregated simulated runoff
2383 real(dp), dimension(:), allocatable :: runoff_agg
2384
2385 ! measured runoff
2386 real(dp), dimension(:), allocatable :: runoff_obs
2387
2388 ! mask for measured runoff
2389 logical, dimension(:), allocatable :: runoff_obs_mask
2390
2391 ! run mHM
2392 allocate(bfi(domainmeta%nDomains))
2393 call eval(parameterset, runoff = runoff, bfi = bfi)
2394
2395 ! -----------------------------
2396 ! BFI
2397 ! -----------------------------
2398
2399 ! initialize some variables
2400 objective_bfi = 0.0_dp
2401
2402 if ( any(bfi_obs < 0.0_dp) ) then
2403 allocate(domain_ids(domainmeta%nDomains))
2404 allocate(domain_ids_pack(count(bfi_obs < 0.0_dp)))
2405 domain_ids = [(i, i=1,size(domain_ids))]
2406 domain_ids_pack = pack(domain_ids, mask=(bfi_obs < 0.0_dp))
2407 call error_message( &
2408 "objective_kge_q_BFI: missing BFI values for domain ", &
2409 trim(adjustl(num2str(domain_ids_pack(1)))) &
2410 )
2411 end if
2412
2413 ! loop over domain - for applying power law later on
2414 do idomain = 1, domainmeta%nDomains
2415 objective_bfi = objective_bfi + abs(bfi(idomain) - bfi_obs(idomain)) / domainmeta%nDomains
2416 end do
2417
2418 ! -----------------------------
2419 ! RUNOFF
2420 ! -----------------------------
2421 objective_q = 0.0_dp
2422 ngaugestotal = size(runoff, dim = 2)
2423
2424 do gg = 1, ngaugestotal
2425
2426 ! extract runoff
2427 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2428
2429 ! KGE
2430 objective_q = objective_q + &
2431 ((1.0_dp - kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)) / real(ngaugestotal, dp))**6
2432
2433 end do
2434
2435 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
2436
2437 ! compromise solution - sixth root
2438 objective_q = objective_q**onesixth
2439
2440 objective_kge_q_bfi = (objective_bfi + 1._dp)*objective_q
2441 call message(' objective_kge_q_BFI = ', num2str(objective_kge_q_bfi, '(F9.5)'))
2442
2443 END FUNCTION objective_kge_q_bfi
2444
2445
2446 ! -----------------------------------------------------------------
2447
2448 ! NAME
2449 ! objective_kge_q_rmse_et
2450
2451 ! PURPOSE
2452 !> \brief Objective function of KGE for runoff and RMSE for domain_avg ET (standarized scores)
2453
2454 !> \details Objective function of KGE for runoff and RMSE for domain_avg ET (standarized scores)
2455
2456 ! INTENT(IN)
2457 !> \param[in] "real(dp), dimension(:) :: parameterset"
2458 !> \param[in] "procedure(eval_interface) :: eval"
2459
2460 ! RETURN
2461 !> \return real(dp) :: objective_kge_q_rmse_et &mdash; objective function value
2462 !> (which will be e.g. minimized by an optimization routine like DDS)
2463
2464 ! HISTORY
2465 !> \authors Johannes Brenner
2466
2467 !> \date July 2017
2468
2469 ! Modifications:
2470 ! Robert Schweppe Jun 2018 - refactoring and reformatting
2471
2472 FUNCTION objective_kge_q_rmse_et(parameterset, eval)
2473
2474 use mo_optimization_types, only : optidata_sim
2475 use mo_common_constants, only : eps_dp, nodata_dp
2476 use mo_common_mhm_mrm_variables, only : evalper
2477 use mo_common_variables, only : domainmeta
2478 use mo_errormeasures, only : rmse
2479 use mo_global_variables, only : l1_etobs
2480 use mo_julian, only : caldat
2481 use mo_moment, only : average, mean
2482 use mo_standard_score, only : classified_standard_score
2483 use mo_string_utils, only : num2str
2484 use mo_temporal_aggregation, only : day2mon_average
2485 use mo_errormeasures, only : kge
2486 use mo_mrm_objective_function_runoff, only : extract_runoff
2487
2488 implicit none
2489
2490 real(dp), dimension(:), intent(in) :: parameterset
2491
2492 procedure(eval_interface), INTENT(IN), POINTER :: eval
2493
2494 real(dp) :: objective_kge_q_rmse_et
2495
2496 ! modelled runoff for a given parameter set
2497 ! dim1=nTimeSteps, dim2=nGauges
2498 real(dp), allocatable, dimension(:, :) :: runoff
2499
2500 !> simulated et
2501 type(optidata_sim), dimension(:), allocatable :: etoptisim
2502
2503 ! time loop counter
2504 integer(i4) :: itime
2505
2506 ! domain counter, month counters
2507 integer(i4) :: idomain, pp, mmm
2508
2509 integer(i4) :: year, month, day
2510
2511 real(dp), dimension(domainMeta%nDomains) :: inittime
2512
2513 ! total number of months
2514 integer(i4) :: nmonths
2515
2516 ! vector with months' classes
2517 integer(i4), dimension(:), allocatable :: month_classes
2518
2519 ! monthly values original time series
2520 real(dp), dimension(:), allocatable :: et_sim_m, et_obs_m
2521
2522 ! monthly values anomaly time series
2523 real(dp), dimension(:), allocatable :: et_sim_m_anom, et_obs_m_anom
2524
2525 logical, dimension(:), allocatable :: et_obs_m_mask
2526
2527 ! kge_q(nGaugesTotal)
2528 real(dp), dimension(:), allocatable :: rmse_et
2529
2530 ! obj. functions
2531 real(dp) :: rmse_et_avg, kge_q_avg
2532
2533 ! spatial average of observed et
2534 real(dp), dimension(:), allocatable :: et_catch_avg_domain
2535
2536 ! spatial avergae of modeled et
2537 real(dp), dimension(:), allocatable :: et_opti_catch_avg_domain
2538
2539 ! mask for valid et catchment avg time steps
2540 logical, dimension(:), allocatable :: mask_times
2541
2542 ! rmse_et(domainMeta%nDomains)
2543 real(dp), dimension(:), allocatable :: kge_q
2544
2545 ! gauges counter
2546 integer(i4) :: gg
2547
2548 integer(i4) :: ngaugestotal
2549
2550 ! aggregated simulated runoff
2551 real(dp), dimension(:), allocatable :: runoff_agg
2552
2553 ! measured runoff
2554 real(dp), dimension(:), allocatable :: runoff_obs
2555
2556 ! mask for measured runoff
2557 logical, dimension(:), allocatable :: runoff_obs_mask
2558
2559 ! obtain simulation values of runoff (hourly) and ET
2560 ! for ET only valid cells (domains concatenated)
2561 ! et_opti: aggregate ET to needed time step for optimization (see timeStep_et_input)
2562 allocate(etoptisim(domainmeta%nDomains))
2563 call eval(parameterset, runoff = runoff, etoptisim = etoptisim)
2564
2565 !--------------------------------------------
2566 !! EVAPOTRANSPIRATION
2567 !--------------------------------------------
2568
2569 ! allocate
2570 allocate(rmse_et(domainmeta%nDomains))
2571 rmse_et(:) = nodata_dp
2572
2573 do idomain = 1, domainmeta%nDomains
2574
2575 ! allocate
2576 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
2577 allocate(et_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
2578 allocate(et_opti_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
2579
2580 ! initalize
2581 mask_times = .true.
2582 et_catch_avg_domain = nodata_dp
2583 et_opti_catch_avg_domain = nodata_dp
2584
2585 ! calculate catchment average evapotranspiration
2586 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
2587 ! check for enough data points in time for correlation
2588 if (all(.NOT. l1_etobs(idomain)%maskObs(:, itime))) then
2589 !write (*,*) 'WARNING: et data at time ', iTime, ' is empty.'
2590 !call message('WARNING: objective_et_kge_catchment_avg: ignored current time step since less than')
2591 !call message(' 10 valid cells available in evapotranspiration observation')
2592 mask_times(itime) = .false.
2593 cycle
2594 end if
2595 ! spatial average of observed ET
2596 et_catch_avg_domain(itime) = average(l1_etobs(idomain)%dataObs(:, itime), &
2597 mask = l1_etobs(idomain)%maskObs(:, itime))
2598 ! spatial avergae of modeled ET
2599 et_opti_catch_avg_domain(itime) = average(etoptisim(idomain)%dataSim(:, itime), &
2600 mask = l1_etobs(idomain)%maskObs(:, itime))
2601 end do
2602
2603 ! get initial time of the evaluation period
2604 inittime(idomain) = real(evalper(idomain)%julStart, dp)
2605
2606 ! get calendar days, months, year
2607 call caldat(int(inittime(idomain)), yy = year, mm = month, dd = day)
2608
2609 ! if evapotranspiration input daily
2610 select case(l1_etobs(idomain)%timeStepInput)
2611 ! JBJBJB
2612 ! daily: aggregate to monthly mean
2613 case(-1)
2614 ! calculate monthly averages from daily values of the model
2615 call day2mon_average(et_opti_catch_avg_domain, year, month, day, et_sim_m, misval = nodata_dp)
2616 ! calculate monthly averages from daily values of the observations
2617 call day2mon_average(et_catch_avg_domain, year, month, day, et_obs_m, misval = nodata_dp)
2618 !
2619 ! monthly: proceed without action
2620 case(-2)
2621 ! simulation
2622 allocate(et_sim_m(size(etoptisim(idomain)%dataSim, dim = 2)))
2623 et_sim_m = et_opti_catch_avg_domain
2624 ! observation
2625 allocate(et_obs_m(size(etoptisim(idomain)%dataSim, dim = 2)))
2626 et_obs_m = et_catch_avg_domain
2627
2628 ! yearly: ERROR stop program
2629 case(-3)
2630 call error_message('***ERROR: objective_kge_q_rmse_et: time step of evapotranspiration yearly.')
2631 end select
2632 ! remove mean from modelled time series
2633 et_sim_m(:) = et_sim_m(:) - mean(et_sim_m(:))
2634 ! remove mean from observed time series
2635 et_obs_m(:) = et_obs_m(:) - mean(et_obs_m(:))
2636 ! get number of months for given domain
2637 nmonths = size(et_obs_m)
2638 allocate (month_classes(nmonths))
2639 allocate (et_obs_m_mask(nmonths))
2640 allocate (et_obs_m_anom(nmonths))
2641 allocate (et_sim_m_anom(nmonths))
2642
2643 month_classes(:) = 0
2644 et_obs_m_mask(:) = .true.
2645 et_obs_m_anom(:) = nodata_dp
2646 et_sim_m_anom(:) = nodata_dp
2647
2648 ! define months' classes
2649 mmm = month
2650 do pp = 1, nmonths
2651 month_classes(pp) = mmm
2652 if (mmm .LT. 12) then
2653 mmm = mmm + 1
2654 else
2655 mmm = 1
2656 end if
2657 end do
2658 ! double check missing data
2659 ! define mask for missing data in observations (there are always data for simulations)
2660 where(abs(et_obs_m - nodata_dp) .lt. eps_dp) et_obs_m_mask = .false.
2661
2662 ! calculate standard score
2663 et_obs_m_anom = classified_standard_score(et_obs_m, month_classes, mask = et_obs_m_mask)
2664 et_sim_m_anom = classified_standard_score(et_sim_m, month_classes, mask = et_obs_m_mask)
2665 rmse_et(idomain) = rmse(et_sim_m_anom, et_obs_m_anom, mask = et_obs_m_mask)
2666
2667 deallocate (month_classes)
2668 deallocate (et_obs_m)
2669 deallocate (et_sim_m)
2670 deallocate (et_obs_m_mask)
2671 deallocate (et_sim_m_anom)
2672 deallocate (et_obs_m_anom)
2673 !end if
2674
2675 call etoptisim(idomain)%destroy()
2676 end do
2677
2678 rmse_et_avg = sum(rmse_et(:), abs(rmse_et - nodata_dp) .gt. eps_dp) / &
2679 real(count(abs(rmse_et - nodata_dp) .gt. eps_dp), dp)
2680 deallocate(rmse_et)
2681 deallocate(etoptisim)
2682
2683 !--------------------------------------------
2684 !! RUNOFF
2685 !--------------------------------------------
2686 kge_q_avg = 0_dp
2687 ngaugestotal = size(runoff, dim = 2)
2688 allocate(kge_q(ngaugestotal))
2689 kge_q(:) = nodata_dp
2690
2691 do gg = 1, ngaugestotal
2692
2693 ! extract runoff
2694 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
2695
2696 ! check for whether to proceed with this domain or not
2697 ! potentially 3 years of data
2698 !pp = count(runoff_agg .ge. 0.0_dp )
2699 !if (pp .lt. 365*3 ) then
2700 ! deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
2701 ! cycle
2702 ! else
2703 ! calculate KGE for each domain:
2704 kge_q(gg) = kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
2705 deallocate (runoff_agg, runoff_obs, runoff_obs_mask)
2706 ! end if
2707
2708 end do
2709
2710 ! calculate average KGE value for runoff
2711 kge_q_avg = sum(kge_q(:), abs(kge_q - nodata_dp) .gt. eps_dp) / &
2712 real(count(abs(kge_q - nodata_dp) .gt. eps_dp), dp)
2713 deallocate(kge_q)
2714
2715 !
2716 objective_kge_q_rmse_et = rmse_et_avg * (1._dp - kge_q_avg)
2717
2718 call message(' objective_kge_q_rmse_et = ', num2str(objective_kge_q_rmse_et, '(F9.5)'))
2719
2720 END FUNCTION objective_kge_q_rmse_et
2721
2722 ! ------------------------------------------------------------------
2723
2724 ! NAME
2725 ! objective_spaef_sm
2726
2727 ! PURPOSE
2728 !> \brief Objective function for soil moisture.
2729
2730 !> \details The Spatial Efficiency (SPAEF) of observed and modeled soil
2731 !> moisture fields is calculated - aim: matching spatial patterns
2732
2733 ! INTENT(IN)
2734 !> \param[in] "real(dp), dimension(:) :: parameterset"
2735 !> \param[in] "procedure(eval_interface) :: eval"
2736
2737 ! RETURN
2738 !> \return real(dp) :: objective_spaef_sm &mdash; objective function value
2739 !> (which will be e.g. minimized by an optimization routine like DDS)
2740
2741 ! LITERATURE
2742 ! Koch, J., Demirel, M. C., & Stisen, S. (2018). The SPAtial EFficiency metric (SPAEF):
2743 ! Multiple-component evaluation of spatial patterns for optimization of hydrological models.
2744 ! Geoscientific Model Development, 11(5), 1873–1886. https://doi.org/10.5194/gmd-11-1873-2018
2745
2746
2747 ! HISTORY
2748 !> \authors Pallav Shrestha
2749
2750 !> \date Dec 2024
2751
2752 ! Modifications:
2753
2754
2755 FUNCTION objective_spaef_sm(parameterset, eval)
2756
2757 use mo_optimization_types, only : optidata_sim
2758 use mo_common_constants, only : nodata_dp
2759 use mo_common_variables, only : level1, domainmeta
2760 use mo_global_variables, only : l1_smobs
2761 use mo_message, only : message
2762 use mo_errormeasures, only : spaef
2763 use mo_string_utils, only : num2str
2764
2765 implicit none
2766
2767 real(dp), dimension(:), intent(in) :: parameterset
2768
2769 procedure(eval_interface), INTENT(IN), POINTER :: eval
2770
2771 ! objective function value
2772 real(dp) :: objective_spaef_sm
2773
2774 ! domain loop counter
2775 integer(i4) :: idomain
2776
2777 ! time loop counter
2778 integer(i4) :: itime
2779
2780 ! for sixth root
2781#ifndef MPI
2782 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2783#endif
2784
2785 ! vectorized arrays of SM
2786 real(dp), dimension(:), allocatable :: vec1, vec2
2787
2788 ! spaef at every time step
2789 real(dp), dimension(:), allocatable :: spaef_time_series
2790
2791 ! simulated soil moisture
2792 type(optidata_sim), dimension(:), allocatable :: smoptisim
2793
2794 ! mask of valid sm cells
2795 logical, dimension(:), allocatable :: mask_sm
2796
2797 ! mask for valid sm time steps
2798 logical, dimension(:), allocatable :: mask_times
2799
2800
2801 allocate(smoptisim(domainmeta%nDomains))
2802 call eval(parameterset, smoptisim = smoptisim)
2803
2804 ! initialize some variables
2805 objective_spaef_sm = 0.0_dp
2806
2807 ! loop over domain - for applying power law later on
2808 do idomain = 1, domainmeta%nDomains
2809
2810 ! allocate
2811 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
2812 allocate(spaef_time_series(size(smoptisim(idomain)%dataSim, dim = 2)))
2813 allocate(vec1(size(smoptisim(idomain)%dataSim, dim = 1)))
2814 allocate(vec2(size(smoptisim(idomain)%dataSim, dim = 1)))
2815 allocate(mask_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
2816
2817 ! initalize
2818 mask_times = .false.
2819 spaef_time_series = 0.0_dp
2820
2821 ! calculate the SPAEF criterion across each time step
2822 do itime = 1, size(smoptisim(idomain)%dataSim, dim = 2)
2823 vec1 = l1_smobs(idomain)%dataObs(:, itime)
2824 vec2 = smoptisim(idomain)%dataSim(:, itime)
2825 mask_sm = l1_smobs(idomain)%maskObs(:, itime)
2826 if (count(mask_sm) > 1_i4) then ! sample size must be at least 2
2827 mask_times(itime) = .true.
2828 spaef_time_series(itime) = spaef(vec1, vec2, mask = mask_sm)
2829 end if
2830 end do
2831
2832 if (count(mask_times) > 0_i4) then
2833 ! calculate the temporal average of SPAEF, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
2834 objective_spaef_sm = objective_spaef_sm + &
2835 ((1.0_dp - sum(spaef_time_series, mask = mask_times) / real(count(mask_times), dp)) / &
2836 real(domainmeta%overallnumberofdomains, dp))**6
2837 else
2838 call message('***ERROR: mo_objective_funtion: objective_spaef_sm: No soil moisture observations available!')
2839 stop
2840 end if
2841
2842 ! deallocate
2843 deallocate(mask_times)
2844 deallocate(spaef_time_series)
2845 deallocate(vec1)
2846 deallocate(vec2)
2847 deallocate(mask_sm)
2848 call smoptisim(idomain)%destroy()
2849 end do
2850 deallocate(smoptisim)
2851
2852#ifndef MPI
2853 objective_spaef_sm = objective_spaef_sm**onesixth
2854
2855 call message(' objective_spaef_sm = ', num2str(objective_spaef_sm, '(F9.5)'))
2856 call message(' spaef_sm = ', num2str(1 - objective_spaef_sm, '(F9.5)'))
2857#endif
2858
2859 END FUNCTION objective_spaef_sm
2860
2861 ! ------------------------------------------------------------------
2862
2863 ! NAME
2864 ! objective_spaef_et
2865
2866 ! PURPOSE
2867 !> \brief Objective function for ET.
2868
2869 !> \details The Spatial Efficiency (SPAEF) of observed and modeled ET
2870 !> fields is calculated - aim: matching spatial patterns
2871
2872 ! INTENT(IN)
2873 !> \param[in] "real(dp), dimension(:) :: parameterset"
2874 !> \param[in] "procedure(eval_interface) :: eval"
2875
2876 ! RETURN
2877 !> \return real(dp) :: objective_spaef_et &mdash; objective function value
2878 !> (which will be e.g. minimized by an optimization routine like DDS)
2879
2880 ! LITERATURE
2881 ! Koch, J., Demirel, M. C., & Stisen, S. (2018). The SPAtial EFficiency metric (SPAEF):
2882 ! Multiple-component evaluation of spatial patterns for optimization of hydrological models.
2883 ! Geoscientific Model Development, 11(5), 1873–1886. https://doi.org/10.5194/gmd-11-1873-2018
2884
2885
2886 ! HISTORY
2887 !> \authors Pallav Shrestha
2888
2889 !> \date Dec 2024
2890
2891 ! Modifications:
2892
2893
2894 FUNCTION objective_spaef_et(parameterset, eval)
2895
2896 use mo_optimization_types, only : optidata_sim
2897 use mo_common_constants, only : nodata_dp
2898 use mo_common_variables, only : level1, domainmeta
2899 use mo_global_variables, only : l1_etobs
2900 use mo_message, only : message
2901 use mo_errormeasures, only : spaef
2902 use mo_string_utils, only : num2str
2903
2904 implicit none
2905
2906 real(dp), dimension(:), intent(in) :: parameterset
2907
2908 procedure(eval_interface), INTENT(IN), POINTER :: eval
2909
2910 ! objective function value
2911 real(dp) :: objective_spaef_et
2912
2913 ! domain loop counter
2914 integer(i4) :: idomain
2915
2916 ! time loop counter
2917 integer(i4) :: itime
2918
2919 ! for sixth root
2920#ifndef MPI
2921 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
2922#endif
2923
2924 ! vectorized arrays of ET
2925 real(dp), dimension(:), allocatable :: vec1, vec2
2926
2927 ! spaef at every time step
2928 real(dp), dimension(:), allocatable :: spaef_time_series
2929
2930 ! simulated ET
2931 type(optidata_sim), dimension(:), allocatable :: etoptisim
2932
2933 ! mask of valid ET cells
2934 logical, dimension(:), allocatable :: mask_et
2935
2936 ! mask for valid ET time steps
2937 logical, dimension(:), allocatable :: mask_times
2938
2939
2940 allocate(etoptisim(domainmeta%nDomains))
2941 call eval(parameterset, etoptisim = etoptisim)
2942
2943 ! initialize some variables
2944 objective_spaef_et = 0.0_dp
2945
2946 ! loop over domain - for applying power law later on
2947 do idomain = 1, domainmeta%nDomains
2948
2949 ! allocate
2950 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
2951 allocate(spaef_time_series(size(etoptisim(idomain)%dataSim, dim = 2)))
2952 allocate(vec1(size(etoptisim(idomain)%dataSim, dim = 1)))
2953 allocate(vec2(size(etoptisim(idomain)%dataSim, dim = 1)))
2954 allocate(mask_et(size(etoptisim(idomain)%dataSim, dim = 1)))
2955
2956 ! initalize
2957 mask_times = .false.
2958 spaef_time_series = 0.0_dp
2959
2960 ! calculate the SPAEF criterion across each time step
2961 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
2962 vec1 = l1_etobs(idomain)%dataObs(:, itime)
2963 vec2 = etoptisim(idomain)%dataSim(:, itime)
2964 mask_et = l1_etobs(idomain)%maskObs(:, itime)
2965 if (count(mask_et) > 1_i4) then ! sample size must be at least 2
2966 mask_times(itime) = .true.
2967 spaef_time_series(itime) = spaef(vec1, vec2, mask = mask_et)
2968 end if
2969 end do
2970
2971 if (count(mask_times) > 0_i4) then
2972 ! calculate the temporal average of SPAEF, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
2973 objective_spaef_et = objective_spaef_et + &
2974 ((1.0_dp - sum(spaef_time_series, mask = mask_times) / real(count(mask_times), dp)) / &
2975 real(domainmeta%overallnumberofdomains, dp))**6
2976 else
2977 call message('***ERROR: mo_objective_funtion: objective_spaef_et: No ET observations available!')
2978 stop
2979 end if
2980
2981 ! deallocate
2982 deallocate(mask_times)
2983 deallocate(spaef_time_series)
2984 deallocate(vec1)
2985 deallocate(vec2)
2986 deallocate(mask_et)
2987 call etoptisim(idomain)%destroy()
2988 end do
2989 deallocate(etoptisim)
2990
2991#ifndef MPI
2992 objective_spaef_et = objective_spaef_et**onesixth
2993
2994 call message(' objective_spaef_et = ', num2str(objective_spaef_et, '(F9.5)'))
2995 call message(' spaef_et = ', num2str(1 - objective_spaef_et, '(F9.5)'))
2996#endif
2997
2998 END FUNCTION objective_spaef_et
2999
3000 ! ------------------------------------------------------------------
3001
3002 ! NAME
3003 ! objective_spaef_sm_spaef_et
3004
3005 ! PURPOSE
3006 !> \brief Bivariate Objective function for SM and ET.
3007
3008 !> \details The Spatial Efficiency (SPAEF) of observed and modeled SM and ET
3009 !> fields (equal weights) is calculated - aim: matching spatial patterns
3010
3011 ! INTENT(IN)
3012 !> \param[in] "real(dp), dimension(:) :: parameterset"
3013 !> \param[in] "procedure(eval_interface) :: eval"
3014
3015 ! RETURN
3016 !> \return real(dp) :: objective_spaef_sm_spaef_et &mdash; objective function value
3017 !> (which will be e.g. minimized by an optimization routine like DDS)
3018
3019 ! LITERATURE
3020 ! Koch, J., Demirel, M. C., & Stisen, S. (2018). The SPAtial EFficiency metric (SPAEF):
3021 ! Multiple-component evaluation of spatial patterns for optimization of hydrological models.
3022 ! Geoscientific Model Development, 11(5), 1873–1886. https://doi.org/10.5194/gmd-11-1873-2018
3023
3024
3025 ! HISTORY
3026 !> \authors Pallav Shrestha
3027
3028 !> \date Dec 2024
3029
3030 ! Modifications:
3031
3032
3033 FUNCTION objective_spaef_sm_spaef_et(parameterset, eval)
3034
3035 use mo_optimization_types, only : optidata_sim
3036 use mo_common_constants, only : nodata_dp
3037 use mo_common_variables, only : level1, domainmeta
3038 use mo_global_variables, only : l1_smobs, l1_etobs
3039 use mo_message, only : message
3040 use mo_errormeasures, only : spaef
3041 use mo_string_utils, only : num2str
3042
3043 implicit none
3044
3045 real(dp), dimension(:), intent(in) :: parameterset
3046
3047 procedure(eval_interface), INTENT(IN), POINTER :: eval
3048
3049 ! objective function value
3050 real(dp) :: objective_spaef_sm_spaef_et
3051
3052 ! domain loop counter
3053 integer(i4) :: idomain
3054
3055 ! time loop counter
3056 integer(i4) :: itime
3057
3058 ! for sixth root
3059#ifndef MPI
3060 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
3061#endif
3062
3063 ! vectorized arrays of SM and ET
3064 real(dp), dimension(:), allocatable :: vec1_sm, vec2_sm, vec1_et, vec2_et
3065
3066 ! spaef at every time step
3067 real(dp), dimension(:), allocatable :: spaef_time_series_sm, spaef_time_series_et
3068
3069 ! simulated SM
3070 type(optidata_sim), dimension(:), allocatable :: smoptisim
3071
3072 ! simulated ET
3073 type(optidata_sim), dimension(:), allocatable :: etoptisim
3074
3075 ! mask of valid cells
3076 logical, dimension(:), allocatable :: mask_et, mask_sm
3077
3078 ! mask for valid time steps
3079 logical, dimension(:), allocatable :: mask_times
3080
3081
3082 allocate(smoptisim(domainmeta%nDomains))
3083 allocate(etoptisim(domainmeta%nDomains))
3084 call eval(parameterset, smoptisim = smoptisim, etoptisim = etoptisim)
3085
3086 ! initialize some variables
3087 objective_spaef_sm_spaef_et = 0.0_dp
3088
3089 ! loop over domain - for applying power law later on
3090 do idomain = 1, domainmeta%nDomains
3091
3092 ! check whether the time dimension of SM and ET are same
3093 if (size(etoptisim(idomain)%dataSim, dim = 2) /= size(smoptisim(idomain)%dataSim, dim = 2)) then
3094 call message('***ERROR: mo_objective_funtion: objective_spaef_sm_spaef_et: The time dimension of SM (', &
3095 trim(num2str(size(smoptisim(idomain)%dataSim, dim = 2))),') and ET (', &
3096 trim(num2str(size(etoptisim(idomain)%dataSim, dim = 2))),') are not equal for domain ', &
3097 trim(num2str(idomain)), &
3098 '. Are both of them at same temporal resolution?')
3099 stop
3100 end if
3101
3102 ! allocate
3103 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
3104 allocate(spaef_time_series_sm(size(etoptisim(idomain)%dataSim, dim = 2)))
3105 allocate(spaef_time_series_et(size(etoptisim(idomain)%dataSim, dim = 2)))
3106
3107 allocate(vec1_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3108 allocate(vec2_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3109 allocate(mask_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3110
3111 allocate(vec1_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3112 allocate(vec2_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3113 allocate(mask_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3114
3115 ! initalize
3116 mask_times = .false.
3117 spaef_time_series_sm = 0.0_dp
3118 spaef_time_series_et = 0.0_dp
3119
3120 ! calculate the SPAEF criterion across each time step
3121 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
3122
3123 vec1_et = l1_etobs(idomain)%dataObs(:, itime)
3124 vec2_et = etoptisim(idomain)%dataSim(:, itime)
3125 mask_et = l1_etobs(idomain)%maskObs(:, itime)
3126
3127 vec1_sm = l1_smobs(idomain)%dataObs(:, itime)
3128 vec2_sm = smoptisim(idomain)%dataSim(:, itime)
3129 mask_sm = l1_smobs(idomain)%maskObs(:, itime)
3130
3131 if (count(mask_et) > 1_i4 .and. count(mask_sm) > 1_i4) then ! sample size must be at least 2
3132 mask_times(itime) = .true.
3133 spaef_time_series_sm(itime) = spaef(vec1_sm, vec2_sm, mask = mask_sm)
3134 spaef_time_series_et(itime) = spaef(vec1_et, vec2_et, mask = mask_et)
3135 end if
3136 end do
3137
3138 if (count(mask_times) > 0_i4) then
3139 ! calculate the temporal average of SPAEF, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
3140 objective_spaef_sm_spaef_et = objective_spaef_sm_spaef_et + &
3141 ( &
3142 (1.0_dp - ( &
3143 0.5 * sum(spaef_time_series_sm, mask = mask_times) + &
3144 0.5 * sum(spaef_time_series_et, mask = mask_times) &
3145 ) / real(count(mask_times), dp) &
3146 ) / real(domainmeta%overallNumberOfDomains, dp) &
3147 )**6
3148 call message(' spaef_sm = ', &
3149 & num2str(sum(spaef_time_series_sm, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
3150 &' for domain ', num2str(idomain))
3151 call message(' spaef_et = ', &
3152 & num2str(sum(spaef_time_series_et, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
3153 &' for domain ', num2str(idomain))
3154 else
3155 call message('***ERROR: mo_objective_funtion: objective_spaef_sm_spaef_et: No SM and/or ET observations available!')
3156 stop
3157 end if
3158
3159 ! deallocate
3160 deallocate(mask_times)
3161 deallocate(spaef_time_series_sm, spaef_time_series_et)
3162 deallocate(vec1_sm, vec1_et)
3163 deallocate(vec2_sm, vec2_et)
3164 deallocate(mask_et, mask_sm)
3165 call smoptisim(idomain)%destroy()
3166 call etoptisim(idomain)%destroy()
3167 end do
3168 deallocate(etoptisim)
3169 deallocate(smoptisim)
3170
3171#ifndef MPI
3172 objective_spaef_sm_spaef_et = objective_spaef_sm_spaef_et**onesixth
3173
3174 call message(' objective_spaef_sm_spaef_et = ', num2str(objective_spaef_sm_spaef_et, '(F9.5)'))
3175#endif
3176
3177 END FUNCTION objective_spaef_sm_spaef_et
3178
3179 ! ------------------------------------------------------------------
3180
3181 ! NAME
3182 ! objective_spaef_sm_kge_q
3183
3184 ! PURPOSE
3185 !> \brief Bivariate Objective function for SM and Q.
3186
3187 !> \details The Spatial Efficiency (SPAEF) of observed and modeled SM
3188 !> fields and KGE of streamflow (input wt:1 weights) are calculated - aim: matching spatial patterns
3189
3190 ! INTENT(IN)
3191 !> \param[in] "real(dp), dimension(:) :: parameterset"
3192 !> \param[in] "procedure(eval_interface) :: eval"
3193
3194 ! RETURN
3195 !> \return real(dp) :: objective_spaef_sm_kge_q &mdash; objective function value
3196 !> (which will be e.g. minimized by an optimization routine like DDS)
3197
3198 ! LITERATURE
3199 ! Koch, J., Demirel, M. C., & Stisen, S. (2018). The SPAtial EFficiency metric (SPAEF):
3200 ! Multiple-component evaluation of spatial patterns for optimization of hydrological models.
3201 ! Geoscientific Model Development, 11(5), 1873–1886. https://doi.org/10.5194/gmd-11-1873-2018
3202
3203
3204 ! HISTORY
3205 !> \authors Pallav Shrestha
3206
3207 !> \date Dec 2024
3208
3209 ! Modifications:
3210
3211
3212 FUNCTION objective_spaef_sm_kge_q(parameterset, eval)
3213
3214 use mo_optimization_types, only : optidata_sim
3215 use mo_common_constants, only : nodata_dp
3216 use mo_common_variables, only : level1, domainmeta
3217 use mo_global_variables, only : l1_smobs, wt_for_optional_data
3218 use mo_message, only : message
3219 use mo_errormeasures, only : spaef, kge
3220 use mo_string_utils, only : num2str
3221 use mo_mrm_global_variables, only : gauge
3222 use mo_mrm_objective_function_runoff, only : extract_runoff
3223
3224 implicit none
3225
3226 real(dp), dimension(:), intent(in) :: parameterset
3227
3228 procedure(eval_interface), INTENT(IN), POINTER :: eval
3229
3230 ! objective function value
3231 real(dp) :: objective_spaef_sm_kge_q
3232 real(dp) :: objective_kgestreamflow
3233
3234 ! domain loop counter
3235 integer(i4) :: idomain, jdomain
3236
3237 ! time loop counter
3238 integer(i4) :: itime
3239
3240 ! for sixth root
3241#ifndef MPI
3242 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
3243#endif
3244
3245 ! vectorized arrays of SM
3246 real(dp), dimension(:), allocatable :: vec1_sm, vec2_sm
3247
3248 ! spaef at every time step
3249 real(dp), dimension(:), allocatable :: spaef_time_series_sm
3250
3251 ! simulated streamflow
3252 ! dim1=nTimeSteps, dim2=nGauges
3253 real(dp), allocatable, dimension(:, :) :: runoff
3254
3255 ! simulated SM
3256 type(optidata_sim), dimension(:), allocatable :: smoptisim
3257
3258 ! mask of valid cells
3259 logical, dimension(:), allocatable :: mask_sm
3260
3261 ! mask for valid time steps
3262 logical, dimension(:), allocatable :: mask_times
3263
3264 ! gauges counter
3265 integer(i4) :: gg
3266 integer(i4) :: ngaugestotal
3267
3268 ! aggregated simulated
3269 real(dp), dimension(:), allocatable :: runoff_obs
3270
3271 ! measured
3272 real(dp), dimension(:), allocatable :: runoff_agg
3273
3274 ! mask for measured
3275 logical, dimension(:), allocatable :: runoff_obs_mask
3276
3277
3278
3279 allocate(smoptisim(domainmeta%nDomains))
3280 call eval(parameterset, smoptisim = smoptisim, runoff = runoff)
3281 ngaugestotal = size(runoff, dim = 2)
3282
3283 ! initialize some variables
3284 objective_spaef_sm_kge_q = 0.0_dp
3285
3286 ! loop over domain - for applying power law later on
3287 do idomain = 1, domainmeta%nDomains
3288
3289 ! allocate
3290 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
3291 allocate(spaef_time_series_sm(size(smoptisim(idomain)%dataSim, dim = 2)))
3292
3293 allocate(vec1_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3294 allocate(vec2_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3295 allocate(mask_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3296
3297 ! initalize
3298 mask_times = .false.
3299 spaef_time_series_sm = 0.0_dp
3300 objective_kgestreamflow = 0.0_dp
3301
3302 ! calculate the SPAEF criterion across each time step
3303 do itime = 1, size(smoptisim(idomain)%dataSim, dim = 2)
3304
3305 vec1_sm = l1_smobs(idomain)%dataObs(:, itime)
3306 vec2_sm = smoptisim(idomain)%dataSim(:, itime)
3307 mask_sm = l1_smobs(idomain)%maskObs(:, itime)
3308
3309 if (count(mask_sm) > 1_i4) then ! sample size must be at least 2
3310 mask_times(itime) = .true.
3311 spaef_time_series_sm(itime) = spaef(vec1_sm, vec2_sm, mask = mask_sm)
3312 end if
3313 end do
3314
3315 ! calculate the KGE of streamflow
3316 do gg = 1, ngaugestotal
3317 ! get domain
3318 jdomain = gauge%domainId(gg)
3319 ! only proceed if current domain
3320 if (idomain == jdomain) then
3321 ! extract runoff
3322 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
3323 ! kge
3324 objective_kgestreamflow = objective_kgestreamflow + &
3325 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
3326 end if
3327 end do
3328 objective_kgestreamflow = objective_kgestreamflow / real(ngaugestotal, dp)
3329
3330 ! Combine
3331 if (count(mask_times) > 0_i4) then
3332 ! calculate the temporal average of SPAEF, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
3333 objective_spaef_sm_kge_q = objective_spaef_sm_kge_q + &
3334 ( &
3335 (1.0_dp - ( &
3336 wt_for_optional_data * sum(spaef_time_series_sm, mask = mask_times) / real(count(mask_times), dp) + &
3337 1_i4 * objective_kgestreamflow) / (wt_for_optional_data + 1_i4) &
3338 ) / real(domainmeta%overallNumberOfDomains, dp) &
3339 )**6
3340 call message(' spaef_sm = ', &
3341 & num2str(sum(spaef_time_series_sm, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
3342 &' for domain ', num2str(idomain))
3343 call message(' kge_q = ', num2str(objective_kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
3344 else
3345 call message('***ERROR: mo_objective_funtion: objective_spaef_sm_kge_q: No SM observations available!')
3346 stop
3347 end if
3348
3349 ! deallocate
3350 deallocate(mask_times)
3351 deallocate(spaef_time_series_sm)
3352 deallocate(vec1_sm)
3353 deallocate(vec2_sm)
3354 deallocate(mask_sm)
3355 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
3356 call smoptisim(idomain)%destroy()
3357 end do
3358 deallocate(smoptisim)
3359
3360#ifndef MPI
3361 objective_spaef_sm_kge_q = objective_spaef_sm_kge_q**onesixth
3362
3363 call message(' objective_spaef_sm_kge_q = ', num2str(objective_spaef_sm_kge_q, '(F9.5)'))
3364#endif
3365
3366 END FUNCTION objective_spaef_sm_kge_q
3367
3368
3369
3370 ! ------------------------------------------------------------------
3371
3372 ! NAME
3373 ! objective_spaef_et_kge_q
3374
3375 ! PURPOSE
3376 !> \brief Bivariate Objective function for ET and Q.
3377
3378 !> \details The Spatial Efficiency (SPAEF) of observed and modeled ET
3379 !> fields and KGE of streamflow (5:1 weights) are calculated - aim: matching spatial patterns
3380
3381 ! INTENT(IN)
3382 !> \param[in] "real(dp), dimension(:) :: parameterset"
3383 !> \param[in] "procedure(eval_interface) :: eval"
3384
3385 ! RETURN
3386 !> \return real(dp) :: objective_spaef_et_kge_q &mdash; objective function value
3387 !> (which will be e.g. minimized by an optimization routine like DDS)
3388
3389 ! LITERATURE
3390 ! Koch, J., Demirel, M. C., & Stisen, S. (2018). The SPAtial EFficiency metric (SPAEF):
3391 ! Multiple-component evaluation of spatial patterns for optimization of hydrological models.
3392 ! Geoscientific Model Development, 11(5), 1873–1886. https://doi.org/10.5194/gmd-11-1873-2018
3393
3394
3395 ! HISTORY
3396 !> \authors Pallav Shrestha
3397
3398 !> \date Dec 2024
3399
3400 ! Modifications:
3401
3402
3403 FUNCTION objective_spaef_et_kge_q(parameterset, eval)
3404
3405 use mo_optimization_types, only : optidata_sim
3406 use mo_common_constants, only : nodata_dp
3407 use mo_common_variables, only : level1, domainmeta
3408 use mo_global_variables, only : l1_etobs, wt_for_optional_data
3409 use mo_message, only : message
3410 use mo_errormeasures, only : spaef, kge
3411 use mo_string_utils, only : num2str
3412 use mo_mrm_global_variables, only : gauge
3413 use mo_mrm_objective_function_runoff, only : extract_runoff
3414
3415 implicit none
3416
3417 real(dp), dimension(:), intent(in) :: parameterset
3418
3419 procedure(eval_interface), INTENT(IN), POINTER :: eval
3420
3421 ! objective function value
3422 real(dp) :: objective_spaef_et_kge_q
3423 real(dp) :: objective_kgestreamflow
3424
3425 ! domain loop counter
3426 integer(i4) :: idomain, jdomain
3427
3428 ! time loop counter
3429 integer(i4) :: itime
3430
3431 ! for sixth root
3432#ifndef MPI
3433 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
3434#endif
3435
3436 ! vectorized arrays of ET
3437 real(dp), dimension(:), allocatable :: vec1_et, vec2_et
3438
3439 ! spaef at every time step
3440 real(dp), dimension(:), allocatable :: spaef_time_series_et
3441
3442 ! simulated streamflow
3443 ! dim1=nTimeSteps, dim2=nGauges
3444 real(dp), allocatable, dimension(:, :) :: runoff
3445
3446 ! simulated ET
3447 type(optidata_sim), dimension(:), allocatable :: etoptisim
3448
3449 ! mask of valid cells
3450 logical, dimension(:), allocatable :: mask_et
3451
3452 ! mask for valid time steps
3453 logical, dimension(:), allocatable :: mask_times
3454
3455 ! gauges counter
3456 integer(i4) :: gg
3457 integer(i4) :: ngaugestotal, ngaugesdomain
3458
3459 ! aggregated simulated
3460 real(dp), dimension(:), allocatable :: runoff_obs
3461
3462 ! measured
3463 real(dp), dimension(:), allocatable :: runoff_agg
3464
3465 ! mask for measured
3466 logical, dimension(:), allocatable :: runoff_obs_mask
3467
3468
3469
3470 allocate(etoptisim(domainmeta%nDomains))
3471 call eval(parameterset, etoptisim = etoptisim, runoff = runoff)
3472 ngaugestotal = size(runoff, dim = 2)
3473
3474 ! initialize some variables
3475 objective_spaef_et_kge_q = 0.0_dp
3476
3477 ! loop over domain - for applying power law later on
3478 do idomain = 1, domainmeta%nDomains
3479
3480 ! allocate
3481 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
3482 allocate(spaef_time_series_et(size(etoptisim(idomain)%dataSim, dim = 2)))
3483
3484 allocate(vec1_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3485 allocate(vec2_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3486 allocate(mask_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3487
3488 ! initalize
3489 mask_times = .false.
3490 spaef_time_series_et = 0.0_dp
3491 objective_kgestreamflow = 0.0_dp
3492 ngaugesdomain = 0._i4
3493
3494 ! calculate the SPAEF criterion across each time step
3495 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
3496
3497 vec1_et = l1_etobs(idomain)%dataObs(:, itime)
3498 vec2_et = etoptisim(idomain)%dataSim(:, itime)
3499 mask_et = l1_etobs(idomain)%maskObs(:, itime)
3500
3501 if (count(mask_et) > 1_i4) then ! sample size must be at least 2
3502 mask_times(itime) = .true.
3503 spaef_time_series_et(itime) = spaef(vec1_et, vec2_et, mask = mask_et)
3504 end if
3505 end do
3506
3507 ! calculate the KGE of streamflow
3508 do gg = 1, ngaugestotal
3509 ! get domain
3510 jdomain = gauge%domainId(gg)
3511 ! only proceed if current domain
3512 if (idomain == jdomain) then
3513 ! count
3514 ngaugesdomain = ngaugesdomain + 1_i4
3515 ! extract runoff
3516 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
3517 ! kge
3518 objective_kgestreamflow = objective_kgestreamflow + &
3519 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
3520 end if
3521 end do
3522 objective_kgestreamflow = objective_kgestreamflow / real(ngaugesdomain, dp)
3523
3524 ! Combine
3525 if (count(mask_times) > 0_i4) then
3526 ! calculate the temporal average of SPAEF, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
3527 objective_spaef_et_kge_q = objective_spaef_et_kge_q + &
3528 ( &
3529 (1.0_dp - ( &
3530 wt_for_optional_data * sum(spaef_time_series_et, mask = mask_times) / real(count(mask_times), dp) + &
3531 1 * objective_kgestreamflow) / (wt_for_optional_data + 1_i4) &
3532 ) / real(domainmeta%overallNumberOfDomains, dp) &
3533 )**6
3534 call message(' spaef_et = ', &
3535 & num2str(sum(spaef_time_series_et, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
3536 &' for domain ', num2str(idomain))
3537 call message(' kge_q = ', num2str(objective_kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
3538 else
3539 call message('***ERROR: mo_objective_funtion: objective_spaef_et_kge_q: No ET observations available!')
3540 stop
3541 end if
3542
3543 ! deallocate
3544 deallocate(mask_times)
3545 deallocate(spaef_time_series_et)
3546 deallocate(vec1_et)
3547 deallocate(vec2_et)
3548 deallocate(mask_et)
3549 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
3550 call etoptisim(idomain)%destroy()
3551 end do
3552 deallocate(etoptisim)
3553
3554#ifndef MPI
3555 objective_spaef_et_kge_q = objective_spaef_et_kge_q**onesixth
3556
3557 call message(' objective_spaef_et_kge_q = ', num2str(objective_spaef_et_kge_q, '(F9.5)'))
3558#endif
3559
3560 END FUNCTION objective_spaef_et_kge_q
3561
3562
3563 ! ------------------------------------------------------------------
3564
3565 ! NAME
3566 ! objective_spaef_sm_spaef_et_kge_q
3567
3568 ! PURPOSE
3569 !> \brief Trivariate Objective function for SM, ET and Q.
3570
3571 !> \details The Spatial Efficiency (SPAEF) of observed and modeled ET
3572 !> fields and KGE of streamflow (5:1 weights) are calculated - aim: matching spatial patterns
3573
3574 ! INTENT(IN)
3575 !> \param[in] "real(dp), dimension(:) :: parameterset"
3576 !> \param[in] "procedure(eval_interface) :: eval"
3577
3578 ! RETURN
3579 !> \return real(dp) :: objective_spaef_sm_spaef_et_kge_q &mdash; objective function value
3580 !> (which will be e.g. minimized by an optimization routine like DDS)
3581
3582 ! LITERATURE
3583 ! Koch, J., Demirel, M. C., & Stisen, S. (2018). The SPAtial EFficiency metric (SPAEF):
3584 ! Multiple-component evaluation of spatial patterns for optimization of hydrological models.
3585 ! Geoscientific Model Development, 11(5), 1873–1886. https://doi.org/10.5194/gmd-11-1873-2018
3586
3587
3588 ! HISTORY
3589 !> \authors Pallav Shrestha
3590
3591 !> \date Dec 2024
3592
3593 ! Modifications:
3594
3595
3596 FUNCTION objective_spaef_sm_spaef_et_kge_q(parameterset, eval)
3597
3598 use mo_optimization_types, only : optidata_sim
3599 use mo_common_constants, only : nodata_dp
3600 use mo_common_variables, only : level1, domainmeta
3601 use mo_global_variables, only : l1_etobs, l1_smobs, wt_for_optional_data
3602 use mo_message, only : message
3603 use mo_errormeasures, only : spaef, kge
3604 use mo_string_utils, only : num2str
3605 use mo_mrm_global_variables, only : gauge
3606 use mo_mrm_objective_function_runoff, only : extract_runoff
3607
3608 implicit none
3609
3610 real(dp), dimension(:), intent(in) :: parameterset
3611
3612 procedure(eval_interface), INTENT(IN), POINTER :: eval
3613
3614 ! objective function value
3615 real(dp) :: objective_spaef_sm_spaef_et_kge_q
3616 real(dp) :: objective_kgestreamflow
3617
3618 ! domain loop counter
3619 integer(i4) :: idomain, jdomain
3620
3621 ! time loop counter
3622 integer(i4) :: itime
3623
3624 ! for sixth root
3625#ifndef MPI
3626 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
3627#endif
3628
3629 ! vectorized arrays of ET and SM
3630 real(dp), dimension(:), allocatable :: vec1_et, vec2_et, vec1_sm, vec2_sm
3631
3632 ! spaef at every time step
3633 real(dp), dimension(:), allocatable :: spaef_time_series_et, spaef_time_series_sm
3634
3635 ! simulated streamflow
3636 ! dim1=nTimeSteps, dim2=nGauges
3637 real(dp), allocatable, dimension(:, :) :: runoff
3638
3639 ! simulated ET and SM
3640 type(optidata_sim), dimension(:), allocatable :: etoptisim, smoptisim
3641
3642 ! mask of valid cells
3643 logical, dimension(:), allocatable :: mask_et, mask_sm
3644
3645 ! mask for valid time steps
3646 logical, dimension(:), allocatable :: mask_times
3647
3648 ! gauges counter
3649 integer(i4) :: gg
3650 integer(i4) :: ngaugestotal, ngaugesdomain
3651
3652 ! aggregated simulated
3653 real(dp), dimension(:), allocatable :: runoff_obs
3654
3655 ! measured
3656 real(dp), dimension(:), allocatable :: runoff_agg
3657
3658 ! mask for measured
3659 logical, dimension(:), allocatable :: runoff_obs_mask
3660
3661
3662
3663 allocate(etoptisim(domainmeta%nDomains))
3664 allocate(smoptisim(domainmeta%nDomains))
3665 call eval(parameterset, smoptisim = smoptisim, etoptisim = etoptisim, runoff = runoff)
3666 ngaugestotal = size(runoff, dim = 2)
3667
3668 ! initialize some variables
3669 objective_spaef_sm_spaef_et_kge_q = 0.0_dp
3670
3671 ! loop over domain - for applying power law later on
3672 do idomain = 1, domainmeta%nDomains
3673
3674 ! check whether the time dimension of SM and ET are same
3675 if (size(etoptisim(idomain)%dataSim, dim = 2) /= size(smoptisim(idomain)%dataSim, dim = 2)) then
3676 call message('***ERROR: mo_objective_funtion: objective_spaef_sm_spaef_et: The time dimension of SM (', &
3677 trim(num2str(size(smoptisim(idomain)%dataSim, dim = 2))),') and ET (', &
3678 trim(num2str(size(etoptisim(idomain)%dataSim, dim = 2))),') are not equal for domain ', &
3679 trim(num2str(idomain)), &
3680 '. Are both of them at same temporal resolution?')
3681 stop
3682 end if
3683
3684 ! allocate
3685 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
3686 allocate(spaef_time_series_et(size(etoptisim(idomain)%dataSim, dim = 2)))
3687 allocate(spaef_time_series_sm(size(smoptisim(idomain)%dataSim, dim = 2)))
3688
3689 allocate(vec1_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3690 allocate(vec2_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3691 allocate(mask_et(size(etoptisim(idomain)%dataSim, dim = 1)))
3692
3693 allocate(vec1_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3694 allocate(vec2_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3695 allocate(mask_sm(size(smoptisim(idomain)%dataSim, dim = 1)))
3696
3697 ! initalize
3698 mask_times = .false.
3699 spaef_time_series_sm = 0.0_dp
3700 spaef_time_series_et = 0.0_dp
3701 objective_kgestreamflow = 0.0_dp
3702 ngaugesdomain = 0._i4
3703
3704 ! calculate the SPAEF criterion across each time step
3705 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
3706
3707 vec1_sm = l1_smobs(idomain)%dataObs(:, itime)
3708 vec2_sm = smoptisim(idomain)%dataSim(:, itime)
3709 mask_sm = l1_smobs(idomain)%maskObs(:, itime)
3710
3711 vec1_et = l1_etobs(idomain)%dataObs(:, itime)
3712 vec2_et = etoptisim(idomain)%dataSim(:, itime)
3713 mask_et = l1_etobs(idomain)%maskObs(:, itime)
3714
3715 if (count(mask_sm) > 1_i4 .and. count(mask_et) > 1_i4) then ! sample size must be at least 2
3716 mask_times(itime) = .true.
3717 spaef_time_series_sm(itime) = spaef(vec1_sm, vec2_sm, mask = mask_sm)
3718 spaef_time_series_et(itime) = spaef(vec1_et, vec2_et, mask = mask_et)
3719 end if
3720 end do
3721
3722 ! calculate the KGE of streamflow
3723 do gg = 1, ngaugestotal
3724 ! get domain
3725 jdomain = gauge%domainId(gg)
3726 ! only proceed if current domain
3727 if (idomain == jdomain) then
3728 ! count
3729 ngaugesdomain = ngaugesdomain + 1_i4
3730 ! extract runoff
3731 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
3732 ! kge
3733 objective_kgestreamflow = objective_kgestreamflow + &
3734 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
3735 end if
3736 end do
3737 objective_kgestreamflow = objective_kgestreamflow / real(ngaugesdomain, dp)
3738
3739 ! Combine
3740 if (count(mask_times) > 0_i4) then
3741 ! calculate the temporal average of SPAEF, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
3742 objective_spaef_sm_spaef_et_kge_q = objective_spaef_sm_spaef_et_kge_q + &
3743 ( &
3744 (1.0_dp - ( &
3745 0.5 * wt_for_optional_data * sum(spaef_time_series_sm, mask = mask_times) / real(count(mask_times), dp) + &
3746 0.5 * wt_for_optional_data * sum(spaef_time_series_et, mask = mask_times) / real(count(mask_times), dp) + &
3747 1 * objective_kgestreamflow) / (wt_for_optional_data + 1_i4) &
3748 ) / real(domainmeta%overallNumberOfDomains, dp) &
3749 )**6
3750 call message(' spaef_sm = ', &
3751 & num2str(sum(spaef_time_series_sm, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
3752 &' for domain ', num2str(idomain))
3753 call message(' spaef_et = ', &
3754 & num2str(sum(spaef_time_series_et, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
3755 &' for domain ', num2str(idomain))
3756 call message(' kge_q = ', num2str(objective_kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
3757
3758 else
3759 call message('***ERROR: mo_objective_funtion: objective_spaef_sm_spaef_et_kge_q: No SM and/or ET observations available!')
3760 stop
3761 end if
3762
3763 ! deallocate
3764 deallocate(mask_times)
3765 deallocate(spaef_time_series_et, spaef_time_series_sm)
3766 deallocate(vec1_et, vec1_sm)
3767 deallocate(vec2_et, vec2_sm)
3768 deallocate(mask_et, mask_sm)
3769 deallocate(runoff_agg, runoff_obs, runoff_obs_mask)
3770 call smoptisim(idomain)%destroy()
3771 call etoptisim(idomain)%destroy()
3772 end do
3773 deallocate(smoptisim)
3774 deallocate(etoptisim)
3775
3776#ifndef MPI
3777 objective_spaef_sm_spaef_et_kge_q = objective_spaef_sm_spaef_et_kge_q**onesixth
3778
3779 call message(' objective_spaef_sm_spaef_et_kge_q = ', num2str(objective_spaef_sm_spaef_et_kge_q, '(F9.5)'))
3780#endif
3781
3782 END FUNCTION objective_spaef_sm_spaef_et_kge_q
3783
3784
3785 ! ------------------------------------------------------------------
3786
3787 ! NAME
3788 ! objective_pd_et
3789
3790 ! PURPOSE
3791 !> \brief Objective function for ET.
3792
3793 !> \details The objective function only depends on a parameter vector.
3794 !> The model will be called with that parameter vector and
3795 !> the model output is subsequently compared to observed data.
3796
3797 !> Therefore the Pattern Dissimilarity (PD) of observed and modeled ET
3798 !> fields is calculated - aim: matching spatial patters
3799 !> \f[ E(t) = PD\left( ET_{obs}(t), ET_{sim}(t) \right) \f]
3800 !> where
3801 !> \f$ PD \f$ = pattern dissimilarity function,
3802 !> \f$ ET_{obs} \f$ = observed ET,
3803 !> \f$ ET_{sim} \f$ = simulated ET.
3804 !> \f$ E(t) \f$ = pattern dissimilarity at timestep \f$ t \f$.
3805 !> The the pattern dissimilaity (E) is spatially averaged as
3806 !> \f[ \phi_{i} = \frac{1}{T} \cdot \sum_{t=1}^T E_t \f]
3807 !> where \f$ T \f$ denotes the number of time steps.
3808 !> Finally, the overall objective function value \f$ OF \f$ is estimated based on the power-6
3809 !> norm to combine the \f$ \phi_{i} \f$ from all domains \f$ N \f$.
3810 !> \f[ OF = \sqrt[6]{\sum((1.0 - \phi_{i})/N)^6 } . \f]
3811 !> The observed data L1_et, L1_et_mask are global in this module.
3812 !> The observed data L1_et, L1_et_mask are global in this module.
3813
3814 ! INTENT(IN)
3815 !> \param[in] "real(dp), dimension(:) :: parameterset"
3816 !> \param[in] "procedure(eval_interface) :: eval"
3817
3818 ! RETURN
3819 !> \return real(dp) :: objecive_sm_pd &mdash; objective function value
3820 !> (which will be e.g. minimized by an optimization routine like DDS)
3821
3822 ! HISTORY
3823 !> \authors Matthias Zink
3824
3825 !> \date May 2015
3826
3827 ! Modifications:
3828 ! Robert Schweppe Jun 2018 - refactoring and reformatting
3829 ! Pallav Shrestha Jan 2025 - adapting for ET
3830
3831 FUNCTION objective_pd_et(parameterset, eval)
3832
3833 use mo_optimization_types, only : optidata_sim
3834 use mo_common_constants, only : nodata_dp
3835 use mo_common_variables, only : level1, domainmeta
3836 use mo_global_variables, only : l1_etobs
3837 use mo_message, only : message
3838 use mo_spatialsimilarity, only : pd
3839 use mo_string_utils, only : num2str
3840
3841 implicit none
3842
3843 real(dp), dimension(:), intent(in) :: parameterset
3844
3845 procedure(eval_interface), INTENT(IN), POINTER :: eval
3846
3847 ! objective function value
3848 real(dp) :: objective_pd_et
3849
3850 ! domain loop counter
3851 integer(i4) :: idomain
3852
3853 ! time loop counter
3854 integer(i4) :: itime
3855
3856 ! level 1 number of columns and rows
3857 integer(i4) :: nrows1, ncols1
3858
3859 ! for sixth root
3860#ifndef MPI
3861 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
3862#endif
3863
3864 ! matrices of SM from vectorized arrays
3865 real(dp), dimension(:, :), allocatable :: mat1, mat2
3866
3867 ! pattern dissimilarity (pd) at every time step
3868 real(dp), dimension(:), allocatable :: pd_time_series
3869
3870 ! simulated ET
3871 type(optidata_sim), dimension(:), allocatable :: etoptisim
3872
3873 ! mask of valid cells at level1
3874 logical, dimension(:, :), allocatable :: mask1
3875
3876 ! mask of valid sm cells
3877 logical, dimension(:, :), allocatable :: mask_et
3878
3879 ! mask for valid sm catchment avg time steps
3880 logical, dimension(:), allocatable :: mask_times
3881
3882
3883 allocate(etoptisim(domainmeta%nDomains))
3884 call eval(parameterset, etoptisim = etoptisim)
3885
3886 ! initialize some variables
3887 objective_pd_et = 0.0_dp
3888
3889 ! loop over domain - for applying power law later on
3890 do idomain = 1, domainmeta%nDomains
3891
3892 ! get domain information
3893 mask1 = level1(idomain)%mask
3894 ncols1 = level1(idomain)%ncols
3895 nrows1 = level1(idomain)%nrows
3896
3897 ! allocate
3898 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
3899 allocate(pd_time_series(size(etoptisim(idomain)%dataSim, dim = 2)))
3900 allocate(mat1(nrows1, ncols1))
3901 allocate(mat2(nrows1, ncols1))
3902 allocate(mask_et(nrows1, ncols1))
3903
3904 ! initalize
3905 mask_times = .false.
3906 pd_time_series = 0.0_dp
3907
3908 ! calculate pattern similarity criterion
3909 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
3910 mat1 = unpack(l1_etobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
3911 mat2 = unpack(etoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
3912 mask_et = unpack(l1_etobs(idomain)%maskObs(:, itime), mask1, .false.)
3913 pd_time_series = pd(mat1, mat2, mask = mask_et, valid = mask_times(itime))
3914 end do
3915
3916 if (count(mask_times) > 0_i4) then
3917 ! calculate avergae PD over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
3918 objective_pd_et = objective_pd_et + &
3919 ((1.0_dp - sum(pd_time_series, mask = mask_times) / real(count(mask_times), dp)) / &
3920 real(domainmeta%overallnumberofdomains, dp))**6
3921 else
3922 call message('***ERROR: mo_objective_funtion: objective_pd_et: No ET observations available!')
3923 stop
3924 end if
3925
3926 ! deallocate
3927 deallocate(mask_times)
3928 deallocate(pd_time_series)
3929 deallocate(mat1)
3930 deallocate(mat2)
3931 deallocate(mask_et)
3932 call etoptisim(idomain)%destroy()
3933 end do
3934 deallocate(etoptisim)
3935
3936#ifndef MPI
3937 objective_pd_et = objective_pd_et**onesixth
3938
3939 call message(' objective_pd_et = ', num2str(objective_pd_et, '(F9.5)'))
3940 call message(' pd_et = ', num2str(1 - objective_pd_et, '(F9.5)'))
3941#endif
3942
3943 END FUNCTION objective_pd_et
3944
3945
3946 ! ------------------------------------------------------------------
3947
3948 ! NAME
3949 ! objective_pd_sm_pd_et
3950
3951 ! PURPOSE
3952 !> \brief Bivariate Objective function for SM and ET.
3953
3954 !> \details The objective function only depends on a parameter vector.
3955 !> The model will be called with that parameter vector and
3956 !> the model output is subsequently compared to observed data.
3957
3958 ! INTENT(IN)
3959 !> \param[in] "real(dp), dimension(:) :: parameterset"
3960 !> \param[in] "procedure(eval_interface) :: eval"
3961
3962 ! RETURN
3963 !> \return real(dp) :: objecive_sm_pd &mdash; objective function value
3964 !> (which will be e.g. minimized by an optimization routine like DDS)
3965
3966 ! HISTORY
3967 !> \authors Matthias Zink
3968
3969 !> \date May 2015
3970
3971 ! Modifications:
3972 ! Robert Schweppe Jun 2018 - refactoring and reformatting
3973 ! Pallav Shrestha Jan 2025 - adapting for ET + SM
3974
3975 FUNCTION objective_pd_sm_pd_et(parameterset, eval)
3976
3977 use mo_optimization_types, only : optidata_sim
3978 use mo_common_constants, only : nodata_dp
3979 use mo_common_variables, only : level1, domainmeta
3980 use mo_global_variables, only : l1_smobs, l1_etobs
3981 use mo_message, only : message
3982 use mo_spatialsimilarity, only : pd
3983 use mo_string_utils, only : num2str
3984
3985 implicit none
3986
3987 real(dp), dimension(:), intent(in) :: parameterset
3988
3989 procedure(eval_interface), INTENT(IN), POINTER :: eval
3990
3991 ! objective function value
3992 real(dp) :: objective_pd_sm_pd_et
3993
3994 ! domain loop counter
3995 integer(i4) :: idomain
3996
3997 ! time loop counter
3998 integer(i4) :: itime
3999
4000 ! level 1 number of columns and rows
4001 integer(i4) :: nrows1, ncols1
4002
4003 ! for sixth root
4004#ifndef MPI
4005 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
4006#endif
4007
4008 ! matrices from vectorized arrays
4009 real(dp), dimension(:, :), allocatable :: mat1_sm, mat2_sm, mat1_et, mat2_et
4010
4011 ! pattern dissimilarity (pd) at every time step
4012 real(dp), dimension(:), allocatable :: pd_time_series_sm, pd_time_series_et
4013
4014 ! simulated SM
4015 type(optidata_sim), dimension(:), allocatable :: smoptisim
4016
4017 ! simulated ET
4018 type(optidata_sim), dimension(:), allocatable :: etoptisim
4019
4020 ! mask of valid cells at level1
4021 logical, dimension(:, :), allocatable :: mask1
4022
4023 ! mask of valid cells
4024 logical, dimension(:, :), allocatable :: mask_et, mask_sm
4025
4026 ! mask for valid catchment avg time steps
4027 logical, dimension(:), allocatable :: mask_times
4028
4029
4030 allocate(smoptisim(domainmeta%nDomains))
4031 allocate(etoptisim(domainmeta%nDomains))
4032 call eval(parameterset, smoptisim = smoptisim, etoptisim = etoptisim)
4033
4034 ! initialize some variables
4035 objective_pd_sm_pd_et = 0.0_dp
4036
4037 ! loop over domain - for applying power law later on
4038 do idomain = 1, domainmeta%nDomains
4039
4040 ! check whether the time dimension of SM and ET are same
4041 if (size(etoptisim(idomain)%dataSim, dim = 2) /= size(smoptisim(idomain)%dataSim, dim = 2)) then
4042 call message('***ERROR: mo_objective_funtion: objective_pd_sm_pd_et: The time dimension of SM (', &
4043 trim(num2str(size(smoptisim(idomain)%dataSim, dim = 2))),') and ET (', &
4044 trim(num2str(size(etoptisim(idomain)%dataSim, dim = 2))),') are not equal for domain ', &
4045 trim(num2str(idomain)), &
4046 '. Are both of them at same temporal resolution?')
4047 stop
4048 end if
4049
4050 ! get domain information
4051 mask1 = level1(idomain)%mask
4052 ncols1 = level1(idomain)%ncols
4053 nrows1 = level1(idomain)%nrows
4054
4055 ! allocate
4056 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
4057 allocate(pd_time_series_sm(size(smoptisim(idomain)%dataSim, dim = 2)))
4058 allocate(pd_time_series_et(size(etoptisim(idomain)%dataSim, dim = 2)))
4059
4060 allocate(mat1_sm(nrows1, ncols1))
4061 allocate(mat2_sm(nrows1, ncols1))
4062 allocate(mask_sm(nrows1, ncols1))
4063
4064 allocate(mat1_et(nrows1, ncols1))
4065 allocate(mat2_et(nrows1, ncols1))
4066 allocate(mask_et(nrows1, ncols1))
4067
4068 ! initalize
4069 mask_times = .false.
4070 pd_time_series_sm = 0.0_dp
4071 pd_time_series_et = 0.0_dp
4072
4073 ! calculate pattern similarity criterion
4074 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
4075
4076 mat1_sm = unpack(l1_smobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
4077 mat2_sm = unpack(smoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
4078 mask_sm = unpack(l1_smobs(idomain)%maskObs(:, itime), mask1, .false.)
4079 pd_time_series_sm = pd(mat1_sm, mat2_sm, mask = mask_sm, valid = mask_times(itime))
4080
4081 mat1_et = unpack(l1_etobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
4082 mat2_et = unpack(etoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
4083 mask_et = unpack(l1_etobs(idomain)%maskObs(:, itime), mask1, .false.)
4084 pd_time_series_et = pd(mat1_et, mat2_et, mask = mask_et, valid = mask_times(itime))
4085
4086 end do
4087
4088 if (count(mask_times) > 0_i4) then
4089 ! calculate avergae PD over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
4090 objective_pd_sm_pd_et = objective_pd_sm_pd_et + &
4091 ( &
4092 (1.0_dp - ( &
4093 0.5 * sum(pd_time_series_sm, mask = mask_times) + &
4094 0.5 * sum(pd_time_series_et, mask = mask_times) &
4095 ) / real(count(mask_times), dp) &
4096 ) / real(domainmeta%overallNumberOfDomains, dp) &
4097 )**6
4098 call message(' pd_sm = ', &
4099 & num2str(sum(pd_time_series_sm, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
4100 &' for domain ', num2str(idomain))
4101 call message(' pd_et = ', &
4102 & num2str(sum(pd_time_series_et, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
4103 &' for domain ', num2str(idomain))
4104 else
4105 call message('***ERROR: mo_objective_funtion: objective_pd_sm_pd_et: No SM and/or ET observations available!')
4106 stop
4107 end if
4108
4109 ! deallocate
4110 deallocate(mask_times)
4111 deallocate(pd_time_series_sm, pd_time_series_et)
4112 deallocate(mat1_sm, mat1_et)
4113 deallocate(mat2_sm, mat2_et)
4114 deallocate(mask_sm, mask_et)
4115 call etoptisim(idomain)%destroy()
4116 end do
4117 deallocate(smoptisim)
4118 deallocate(etoptisim)
4119
4120#ifndef MPI
4121 objective_pd_sm_pd_et = objective_pd_sm_pd_et**onesixth
4122
4123 call message(' objective_pd_sm_pd_et = ', num2str(objective_pd_sm_pd_et, '(F9.5)'))
4124 call message(' pd_et = ', num2str(1 - objective_pd_sm_pd_et, '(F9.5)'))
4125#endif
4126
4127 END FUNCTION objective_pd_sm_pd_et
4128
4129
4130
4131 ! ------------------------------------------------------------------
4132
4133 ! NAME
4134 ! objective_pd_sm_kge_q
4135
4136 ! PURPOSE
4137 !> \brief Bivariate Objective function for SM and Q.
4138
4139 !> \details The Pattern Dissimilarity (PD) of observed and modeled SM
4140 !> fields and KGE of streamflow (input wt:1 weights) are calculated - aim: matching spatial patterns
4141
4142 ! INTENT(IN)
4143 !> \param[in] "real(dp), dimension(:) :: parameterset"
4144 !> \param[in] "procedure(eval_interface) :: eval"
4145
4146 ! RETURN
4147 !> \return real(dp) :: objecive_sm_pd &mdash; objective function value
4148 !> (which will be e.g. minimized by an optimization routine like DDS)
4149
4150 ! HISTORY
4151 !> \authors Matthias Zink
4152
4153 !> \date May 2015
4154
4155 ! Modifications:
4156 ! Robert Schweppe Jun 2018 - refactoring and reformatting
4157 ! Pallav Shrestha Jan 2025 - adapting for PD of SM and KGE of Q
4158
4159 FUNCTION objective_pd_sm_kge_q(parameterset, eval)
4160
4161 use mo_optimization_types, only : optidata_sim
4162 use mo_common_constants, only : nodata_dp
4163 use mo_common_variables, only : level1, domainmeta
4164 use mo_global_variables, only : l1_smobs, wt_for_optional_data
4165 use mo_message, only : message
4166 use mo_errormeasures, only : kge
4167 use mo_spatialsimilarity, only : pd
4168 use mo_string_utils, only : num2str
4169 use mo_mrm_global_variables, only : gauge
4170 use mo_mrm_objective_function_runoff, only : extract_runoff
4171
4172 implicit none
4173
4174 real(dp), dimension(:), intent(in) :: parameterset
4175
4176 procedure(eval_interface), INTENT(IN), POINTER :: eval
4177
4178 ! objective function value
4179 real(dp) :: objective_pd_sm_kge_q
4180 real(dp) :: objective_kgestreamflow
4181
4182 ! domain loop counter
4183 integer(i4) :: idomain, jdomain
4184
4185 ! time loop counter
4186 integer(i4) :: itime
4187
4188 ! level 1 number of columns and rows
4189 integer(i4) :: nrows1, ncols1
4190
4191 ! for sixth root
4192#ifndef MPI
4193 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
4194#endif
4195
4196 ! matrices of SM from vectorized arrays
4197 real(dp), dimension(:, :), allocatable :: mat1_sm, mat2_sm
4198
4199 ! pattern dissimilarity (pd) at every time step
4200 real(dp), dimension(:), allocatable :: pd_time_series_sm
4201
4202 ! simulated streamflow
4203 ! dim1=nTimeSteps, dim2=nGauges
4204 real(dp), allocatable, dimension(:, :) :: runoff
4205
4206 ! simulated soil moisture
4207 type(optidata_sim), dimension(:), allocatable :: smoptisim
4208
4209 ! mask of valid cells at level1
4210 logical, dimension(:, :), allocatable :: mask1
4211
4212 ! mask of valid sm cells
4213 logical, dimension(:, :), allocatable :: mask_sm
4214
4215 ! mask for valid sm catchment avg time steps
4216 logical, dimension(:), allocatable :: mask_times
4217
4218 ! gauges counter
4219 integer(i4) :: gg
4220 integer(i4) :: ngaugestotal, ngaugesdomain
4221
4222 ! aggregated simulated
4223 real(dp), dimension(:), allocatable :: runoff_obs
4224
4225 ! measured
4226 real(dp), dimension(:), allocatable :: runoff_agg
4227
4228 ! mask for measured
4229 logical, dimension(:), allocatable :: runoff_obs_mask
4230
4231
4232 allocate(smoptisim(domainmeta%nDomains))
4233 call eval(parameterset, smoptisim = smoptisim, runoff = runoff)
4234 ngaugestotal = size(runoff, dim = 2)
4235
4236 ! initialize some variables
4237 objective_pd_sm_kge_q = 0.0_dp
4238
4239 ! loop over domain - for applying power law later on
4240 do idomain = 1, domainmeta%nDomains
4241
4242 ! get domain information
4243 mask1 = level1(idomain)%mask
4244 ncols1 = level1(idomain)%ncols
4245 nrows1 = level1(idomain)%nrows
4246
4247 ! allocate
4248 allocate(mask_times(size(smoptisim(idomain)%dataSim, dim = 2)))
4249 allocate(pd_time_series_sm(size(smoptisim(idomain)%dataSim, dim = 2)))
4250 allocate(mat1_sm(nrows1, ncols1))
4251 allocate(mat2_sm(nrows1, ncols1))
4252 allocate(mask_sm(nrows1, ncols1))
4253
4254 ! initalize
4255 mask_times = .false.
4256 pd_time_series_sm = 0.0_dp
4257 objective_kgestreamflow = 0.0_dp
4258 ngaugesdomain = 0._i4
4259
4260 ! calculate pattern similarity criterion
4261 do itime = 1, size(smoptisim(idomain)%dataSim, dim = 2)
4262
4263 mat1_sm = unpack(l1_smobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
4264 mat2_sm = unpack(smoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
4265 mask_sm = unpack(l1_smobs(idomain)%maskObs(:, itime), mask1, .false.)
4266 pd_time_series_sm = pd(mat1_sm, mat2_sm, mask = mask_sm, valid = mask_times(itime))
4267
4268 end do
4269
4270 ! calculate the KGE of streamflow
4271 do gg = 1, ngaugestotal
4272 ! get domain
4273 jdomain = gauge%domainId(gg)
4274 ! only proceed if current domain
4275 if (idomain == jdomain) then
4276 ! count
4277 ngaugesdomain = ngaugesdomain + 1_i4
4278 ! extract runoff
4279 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
4280 ! kge
4281 objective_kgestreamflow = objective_kgestreamflow + &
4282 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
4283 end if
4284 end do
4285 objective_kgestreamflow = objective_kgestreamflow / real(ngaugesdomain, dp)
4286
4287 ! Combine
4288 if (count(mask_times) > 0_i4) then
4289 ! calculate the temporal average of PD, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
4290 objective_pd_sm_kge_q = objective_pd_sm_kge_q + &
4291 ( &
4292 (1.0_dp - ( &
4293 wt_for_optional_data * sum(pd_time_series_sm, mask = mask_times) / real(count(mask_times), dp) + &
4294 1_i4 * objective_kgestreamflow) / (wt_for_optional_data + 1_i4) &
4295 ) / real(domainmeta%overallNumberOfDomains, dp) &
4296 )**6
4297 call message(' pd_sm = ', &
4298 & num2str(sum(pd_time_series_sm, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
4299 &' for domain ', num2str(idomain))
4300 call message(' kge_q = ', num2str(objective_kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
4301 else
4302 call message('***ERROR: mo_objective_funtion: objective_pd_sm_kge_q: No SM observations available!')
4303 stop
4304 end if
4305
4306 ! deallocate
4307 deallocate(mask_times)
4308 deallocate(pd_time_series_sm)
4309 deallocate(mat1_sm)
4310 deallocate(mat2_sm)
4311 deallocate(mask_sm)
4312 call smoptisim(idomain)%destroy()
4313 end do
4314 deallocate(smoptisim)
4315
4316#ifndef MPI
4317 objective_pd_sm_kge_q = objective_pd_sm_kge_q**onesixth
4318
4319 call message(' objective_pd_sm_kge_q = ', num2str(objective_pd_sm_kge_q, '(F9.5)'))
4320#endif
4321
4322 END FUNCTION objective_pd_sm_kge_q
4323
4324
4325
4326 ! ------------------------------------------------------------------
4327
4328 ! NAME
4329 ! objective_pd_et_kge_q
4330
4331 ! PURPOSE
4332 !> \brief Bivariate Objective function for ET and Q.
4333
4334 !> \details The Pattern Dissimilarity (PD) of observed and modeled ET
4335 !> fields and KGE of streamflow (input wt:1 weights) are calculated - aim: matching spatial patterns
4336
4337 ! INTENT(IN)
4338 !> \param[in] "real(dp), dimension(:) :: parameterset"
4339 !> \param[in] "procedure(eval_interface) :: eval"
4340
4341 ! RETURN
4342 !> \return real(dp) :: objecive_sm_pd &mdash; objective function value
4343 !> (which will be e.g. minimized by an optimization routine like DDS)
4344
4345 ! HISTORY
4346 !> \authors Matthias Zink
4347
4348 !> \date May 2015
4349
4350 ! Modifications:
4351 ! Robert Schweppe Jun 2018 - refactoring and reformatting
4352 ! Pallav Shrestha Jan 2025 - adapting for PD of ET and KGE of Q
4353
4354 FUNCTION objective_pd_et_kge_q(parameterset, eval)
4355
4356 use mo_optimization_types, only : optidata_sim
4357 use mo_common_constants, only : nodata_dp
4358 use mo_common_variables, only : level1, domainmeta
4359 use mo_global_variables, only : l1_etobs, wt_for_optional_data
4360 use mo_message, only : message
4361 use mo_errormeasures, only : kge
4362 use mo_spatialsimilarity, only : pd
4363 use mo_string_utils, only : num2str
4364 use mo_mrm_global_variables, only : gauge
4365 use mo_mrm_objective_function_runoff, only : extract_runoff
4366
4367 implicit none
4368
4369 real(dp), dimension(:), intent(in) :: parameterset
4370
4371 procedure(eval_interface), INTENT(IN), POINTER :: eval
4372
4373 ! objective function value
4374 real(dp) :: objective_pd_et_kge_q
4375 real(dp) :: objective_kgestreamflow
4376
4377 ! domain loop counter
4378 integer(i4) :: idomain, jdomain
4379
4380 ! time loop counter
4381 integer(i4) :: itime
4382
4383 ! level 1 number of columns and rows
4384 integer(i4) :: nrows1, ncols1
4385
4386 ! for sixth root
4387#ifndef MPI
4388 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
4389#endif
4390
4391 ! matrices of SM from vectorized arrays
4392 real(dp), dimension(:, :), allocatable :: mat1_et, mat2_et
4393
4394 ! pattern dissimilarity (pd) at every time step
4395 real(dp), dimension(:), allocatable :: pd_time_series_et
4396
4397 ! simulated streamflow
4398 ! dim1=nTimeSteps, dim2=nGauges
4399 real(dp), allocatable, dimension(:, :) :: runoff
4400
4401 ! simulated soil moisture
4402 type(optidata_sim), dimension(:), allocatable :: etoptisim
4403
4404 ! mask of valid cells at level1
4405 logical, dimension(:, :), allocatable :: mask1
4406
4407 ! mask of valid sm cells
4408 logical, dimension(:, :), allocatable :: mask_et
4409
4410 ! mask for valid sm catchment avg time steps
4411 logical, dimension(:), allocatable :: mask_times
4412
4413 ! gauges counter
4414 integer(i4) :: gg
4415 integer(i4) :: ngaugestotal, ngaugesdomain
4416
4417 ! aggregated simulated
4418 real(dp), dimension(:), allocatable :: runoff_obs
4419
4420 ! measured
4421 real(dp), dimension(:), allocatable :: runoff_agg
4422
4423 ! mask for measured
4424 logical, dimension(:), allocatable :: runoff_obs_mask
4425
4426
4427 allocate(etoptisim(domainmeta%nDomains))
4428 call eval(parameterset, etoptisim = etoptisim, runoff = runoff)
4429 ngaugestotal = size(runoff, dim = 2)
4430
4431 ! initialize some variables
4432 objective_pd_et_kge_q = 0.0_dp
4433
4434 ! loop over domain - for applying power law later on
4435 do idomain = 1, domainmeta%nDomains
4436
4437 ! get domain information
4438 mask1 = level1(idomain)%mask
4439 ncols1 = level1(idomain)%ncols
4440 nrows1 = level1(idomain)%nrows
4441
4442 ! allocate
4443 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
4444 allocate(pd_time_series_et(size(etoptisim(idomain)%dataSim, dim = 2)))
4445 allocate(mat1_et(nrows1, ncols1))
4446 allocate(mat2_et(nrows1, ncols1))
4447 allocate(mask_et(nrows1, ncols1))
4448
4449 ! initalize
4450 mask_times = .false.
4451 pd_time_series_et = 0.0_dp
4452 objective_kgestreamflow = 0.0_dp
4453 ngaugesdomain = 0._i4
4454
4455 ! calculate pattern similarity criterion
4456 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
4457
4458 mat1_et = unpack(l1_etobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
4459 mat2_et = unpack(etoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
4460 mask_et = unpack(l1_etobs(idomain)%maskObs(:, itime), mask1, .false.)
4461 pd_time_series_et = pd(mat1_et, mat2_et, mask = mask_et, valid = mask_times(itime))
4462
4463 end do
4464
4465 ! calculate the KGE of streamflow
4466 do gg = 1, ngaugestotal
4467 ! get domain
4468 jdomain = gauge%domainId(gg)
4469 ! only proceed if current domain
4470 if (idomain == jdomain) then
4471 ! count
4472 ngaugesdomain = ngaugesdomain + 1_i4
4473 ! extract runoff
4474 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
4475 ! kge
4476 objective_kgestreamflow = objective_kgestreamflow + &
4477 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
4478 end if
4479 end do
4480 objective_kgestreamflow = objective_kgestreamflow / real(ngaugesdomain, dp)
4481
4482 ! Combine
4483 if (count(mask_times) > 0_i4) then
4484 ! calculate the temporal average of PD, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
4485 objective_pd_et_kge_q = objective_pd_et_kge_q + &
4486 ( &
4487 (1.0_dp - ( &
4488 wt_for_optional_data * sum(pd_time_series_et, mask = mask_times) / real(count(mask_times), dp) + &
4489 1_i4 * objective_kgestreamflow) / (wt_for_optional_data + 1_i4) &
4490 ) / real(domainmeta%overallNumberOfDomains, dp) &
4491 )**6
4492 call message(' pd_et = ', &
4493 & num2str(sum(pd_time_series_et, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
4494 &' for domain ', num2str(idomain))
4495 call message(' kge_q = ', num2str(objective_kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
4496 else
4497 call message('***ERROR: mo_objective_funtion: objective_pd_et_kge_q: No ET observations available!')
4498 stop
4499 end if
4500
4501 ! deallocate
4502 deallocate(mask_times)
4503 deallocate(pd_time_series_et)
4504 deallocate(mat1_et)
4505 deallocate(mat2_et)
4506 deallocate(mask_et)
4507 call etoptisim(idomain)%destroy()
4508 end do
4509 deallocate(etoptisim)
4510
4511#ifndef MPI
4512 objective_pd_et_kge_q = objective_pd_et_kge_q**onesixth
4513
4514 call message(' objective_pd_et_kge_q = ', num2str(objective_pd_et_kge_q, '(F9.5)'))
4515#endif
4516
4517 END FUNCTION objective_pd_et_kge_q
4518
4519
4520
4521
4522
4523 ! ------------------------------------------------------------------
4524
4525 ! NAME
4526 ! objective_pd_sm_pd_et_kge_q
4527
4528 ! PURPOSE
4529 !> \brief Trivariate Objective function for SM, ET and Q.
4530
4531 !> \details The Pattern Dissimilarity (PD) of observed and modeled ET and SM
4532 !> fields and KGE of streamflow (input wt:1 weights) are calculated - aim: matching spatial patterns
4533
4534 ! INTENT(IN)
4535 !> \param[in] "real(dp), dimension(:) :: parameterset"
4536 !> \param[in] "procedure(eval_interface) :: eval"
4537
4538 ! RETURN
4539 !> \return real(dp) :: objecive_sm_pd &mdash; objective function value
4540 !> (which will be e.g. minimized by an optimization routine like DDS)
4541
4542 ! HISTORY
4543 !> \authors Matthias Zink
4544
4545 !> \date May 2015
4546
4547 ! Modifications:
4548 ! Robert Schweppe Jun 2018 - refactoring and reformatting
4549 ! Pallav Shrestha Jan 2025 - adapting for PD of ET, SM and KGE of Q
4550
4551 FUNCTION objective_pd_sm_pd_et_kge_q(parameterset, eval)
4552
4553 use mo_optimization_types, only : optidata_sim
4554 use mo_common_constants, only : nodata_dp
4555 use mo_common_variables, only : level1, domainmeta
4556 use mo_global_variables, only : l1_etobs, l1_smobs, wt_for_optional_data
4557 use mo_message, only : message
4558 use mo_errormeasures, only : kge
4559 use mo_spatialsimilarity, only : pd
4560 use mo_string_utils, only : num2str
4561 use mo_mrm_global_variables, only : gauge
4562 use mo_mrm_objective_function_runoff, only : extract_runoff
4563
4564 implicit none
4565
4566 real(dp), dimension(:), intent(in) :: parameterset
4567
4568 procedure(eval_interface), INTENT(IN), POINTER :: eval
4569
4570 ! objective function value
4571 real(dp) :: objective_pd_sm_pd_et_kge_q
4572 real(dp) :: objective_kgestreamflow
4573
4574 ! domain loop counter
4575 integer(i4) :: idomain, jdomain
4576
4577 ! time loop counter
4578 integer(i4) :: itime
4579
4580 ! level 1 number of columns and rows
4581 integer(i4) :: nrows1, ncols1
4582
4583 ! for sixth root
4584#ifndef MPI
4585 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
4586#endif
4587
4588 ! matrices of SM from vectorized arrays
4589 real(dp), dimension(:, :), allocatable :: mat1_sm, mat2_sm, mat1_et, mat2_et
4590
4591 ! pattern dissimilarity (pd) at every time step
4592 real(dp), dimension(:), allocatable :: pd_time_series_sm, pd_time_series_et
4593
4594 ! simulated streamflow
4595 ! dim1=nTimeSteps, dim2=nGauges
4596 real(dp), allocatable, dimension(:, :) :: runoff
4597
4598 ! simulated soil moisture
4599 type(optidata_sim), dimension(:), allocatable :: etoptisim, smoptisim
4600
4601 ! mask of valid cells at level1
4602 logical, dimension(:, :), allocatable :: mask1
4603
4604 ! mask of valid sm cells
4605 logical, dimension(:, :), allocatable :: mask_sm, mask_et
4606
4607 ! mask for valid sm catchment avg time steps
4608 logical, dimension(:), allocatable :: mask_times
4609
4610 ! gauges counter
4611 integer(i4) :: gg
4612 integer(i4) :: ngaugestotal, ngaugesdomain
4613
4614 ! aggregated simulated
4615 real(dp), dimension(:), allocatable :: runoff_obs
4616
4617 ! measured
4618 real(dp), dimension(:), allocatable :: runoff_agg
4619
4620 ! mask for measured
4621 logical, dimension(:), allocatable :: runoff_obs_mask
4622
4623
4624 allocate(etoptisim(domainmeta%nDomains))
4625 allocate(smoptisim(domainmeta%nDomains))
4626 call eval(parameterset, smoptisim = smoptisim, etoptisim = etoptisim, runoff = runoff)
4627 ngaugestotal = size(runoff, dim = 2)
4628
4629 ! initialize some variables
4630 objective_pd_sm_pd_et_kge_q = 0.0_dp
4631
4632 ! loop over domain - for applying power law later on
4633 do idomain = 1, domainmeta%nDomains
4634
4635 ! check whether the time dimension of SM and ET are same
4636 if (size(etoptisim(idomain)%dataSim, dim = 2) /= size(smoptisim(idomain)%dataSim, dim = 2)) then
4637 call message('***ERROR: mo_objective_funtion: objective_spaef_sm_spaef_et: The time dimension of SM (', &
4638 trim(num2str(size(smoptisim(idomain)%dataSim, dim = 2))),') and ET (', &
4639 trim(num2str(size(etoptisim(idomain)%dataSim, dim = 2))),') are not equal for domain ', &
4640 trim(num2str(idomain)), &
4641 '. Are both of them at same temporal resolution?')
4642 stop
4643 end if
4644
4645 ! get domain information
4646 mask1 = level1(idomain)%mask
4647 ncols1 = level1(idomain)%ncols
4648 nrows1 = level1(idomain)%nrows
4649
4650 ! allocate
4651 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
4652 allocate(pd_time_series_sm(size(smoptisim(idomain)%dataSim, dim = 2)))
4653 allocate(pd_time_series_et(size(etoptisim(idomain)%dataSim, dim = 2)))
4654
4655 allocate(mat1_et(nrows1, ncols1))
4656 allocate(mat2_et(nrows1, ncols1))
4657 allocate(mask_et(nrows1, ncols1))
4658
4659 allocate(mat1_sm(nrows1, ncols1))
4660 allocate(mat2_sm(nrows1, ncols1))
4661 allocate(mask_sm(nrows1, ncols1))
4662
4663 ! initalize
4664 mask_times = .false.
4665 pd_time_series_sm = 0.0_dp
4666 pd_time_series_et = 0.0_dp
4667 objective_kgestreamflow = 0.0_dp
4668 ngaugesdomain = 0._i4
4669
4670 ! calculate pattern similarity criterion
4671 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
4672
4673 mat1_et = unpack(l1_etobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
4674 mat2_et = unpack(etoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
4675 mask_et = unpack(l1_etobs(idomain)%maskObs(:, itime), mask1, .false.)
4676 pd_time_series_et = pd(mat1_et, mat2_et, mask = mask_et, valid = mask_times(itime))
4677
4678 mat1_sm = unpack(l1_smobs(idomain)%dataObs(:, itime), mask1, nodata_dp)
4679 mat2_sm = unpack(smoptisim(idomain)%dataSim(:, itime), mask1, nodata_dp)
4680 mask_sm = unpack(l1_smobs(idomain)%maskObs(:, itime), mask1, .false.)
4681 pd_time_series_sm = pd(mat1_sm, mat2_sm, mask = mask_sm, valid = mask_times(itime))
4682
4683 end do
4684
4685 ! calculate the KGE of streamflow
4686 do gg = 1, ngaugestotal
4687 ! get domain
4688 jdomain = gauge%domainId(gg)
4689 ! only proceed if current domain
4690 if (idomain == jdomain) then
4691 ! count
4692 ngaugesdomain = ngaugesdomain + 1_i4
4693 ! extract runoff
4694 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
4695 ! kge
4696 objective_kgestreamflow = objective_kgestreamflow + &
4697 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
4698 end if
4699 end do
4700 objective_kgestreamflow = objective_kgestreamflow / real(ngaugesdomain, dp)
4701
4702 ! Combine
4703 if (count(mask_times) > 0_i4) then
4704 ! calculate the temporal average of PD, over all domains with power law -domains are weighted equally ( 1 / real(domainMeta%overallNumberOfDomains,dp))**6
4705 objective_pd_sm_pd_et_kge_q = objective_pd_sm_pd_et_kge_q + &
4706 ( &
4707 (1.0_dp - ( &
4708 0.5 * wt_for_optional_data * sum(pd_time_series_sm, mask = mask_times) / real(count(mask_times), dp) + &
4709 0.5 * wt_for_optional_data * sum(pd_time_series_et, mask = mask_times) / real(count(mask_times), dp) + &
4710 1_i4 * objective_kgestreamflow) / (wt_for_optional_data + 1_i4) &
4711 ) / real(domainmeta%overallNumberOfDomains, dp) &
4712 )**6
4713 call message(' pd_sm = ', &
4714 & num2str(sum(pd_time_series_sm, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
4715 &' for domain ', num2str(idomain))
4716 call message(' pd_et = ', &
4717 & num2str(sum(pd_time_series_et, mask = mask_times)/real(count(mask_times), dp), '(F9.5)'), &
4718 &' for domain ', num2str(idomain))
4719 call message(' kge_q = ', num2str(objective_kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
4720 else
4721 call message('***ERROR: mo_objective_funtion: objective_pd_sm_pd_et_kge_q: No ET observations available!')
4722 stop
4723 end if
4724
4725 ! deallocate
4726 deallocate(mask_times)
4727 deallocate(pd_time_series_sm, pd_time_series_et)
4728 deallocate(mat1_sm, mat1_et)
4729 deallocate(mat2_sm, mat2_et)
4730 deallocate(mask_sm, mask_et)
4731 call smoptisim(idomain)%destroy()
4732 call etoptisim(idomain)%destroy()
4733 end do
4734 deallocate(smoptisim)
4735 deallocate(etoptisim)
4736
4737#ifndef MPI
4738 objective_pd_sm_pd_et_kge_q = objective_pd_sm_pd_et_kge_q**onesixth
4739
4740 call message(' objective_pd_sm_pd_et_kge_q = ', num2str(objective_pd_sm_pd_et_kge_q, '(F9.5)'))
4741#endif
4742
4743 END FUNCTION objective_pd_sm_pd_et_kge_q
4744
4745
4746
4747
4748
4749 ! ------------------------------------------------------------------
4750
4751 ! NAME
4752 ! objective_ca_spf
4753
4754 ! PURPOSE
4755 !> \brief Objective function for matching SPF.
4756
4757 !> \details The Classification accuracy (CA) of observed and modeled snow presence flag (SPF)
4758 !> fields are calculated
4759
4760 ! INTENT(IN)
4761 !> \param[in] "real(dp), dimension(:) :: parameterset"
4762 !> \param[in] "procedure(eval_interface) :: eval"
4763
4764 ! RETURN
4765 !> \return real(dp) :: objective_ca_spf &mdash; objective function value
4766 !> (which will be e.g. minimized by an optimization routine like DDS)
4767
4768 ! HISTORY
4769 !> \authors Pallav Shrestha, Ehsan Modiri
4770
4771 !> \date July 2025
4772
4773 ! Modifications:
4774
4775 FUNCTION objective_ca_spf(parameterset, eval)
4776
4777 use mo_optimization_types, only : optidata_sim
4778 use mo_common_variables, only : domainmeta
4779 use mo_global_variables, only : l1_spfobs
4780 use mo_errormeasures, only : ca
4781 use mo_message, only : message
4782 use mo_string_utils, only : num2str
4783
4784 implicit none
4785
4786 real(dp), dimension(:), intent(in) :: parameterset
4787
4788 procedure(eval_interface), INTENT(IN), pointer :: eval
4789
4790 real(dp) :: objective_ca_spf, classification_accuracy
4791
4792 integer(i4) :: idomain
4793
4794 ! simulated snow water equivalent
4795 type(optidata_sim), dimension(:), allocatable :: sweoptisim
4796
4797 ! converted snow presence flag
4798 type(optidata_sim), dimension(:), allocatable :: spfoptisim
4799
4800 ! for sixth root
4801#ifndef MPI
4802 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
4803#endif
4804
4805
4806 allocate(sweoptisim(domainmeta%nDomains))
4807 allocate(spfoptisim(domainmeta%nDomains))
4808 call eval(parameterset, sweoptisim = sweoptisim)
4809
4810 ! initialize some variables
4811 objective_ca_spf = 0.0_dp
4812
4813
4814 ! loop over domain - for applying power law later on
4815 do idomain = 1, domainmeta%nDomains
4816
4817 ! initalize
4818 classification_accuracy = 0.0_dp
4819
4820 ! convert SWE to SPF
4821 call convert_swe_to_spf(sweoptisim(idomain), spfoptisim(idomain))
4822
4823 ! Evaluate classification accuracy
4824 classification_accuracy = ca(l1_spfobs(idomain)%dataObs, spfoptisim(idomain)%dataSim, mask = l1_spfobs(idomain)%maskObs)
4825
4826 objective_ca_spf = objective_ca_spf + &
4827 ( ((1.0_dp - classification_accuracy )) / real(domainmeta%overallNumberOfDomains, dp) )**6
4828
4829
4830 ! deallocate
4831 call sweoptisim(idomain)%destroy()
4832 call spfoptisim(idomain)%destroy()
4833 end do
4834 deallocate(sweoptisim)
4835 deallocate(spfoptisim)
4836
4837#ifndef MPI
4838 objective_ca_spf = objective_ca_spf**onesixth
4839
4840 call message(' objective_ca_spf = ', num2str(objective_ca_spf, '(F9.5)'))
4841 call message(' classification_accuracy_spf = ', num2str(1 - objective_ca_spf, '(F9.5)'))
4842#endif
4843
4844 END FUNCTION objective_ca_spf
4845
4846
4847
4848
4849
4850 ! ------------------------------------------------------------------
4851
4852 ! NAME
4853 ! objective_ca_spf_kge_q
4854
4855 ! PURPOSE
4856 !> \brief Bivaiate objective function for matching SPF and KGE of streamflow.
4857
4858 !> \details The Classification accuracy (CA) of observed and modeled snow presence flag (SPF)
4859 !> fields, and KGE of streamflow simulations at observation locations are calculated
4860
4861 ! INTENT(IN)
4862 !> \param[in] "real(dp), dimension(:) :: parameterset"
4863 !> \param[in] "procedure(eval_interface) :: eval"
4864
4865 ! RETURN
4866 !> \return real(dp) :: objective_ca_spf_kge_q &mdash; objective function value
4867 !> (which will be e.g. minimized by an optimization routine like DDS)
4868
4869 ! HISTORY
4870 !> \authors Pallav Shrestha, Ehsan Modiri
4871
4872 !> \date July 2025
4873
4874 ! Modifications:
4875
4876 FUNCTION objective_ca_spf_kge_q(parameterset, eval)
4877
4878 use mo_optimization_types, only : optidata_sim
4879 use mo_common_variables, only : domainmeta
4880 use mo_global_variables, only : l1_spfobs, wt_for_optional_data
4881 use mo_errormeasures, only : ca, kge
4882 use mo_message, only : message
4883 use mo_string_utils, only : num2str
4884 use mo_mrm_global_variables, only : gauge
4885 use mo_mrm_objective_function_runoff, only : extract_runoff
4886
4887 implicit none
4888
4889 real(dp), dimension(:), intent(in) :: parameterset
4890
4891 procedure(eval_interface), INTENT(IN), pointer :: eval
4892
4893 real(dp) :: objective_ca_spf_kge_q, classification_accuracy, kgestreamflow
4894
4895 integer(i4) :: idomain, jdomain
4896
4897 ! simulated snow water equivalent
4898 type(optidata_sim), dimension(:), allocatable :: sweoptisim
4899
4900 ! converted snow presence flag
4901 type(optidata_sim), dimension(:), allocatable :: spfoptisim
4902
4903 ! simulated streamflow
4904 ! dim1=nTimeSteps, dim2=nGauges
4905 real(dp), allocatable, dimension(:, :) :: runoff
4906
4907 ! gauges counter
4908 integer(i4) :: gg
4909 integer(i4) :: ngaugestotal, ngaugesdomain
4910
4911 ! aggregated simulated
4912 real(dp), dimension(:), allocatable :: runoff_obs
4913
4914 ! measured
4915 real(dp), dimension(:), allocatable :: runoff_agg
4916
4917 ! mask for measured
4918 logical, dimension(:), allocatable :: runoff_obs_mask
4919
4920 ! for sixth root
4921#ifndef MPI
4922 real(dp), parameter :: onesixth = 1.0_dp / 6.0_dp
4923#endif
4924
4925
4926 allocate(sweoptisim(domainmeta%nDomains))
4927 allocate(spfoptisim(domainmeta%nDomains))
4928 call eval(parameterset, sweoptisim = sweoptisim, runoff = runoff)
4929 ngaugestotal = size(runoff, dim = 2)
4930
4931 ! initialize some variables
4932 objective_ca_spf_kge_q = 0.0_dp
4933
4934
4935 ! loop over domain - for applying power law later on
4936 do idomain = 1, domainmeta%nDomains
4937
4938 ! initalize
4939 classification_accuracy = 0.0_dp
4940 kgestreamflow = 0.0_dp
4941 ngaugesdomain = 0_i4
4942
4943 !! SNOW
4944
4945 ! convert SWE to SPF
4946 call convert_swe_to_spf(sweoptisim(idomain), spfoptisim(idomain))
4947
4948 ! Evaluate classification accuracy
4949 classification_accuracy = ca(l1_spfobs(idomain)%dataObs, spfoptisim(idomain)%dataSim, mask = l1_spfobs(idomain)%maskObs)
4950
4951 !! STREAMFLOW
4952
4953 ! calculate the KGE of streamflow
4954 do gg = 1, ngaugestotal
4955 ! get domain
4956 jdomain = gauge%domainId(gg)
4957 ! only proceed if current domain
4958 if (idomain == jdomain) then
4959 ! count
4960 ngaugesdomain = ngaugesdomain + 1_i4
4961 ! extract runoff
4962 call extract_runoff(gg, runoff, runoff_agg, runoff_obs, runoff_obs_mask)
4963 ! kge
4964 kgestreamflow = kgestreamflow + &
4965 kge(runoff_obs, runoff_agg, mask = runoff_obs_mask)
4966 end if
4967 end do
4968 kgestreamflow = kgestreamflow / real(ngaugesdomain, dp)
4969
4970
4971 !! COMBINE
4972 objective_ca_spf_kge_q = objective_ca_spf_kge_q + &
4973 ( &
4974 (1.0_dp - ( &
4975 wt_for_optional_data * classification_accuracy + 1_i4 * kgestreamflow) / (wt_for_optional_data + 1_i4) &
4976 ) / real(domainmeta%overallNumberOfDomains, dp) &
4977 )**6
4978 call message(' ca_spf = ', num2str(classification_accuracy, '(F9.5)'), ' for domain ', num2str(idomain))
4979 call message(' kge_q = ', num2str(kgestreamflow, '(F9.5)'), ' for domain ', num2str(idomain))
4980
4981
4982 ! deallocate
4983 call sweoptisim(idomain)%destroy()
4984 call spfoptisim(idomain)%destroy()
4985 end do
4986 deallocate(sweoptisim)
4987 deallocate(spfoptisim)
4988
4989#ifndef MPI
4990 objective_ca_spf_kge_q = objective_ca_spf_kge_q**onesixth
4991
4992 call message(' objective_ca_spf_kge_q = ', num2str(objective_ca_spf_kge_q, '(F9.5)'))
4993#endif
4994
4995 END FUNCTION objective_ca_spf_kge_q
4996
4997
4998
4999
5000
5001
5002 subroutine create_domain_avg_tws(iDomain, twsOptiSim, tws_catch_avg_domain, &
5003 tws_opti_catch_avg_domain, mask_times)
5004 use mo_optimization_types, only : optidata_sim
5005 use mo_common_constants, only : nodata_dp
5006 use mo_global_variables, only : l1_twsaobs
5007 use mo_moment, only : average
5008 ! current domain Id
5009 integer(i4), intent(in) :: iDomain
5010
5011 ! simulated tws
5012 type(optidata_sim), dimension(:), intent(in) :: twsOptiSim
5013
5014 ! aggregated simulated
5015 real(dp), dimension(:), allocatable, intent(out) :: tws_catch_avg_domain
5016
5017 ! extracted measured
5018 real(dp), dimension(:), allocatable, intent(out) :: tws_opti_catch_avg_domain
5019
5020 ! mask of no data values
5021 logical, dimension(:), allocatable, intent(out) :: mask_times
5022
5023 ! local
5024 ! time loop counter
5025 integer(i4) :: iTime
5026
5027 ! allocate
5028 allocate(mask_times(size(twsoptisim(idomain)%dataSim, dim = 2)))
5029 allocate(tws_catch_avg_domain(size(twsoptisim(idomain)%dataSim, dim = 2)))
5030 allocate(tws_opti_catch_avg_domain(size(twsoptisim(idomain)%dataSim, dim = 2)))
5031
5032 ! initalize
5033 mask_times = .true.
5034 tws_catch_avg_domain = nodata_dp
5035 tws_opti_catch_avg_domain = nodata_dp
5036
5037 ! calculate catchment average evapotranspiration
5038 do itime = 1, size(twsoptisim(idomain)%dataSim, dim = 2)
5039
5040 ! check for enough data points in time for correlation
5041 if (all(.NOT. l1_twsaobs(idomain)%maskObs(:, itime))) then
5042 !write (*,*) 'WARNING: et data at time ', iTime, ' is empty.'
5043 !call message('WARNING: objective_et_kge_catchment_avg: ignored current time step since less than')
5044 !call message(' 10 valid cells available in evapotranspiration observation')
5045 mask_times(itime) = .false.
5046 cycle
5047 end if
5048
5049 tws_catch_avg_domain(itime) = average(l1_twsaobs(idomain)%dataObs(:, itime), &
5050 mask = l1_twsaobs(idomain)%maskObs(:, itime))
5051 tws_opti_catch_avg_domain(itime) = average(twsoptisim(idomain)%dataSim(:, itime), &
5052 mask = l1_twsaobs(idomain)%maskObs(:, itime))
5053 end do
5054
5055 end subroutine create_domain_avg_tws
5056
5057 subroutine create_domain_avg_et(iDomain, etOptiSim, et_catch_avg_domain, &
5058 et_opti_catch_avg_domain, mask_times)
5059 use mo_optimization_types, only : optidata_sim
5060 use mo_common_constants, only : nodata_dp
5061 use mo_global_variables, only : l1_etobs
5062 use mo_moment, only : average
5063 ! current domain Id
5064 integer(i4), intent(in) :: iDomain
5065
5066 ! simulated et
5067 type(optidata_sim), dimension(:), intent(in) :: etOptiSim
5068
5069 ! aggregated simulated
5070 real(dp), dimension(:), allocatable, intent(out) :: et_catch_avg_domain
5071
5072 ! extracted measured
5073 real(dp), dimension(:), allocatable, intent(out) :: et_opti_catch_avg_domain
5074
5075 ! mask of no data values
5076 logical, dimension(:), allocatable, intent(out) :: mask_times
5077
5078 ! local
5079 ! time loop counter
5080 integer(i4) :: iTime
5081
5082 ! allocate
5083 allocate(mask_times(size(etoptisim(idomain)%dataSim, dim = 2)))
5084 allocate(et_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
5085 allocate(et_opti_catch_avg_domain(size(etoptisim(idomain)%dataSim, dim = 2)))
5086
5087 ! initalize
5088 mask_times = .true.
5089 et_catch_avg_domain = nodata_dp
5090 et_opti_catch_avg_domain = nodata_dp
5091
5092 ! calculate catchment average evapotranspiration
5093 do itime = 1, size(etoptisim(idomain)%dataSim, dim = 2)
5094
5095 ! check for enough data points in time for correlation
5096 if (all(.NOT. l1_etobs(idomain)%maskObs(:, itime))) then
5097 !write (*,*) 'WARNING: et data at time ', iTime, ' is empty.'
5098 !call message('WARNING: objective_et_kge_catchment_avg: ignored current time step since less than')
5099 !call message(' 10 valid cells available in evapotranspiration observation')
5100 mask_times(itime) = .false.
5101 cycle
5102 end if
5103
5104 et_catch_avg_domain(itime) = average(l1_etobs(idomain)%dataObs(:, itime), &
5105 mask = l1_etobs(idomain)%maskObs(:, itime))
5106 et_opti_catch_avg_domain(itime) = average(etoptisim(idomain)%dataSim(:, itime), &
5107 mask = l1_etobs(idomain)%maskObs(:, itime))
5108 end do
5109
5110 end subroutine create_domain_avg_et
5111
5112 subroutine convert_tws_to_twsa(twsOptiSim, L1_twsaObs, twsaOptiSim)
5113 use mo_optimization_types, only : optidata_sim, optidata
5114 use mo_moment, only : average
5115 ! simulated tws
5116 type(optidata_sim), intent(in) :: twsOptiSim
5117 ! observed twsa
5118 type(optidata), intent(in) :: L1_twsaObs
5119 ! simulated twsa
5120 type(optidata_sim), intent(inout) :: twsaOptiSim
5121
5122 ! local
5123 integer(i4) :: iCell
5124 real(dp) :: twsa_av_cell
5125
5126 allocate(twsaoptisim%dataSim(size(twsoptisim%dataSim(:, :), dim = 1), size(twsoptisim%dataSim(:, :), dim = 2)))
5127
5128 do icell = 1, size(twsoptisim%dataSim(:, :), dim = 1)
5129 twsa_av_cell = average(twsoptisim%dataSim(icell, :), mask = l1_twsaobs%maskObs(icell, :))
5130 twsaoptisim%dataSim(icell, :) = twsoptisim%dataSim(icell, :) - twsa_av_cell
5131 end do
5132
5133 end subroutine convert_tws_to_twsa
5134
5135 subroutine convert_swe_to_spf(sweOptiSim, spfOptiSim)
5136 use mo_optimization_types, only : optidata_sim, optidata
5137 use mo_global_variables, only : swe_threshold_for_spf
5138 ! simulated swe
5139 type(optidata_sim), intent(in) :: sweOptiSim
5140 ! simulated spf
5141 type(optidata_sim), intent(inout) :: spfOptiSim
5142
5143 ! local
5144 integer(i4) :: iCell
5145 real(dp) :: twsa_av_cell
5146
5147 allocate(spfoptisim%dataSim(size(sweoptisim%dataSim(:, :), dim = 1), size(sweoptisim%dataSim(:, :), dim = 2)))
5148 ! initialize with no snow
5149 spfoptisim%dataSim = 0
5150
5151 where (sweoptisim%dataSim > swe_threshold_for_spf)
5152 ! snow
5153 spfoptisim%dataSim = 1
5154 elsewhere
5155 ! no snow
5156 spfoptisim%dataSim = 0
5157 endwhere
5158
5159 end subroutine convert_swe_to_spf
5160
5161END 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::wt_for_optional_data
integer(i4) wt_for_optional_data
Definition mo_global_variables.f90:125
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_spfobs
type(optidata), dimension(:), allocatable, public l1_spfobs
Definition mo_global_variables.f90:123
mo_global_variables::swe_threshold_for_spf
integer(i4) swe_threshold_for_spf
Definition mo_global_variables.f90:126
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:129
mo_mrm_global_variables
Global variables for mRM only.
Definition mo_mrm_global_variables.f90:25
mo_mrm_global_variables::gauge
type(gaugingstation), public gauge
Definition mo_mrm_global_variables.f90:97
mo_mrm_global_variables::ngaugestotal
integer(i4), public ngaugestotal
Definition mo_mrm_global_variables.f90:84
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:46
mo_objective_function::objective_sm_corr
real(dp) function objective_sm_corr(parameterset, eval)
Objective function for soil moisture.
Definition mo_objective_function.F90:1097
mo_objective_function::convert_swe_to_spf
subroutine convert_swe_to_spf(sweoptisim, spfoptisim)
Definition mo_objective_function.F90:5136
mo_objective_function::objective
real(dp) function, public objective(parameterset, eval, arg1, arg2, arg3)
Wrapper for objective functions.
Definition mo_objective_function.F90:107
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:5059
mo_objective_function::objective_spaef_sm_spaef_et
real(dp) function objective_spaef_sm_spaef_et(parameterset, eval)
Bivariate Objective function for SM and ET.
Definition mo_objective_function.F90:3034
mo_objective_function::objective_spaef_sm
real(dp) function objective_spaef_sm(parameterset, eval)
Objective function for soil moisture.
Definition mo_objective_function.F90:2756
mo_objective_function::objective_pd_sm_pd_et
real(dp) function objective_pd_sm_pd_et(parameterset, eval)
Bivariate Objective function for SM and ET.
Definition mo_objective_function.F90:3976
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:1739
mo_objective_function::objective_sm_pd
real(dp) function objective_sm_pd(parameterset, eval)
Objective function for soil moisture.
Definition mo_objective_function.F90:1230
mo_objective_function::convert_tws_to_twsa
subroutine convert_tws_to_twsa(twsoptisim, l1_twsaobs, twsaoptisim)
Definition mo_objective_function.F90:5113
mo_objective_function::objective_ca_spf
real(dp) function objective_ca_spf(parameterset, eval)
Objective function for matching SPF.
Definition mo_objective_function.F90:4776
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:2473
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:1986
mo_objective_function::objective_ca_spf_kge_q
real(dp) function objective_ca_spf_kge_q(parameterset, eval)
Bivaiate objective function for matching SPF and KGE of streamflow.
Definition mo_objective_function.F90:4877
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:5004
mo_objective_function::objective_pd_et
real(dp) function objective_pd_et(parameterset, eval)
Objective function for ET.
Definition mo_objective_function.F90:3832
mo_objective_function::objective_spaef_et
real(dp) function objective_spaef_et(parameterset, eval)
Objective function for ET.
Definition mo_objective_function.F90:2895
mo_objective_function::objective_pd_sm_kge_q
real(dp) function objective_pd_sm_kge_q(parameterset, eval)
Bivariate Objective function for SM and Q.
Definition mo_objective_function.F90:4160
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:1384
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:1503
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:2341
mo_objective_function::objective_spaef_sm_kge_q
real(dp) function objective_spaef_sm_kge_q(parameterset, eval)
Bivariate Objective function for SM and Q.
Definition mo_objective_function.F90:3213
mo_objective_function::objective_pd_et_kge_q
real(dp) function objective_pd_et_kge_q(parameterset, eval)
Bivariate Objective function for ET and Q.
Definition mo_objective_function.F90:4355
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:1020
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:762
mo_objective_function::objective_spaef_et_kge_q
real(dp) function objective_spaef_et_kge_q(parameterset, eval)
Bivariate Objective function for ET and Q.
Definition mo_objective_function.F90:3404
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:618
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:1885
mo_objective_function::objective_pd_sm_pd_et_kge_q
real(dp) function objective_pd_sm_pd_et_kge_q(parameterset, eval)
Trivariate Objective function for SM, ET and Q.
Definition mo_objective_function.F90:4552
mo_objective_function::objective_spaef_sm_spaef_et_kge_q
real(dp) function objective_spaef_sm_spaef_et_kge_q(parameterset, eval)
Trivariate Objective function for SM, ET and Q.
Definition mo_objective_function.F90:3597
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:2165
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