LCOV - code coverage report
Current view: top level - mHM - mo_mhm_interface.F90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 97 108 89.8 %
Date: 2024-04-30 08:53:32 Functions: 4 6 66.7 %

          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

Generated by: LCOV version 1.16