Line data Source code
1 : !> \file mo_mhm_interface.f90
2 : !> \brief \copybrief mo_mhm_interface
3 : !> \details \copydetails mo_mhm_interface
4 :
5 : !> \brief Module providing interfaces for mHM.
6 : !> \details Interfaces to control the mHM workflow from outside (init, run, get infos, etc.).
7 : !> \authors Sebastian Mueller
8 : !> \version 0.1
9 : !> \date Oct 2021
10 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
11 : !! mHM is released under the LGPLv3+ license \license_note
12 : !> \ingroup f_mhm
13 : module mo_mhm_interface
14 :
15 : use mo_kind, only: i4, dp
16 : use mo_message, only: message, error_message
17 : use mo_string_utils, only: num2str
18 :
19 : #ifdef MPI
20 : use mpi_f08
21 : #endif
22 :
23 : implicit none
24 :
25 : private
26 :
27 : public :: mhm_interface_init
28 : public :: mhm_interface_get_parameter
29 : public :: mhm_interface_get_parameter_number
30 : public :: mhm_interface_run
31 : public :: mhm_interface_run_optimization
32 : public :: mhm_interface_finalize
33 :
34 : contains
35 :
36 : !> \brief initialize mHM from given namelist paths.
37 28 : subroutine mhm_interface_init(namelist_mhm, namelist_mhm_param, namelist_mhm_output, namelist_mrm_output, cwd)
38 : use mo_file, only: &
39 : file_namelist_mhm, &
40 : unamelist_mhm, &
41 : file_namelist_mhm_param, &
42 : unamelist_mhm_param, &
43 : file_defOutput
44 : use mo_mrm_file, only: mrm_file_defOutput => file_defOutput
45 : use mo_common_read_config, only: &
46 : common_read_config
47 : use mo_common_mHM_mRM_read_config, only : &
48 : common_mHM_mRM_read_config, &
49 : check_optimization_settings
50 : use mo_mpr_read_config, only: mpr_read_config
51 : use mo_mhm_read_config, only: mhm_read_config
52 : use mo_read_wrapper, only : read_data
53 : use mo_mrm_init, only: &
54 : mrm_init, &
55 : mrm_configuration
56 : use mo_common_variables, only: &
57 : level0, &
58 : level1, &
59 : itimer, &
60 : domainMeta, &
61 : processMatrix
62 : use mo_common_mHM_mRM_variables, only : &
63 : timeStep, &
64 : simPer, &
65 : optimize, &
66 : opti_function, &
67 : mrm_coupling_mode, &
68 : read_restart
69 : use mo_mhm_messages, only: &
70 : startup_message, &
71 : domain_dir_check_message
72 : use mo_timer, only: &
73 : timers_init, &
74 : timer_start, &
75 : timer_stop, &
76 : timer_get
77 : use mo_startup, only: mhm_initialize
78 : use mo_file, only: &
79 : unamelist_mhm, &
80 : unamelist_mhm_param
81 : use mo_global_variables, only: &
82 : couple_cfg, &
83 : meteo_handler, &
84 : L1_twsaObs, &
85 : L1_etObs, &
86 : L1_neutronsObs, &
87 : L1_smObs, &
88 : BFI_calc
89 : use mo_read_optional_data, only: readOptidataObs
90 : use mo_write_ascii, only: write_configfile
91 : use mo_mhm_bfi, only: calculate_BFI
92 : use mo_os, only: change_dir
93 :
94 : implicit none
95 :
96 : character(*), optional, intent(in) :: namelist_mhm !< path to mHM configuration namelist
97 : character(*), optional, intent(in) :: namelist_mhm_param !< path to mHM parameter namelist
98 : character(*), optional, intent(in) :: namelist_mhm_output !< path to mHM output namelist
99 : character(*), optional, intent(in) :: namelist_mrm_output !< path to mRM output namelist
100 : character(*), optional, intent(in) :: cwd !< desired working directory
101 :
102 : integer(i4) :: domainID, iDomain
103 :
104 : #ifdef MPI
105 : integer :: ierror
106 : integer(i4) :: nproc, rank
107 : #endif
108 :
109 : ! reset nml paths if wanted
110 0 : if (present(namelist_mhm)) file_namelist_mhm = namelist_mhm
111 14 : if (present(namelist_mhm_param)) file_namelist_mhm_param = namelist_mhm_param
112 14 : if (present(namelist_mhm_output)) file_defOutput = namelist_mhm_output
113 14 : if (present(namelist_mrm_output)) mrm_file_defOutput = namelist_mrm_output
114 : ! change working directory
115 14 : if (present(cwd)) call change_dir(cwd)
116 :
117 : ! startup message
118 14 : call startup_message()
119 :
120 : ! coupling configuration
121 14 : call couple_cfg%read_config(file_namelist_mhm, unamelist_mhm)
122 : ! read configs
123 14 : call common_read_config(file_namelist_mhm, unamelist_mhm)
124 14 : call mpr_read_config(file_namelist_mhm, unamelist_mhm, file_namelist_mhm_param, unamelist_mhm_param)
125 14 : call common_mHM_mRM_read_config(file_namelist_mhm, unamelist_mhm)
126 14 : call mhm_read_config(file_namelist_mhm, unamelist_mhm)
127 14 : call couple_cfg%check(domainMeta, optimize)
128 14 : call meteo_handler%config(file_namelist_mhm, unamelist_mhm, optimize, domainMeta, processMatrix, timestep, couple_cfg)
129 14 : mrm_coupling_mode = 2_i4 ! TODO: this shouldn't be needed
130 14 : call mrm_configuration(file_namelist_mhm, unamelist_mhm, file_namelist_mhm_param, unamelist_mhm_param)
131 14 : call check_optimization_settings()
132 :
133 : ! Message about input directories
134 14 : call domain_dir_check_message()
135 :
136 : ! Start timings
137 14 : call timers_init
138 :
139 : ! --------------------------------------------------------------------------
140 : ! READ AND INITIALIZE
141 : ! --------------------------------------------------------------------------
142 14 : itimer = 1
143 : #ifdef MPI
144 : call MPI_Comm_size(domainMeta%comMaster, nproc, ierror)
145 : ! find the number the process is referred to, called rank
146 : call MPI_Comm_rank(domainMeta%comMaster, rank, ierror)
147 : ! ComLocal is a communicator, i.e. a group of processes assigned to the same
148 : ! domain, with a master and subprocesses. Only the master processes of these
149 : ! groups need to read the data. The master process with rank 0 only
150 : ! coordinates the other processes and does not need to read the data.
151 : if (rank > 0 .and. domainMeta%isMasterInComLocal) then
152 : #endif
153 14 : call message()
154 :
155 14 : if (.not. read_restart) then
156 13 : call message(' Read data ...')
157 13 : call timer_start(itimer)
158 : ! for DEM, slope, ... define nGvar local
159 : ! read_data has a domain loop inside
160 13 : call read_data(simPer)
161 13 : call timer_stop(itimer)
162 13 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
163 : end if
164 :
165 : ! read data for every domain
166 14 : itimer = itimer + 1
167 14 : call message(' Initialize domains ...')
168 14 : call timer_start(itimer)
169 14 : call mhm_initialize()
170 14 : call meteo_handler%init_level2(level0, level1)
171 14 : call timer_stop(itimer)
172 14 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
173 14 : if (processMatrix(8, 1) > 0) &
174 13 : call mrm_init(file_namelist_mhm, unamelist_mhm, file_namelist_mhm_param, unamelist_mhm_param)
175 :
176 14 : itimer = itimer + 1
177 14 : call message(' Read forcing and optional data ...')
178 14 : call timer_start(itimer)
179 :
180 40 : do iDomain = 1, domainMeta%nDomains
181 26 : domainID = domainMeta%indices(iDomain)
182 : ! read meteorology now, if it should be loaded in one go
183 26 : if (meteo_handler%single_read(iDomain)) call meteo_handler%prepare_data(1, iDomain, level1, simPer)
184 :
185 : ! read optional optional data if necessary
186 40 : if (optimize) then
187 1 : select case (opti_function)
188 : case(10 : 13, 28)
189 : ! read optional spatio-temporal soil mositure data
190 1 : call readOptidataObs(iDomain, domainID, L1_smObs(iDomain))
191 : case(17)
192 : ! read optional spatio-temporal neutrons data
193 1 : call readOptidataObs(iDomain, domainID, L1_neutronsObs(iDomain))
194 : case(27, 29, 30)
195 : ! read optional spatio-temporal evapotranspiration data
196 0 : call readOptidataObs(iDomain, domainID, L1_etObs(iDomain))
197 : case(15)
198 : ! read optional spatio-temporal tws data
199 1 : call readOptidataObs(iDomain, domainID, L1_twsaObs(iDomain))
200 : case(33)
201 : ! read optional spatio-temporal evapotranspiration data
202 6 : if (domainMeta%optidata(iDomain) == 0 .or. domainMeta%optidata(iDomain) == 5 .or. &
203 : domainMeta%optidata(iDomain) == 6 ) then
204 3 : call readOptidataObs(iDomain, domainID, L1_etObs(iDomain))
205 : end if
206 : ! read optional spatio-temporal tws data
207 6 : if (domainMeta%optidata(iDomain) == 0 .or. domainMeta%optidata(iDomain) == 3 .or. &
208 13 : domainMeta%optidata(iDomain) == 6 ) then
209 3 : call readOptidataObs(iDomain, domainID, L1_twsaObs(iDomain))
210 : end if
211 : end select
212 : end if
213 : end do
214 :
215 : ! calculate observed BFI if wanted
216 14 : if ( optimize .and. opti_function==34 .and. BFI_calc ) call calculate_BFI()
217 :
218 14 : call timer_stop(itimer)
219 14 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
220 :
221 : !this call may be moved to another position as it writes the master config out file for all domains
222 : call write_configfile(meteo_handler%dirPrecipitation, meteo_handler%dirReferenceET, meteo_handler%dirTemperature)
223 :
224 : #ifdef MPI
225 : end if
226 : #endif
227 :
228 28 : end subroutine mhm_interface_init
229 :
230 : !> \brief Get current global parameter value of mHM.
231 0 : subroutine mhm_interface_get_parameter(para)
232 14 : use mo_common_variables, only: global_parameters
233 :
234 : implicit none
235 :
236 : real(dp), dimension(:), allocatable, intent(out) :: para !< global parameter values of mHM
237 :
238 0 : allocate(para(size(global_parameters, dim=1)))
239 :
240 0 : para = global_parameters(:, 3)
241 :
242 0 : end subroutine mhm_interface_get_parameter
243 :
244 : !> \brief Get number of current global parameter value of mHM.
245 0 : subroutine mhm_interface_get_parameter_number(n)
246 0 : use mo_common_variables, only: global_parameters
247 :
248 : implicit none
249 :
250 : integer(i4), intent(out) :: n !< number of global parameter values of mHM
251 :
252 0 : n = size(global_parameters, dim=1)
253 :
254 0 : end subroutine mhm_interface_get_parameter_number
255 :
256 : !> \brief Run mHM with current settings.
257 9 : subroutine mhm_interface_run()
258 : use mo_common_variables, only: &
259 : #ifdef MPI
260 : domainMeta, &
261 : #endif
262 0 : itimer, &
263 : global_parameters
264 : use mo_timer, only: &
265 : timer_start, &
266 : timer_stop, &
267 : timer_get
268 : use mo_mhm_eval, only: mhm_eval
269 :
270 : implicit none
271 :
272 : #ifdef MPI
273 : integer :: ierror
274 : integer(i4) :: nproc, rank
275 :
276 : call MPI_Comm_size(domainMeta%comMaster, nproc, ierror)
277 : ! find the number the process is referred to, called rank
278 : call MPI_Comm_rank(domainMeta%comMaster, rank, ierror)
279 : #endif
280 :
281 9 : itimer = itimer + 1
282 :
283 : ! --------------------------------------------------------------------------
284 : ! call mHM
285 : ! get runoff timeseries if possible (i.e. when domainMeta%doRouting,
286 : ! processMatrix(8,1) > 0)
287 : ! get other model outputs (i.e. gridded fields of model output)
288 : ! --------------------------------------------------------------------------
289 :
290 : #ifdef MPI
291 : if (rank > 0 .and. domainMeta%isMasterInComLocal) then
292 : #endif
293 :
294 9 : call message(' Run mHM')
295 9 : call timer_start(itimer)
296 9 : call mhm_eval(global_parameters(:, 3))
297 9 : call timer_stop(itimer)
298 9 : call message(' in ', trim(num2str(timer_get(itimer), '(F12.3)')), ' seconds.')
299 :
300 : #ifdef MPI
301 : endif
302 : #endif
303 :
304 9 : end subroutine mhm_interface_run
305 :
306 : !> \brief Run mHM optimization with current settings.
307 5 : subroutine mhm_interface_run_optimization()
308 : use mo_common_variables, only: &
309 : #ifdef MPI
310 : domainMeta, &
311 : #endif
312 : itimer, &
313 : dirConfigOut, &
314 : global_parameters,&
315 9 : global_parameters_name, &
316 : processMatrix
317 : use mo_common_mHM_mRM_variables, only : &
318 : opti_function
319 : use mo_timer, only: &
320 : timer_start, &
321 : timer_stop, &
322 : timer_get
323 : use mo_mhm_eval, only: mhm_eval
324 : use mo_objective_function, only: &
325 : #ifdef MPI
326 : objective_subprocess, &
327 : objective_master, &
328 : #endif
329 : objective
330 : use mo_optimization, only: optimization
331 : use mo_mrm_objective_function_runoff, only: &
332 : #ifdef MPI
333 : single_objective_runoff_master, &
334 : single_objective_runoff_subprocess, &
335 : #endif
336 : single_objective_runoff
337 : use mo_write_ascii, only: &
338 : write_optifile, & ! Writing optimized parameter set and objective
339 : write_optinamelist ! Writing optimized parameter set to a namelist
340 :
341 : implicit none
342 :
343 : procedure(mhm_eval), pointer :: eval
344 : procedure(objective), pointer :: obj_func
345 :
346 5 : real(dp) :: funcbest ! best objective function achivied during optimization
347 5 : logical, dimension(:), allocatable :: maskpara ! true = parameter will be optimized, = parameter(i,4) = 1
348 : ! ! false = parameter will not be optimized = parameter(i,4) = 0
349 :
350 : #ifdef MPI
351 : integer :: ierror
352 : integer(i4) :: nproc, rank
353 :
354 : call MPI_Comm_size(domainMeta%comMaster, nproc, ierror)
355 : ! find the number the process is referred to, called rank
356 : call MPI_Comm_rank(domainMeta%comMaster, rank, ierror)
357 : #endif
358 :
359 5 : itimer = itimer + 1
360 5 : call message(' Run mHM optimization')
361 5 : call timer_start(itimer)
362 :
363 5 : eval => mhm_eval
364 :
365 1 : select case(opti_function)
366 : case(1 : 9, 14, 31 : 32)
367 : ! call optimization against only runoff (no other variables)
368 1 : obj_func => single_objective_runoff
369 : #ifdef MPI
370 : if (rank == 0 .and. domainMeta%isMasterInComLocal) then
371 : obj_func => single_objective_runoff_master
372 : call optimization(eval, obj_func, dirConfigOut, funcBest, maskpara)
373 : else if (domainMeta%isMasterInComLocal) then
374 : ! In case of a master process from ComLocal, i.e. a master of a group of
375 : ! processes that are assigned to a single domain, this process calls the
376 : ! objective subroutine directly. The master over all processes collects
377 : ! the data and runs the dds/sce/other opti method.
378 : call single_objective_runoff_subprocess(eval)
379 : end if
380 : #else
381 1 : call optimization(eval, obj_func, dirConfigOut, funcBest, maskpara)
382 : #endif
383 :
384 : case(10 : 13, 15, 17, 27, 28, 29, 30, 33, 34)
385 : ! call optimization for other variables
386 4 : obj_func => objective
387 : #ifdef MPI
388 : if (rank == 0 .and. domainMeta%isMasterInComLocal) then
389 : obj_func => objective_master
390 : call optimization(eval, obj_func, dirConfigOut, funcBest, maskpara)
391 : else if (domainMeta%isMasterInComLocal) then
392 : ! In case of a master process from ComLocal, i.e. a master of a group of
393 : ! processes that are assigned to a single domain, this process calls the
394 : ! objective subroutine directly. The master over all processes collects
395 : ! the data and runs the dds/sce/other opti method.
396 : call objective_subprocess(eval)
397 : end if
398 : #else
399 4 : call optimization(eval, obj_func, dirConfigOut, funcBest, maskpara)
400 : #endif
401 :
402 : case default
403 : call error_message('***ERROR: mhm_driver: The given objective function number ', &
404 10 : trim(adjustl(num2str(opti_function))), ' in mhm.nml is not valid!')
405 : end select
406 :
407 : #ifdef MPI
408 : if (rank == 0 .and. domainMeta%isMasterInComLocal) then
409 : #endif
410 :
411 : ! write a file with final objective function and the best parameter set
412 5 : call write_optifile(funcbest, global_parameters(:, 3), global_parameters_name(:))
413 : ! write a file with final best parameter set in a namlist format
414 5 : call write_optinamelist(processMatrix, global_parameters, maskpara, global_parameters_name(:))
415 5 : deallocate(maskpara)
416 :
417 : #ifdef MPI
418 : end if
419 : #endif
420 :
421 5 : call timer_stop(itimer)
422 5 : call message(' in ', trim(num2str(timer_get(itimer), '(F12.3)')), ' seconds.')
423 :
424 5 : end subroutine mhm_interface_run_optimization
425 :
426 : !> \brief Write mHM restart.
427 14 : subroutine mhm_interface_finalize()
428 : use mo_common_variables, only: &
429 : #ifdef MPI
430 : domainMeta, &
431 : #endif
432 : itimer, &
433 : mhmFileRestartOut, &
434 5 : write_restart, &
435 : processMatrix
436 : use mo_common_mHM_mRM_variables, only : &
437 : optimize
438 : use mo_timer, only: &
439 : timer_start, &
440 : timer_stop, &
441 : timer_get
442 : use mo_restart, only: write_restart_files
443 : use mo_mrm_write, only : mrm_write
444 : use mo_mhm_messages, only: finish_message
445 : use mo_clean_up, only: deallocate_global_variables
446 :
447 : implicit none
448 :
449 : #ifdef MPI
450 : integer :: ierror
451 : integer(i4) :: nproc, rank
452 :
453 : call MPI_Comm_size(domainMeta%comMaster, nproc, ierror)
454 : ! find the number the process is referred to, called rank
455 : call MPI_Comm_rank(domainMeta%comMaster, rank, ierror)
456 : if (rank > 0 .and. domainMeta%isMasterInComLocal) then
457 : #endif
458 :
459 : ! --------------------------------------------------------------------------
460 : ! WRITE RESTART files
461 : ! --------------------------------------------------------------------------
462 23 : if (write_restart .AND. (.NOT. optimize)) then
463 9 : itimer = itimer + 1
464 9 : call message()
465 9 : call message(' Write restart file')
466 9 : call timer_start(itimer)
467 : call write_restart_files(mhmFileRestartOut)
468 9 : call timer_stop(itimer)
469 9 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
470 : end if
471 :
472 : ! --------------------------------------------------------------------------
473 : ! WRITE RUNOFF (INCLUDING RESTART FILES, has to be called after mHM restart
474 : ! files are written)
475 : ! --------------------------------------------------------------------------
476 14 : if (processMatrix(8, 1) > 0) call mrm_write()
477 :
478 : #ifdef MPI
479 : end if
480 : #endif
481 :
482 14 : call finish_message()
483 :
484 : ! clean up all allocated variables
485 14 : call deallocate_global_variables()
486 :
487 14 : end subroutine mhm_interface_finalize
488 :
489 : end module mo_mhm_interface
|
