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