mo_mhm_interface_run.f90

Go to the documentation of this file.

!> \file    mo_mhm_interface_run.f90
!> \brief   \copybrief mo_mhm_interface_run
!> \details \copydetails mo_mhm_interface_run
 
!> \brief   Module providing interfaces for running preconfigured mHM.
!> \details Interfaces to control the mHM run from outside (prepare domain, do timestep, ...).
!> \authors Sebastian Mueller, Matthias Kelbling
!> \version 0.1
!> \date    Jan 2022
!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
!! mHM is released under the LGPLv3+ license \license_note
!> \ingroup f_mhm
module mo_mhm_interface_run
 
  ! forces
  use mo_kind, only: i4, dp
  use mo_message, only: message, error_message
  use mo_string_utils, only : num2str
  ! mhm
  use mo_common_run_variables, only : run_cfg
  use mo_optimization_types, only : optidata_sim
  use mo_common_datetime_type, only : datetimeinfo
  use mo_common_mhm_mrm_variables, only : &
    resolutionrouting, &
    lcyearid, &
    mhmfilerestartin, &
    mrmfilerestartin, &
    ntstepday, &
    optimize, &
    read_restart, &
    simper, &
    timestep, &
    warmingdays, &
    c2tstu, &
    restart_reset_fluxes_states
  use mo_common_variables, only : &
    global_parameters, &
    level0, &
    level1, &
    domainmeta, &
    processmatrix
  use mo_global_variables, only : &
    meteo_handler, &
    l1_throughfall, &
    l1_aetcanopy, &
    l1_aetsealed, &
    l1_aetsoil, &
    l1_baseflow, &
    l1_fastrunoff, &
    l1_infilsoil, &
    l1_inter, &
    l1_melt, &
    l1_neutrons, &
    l1_percol, &
    l1_pet_calc, &
    l1_temp_calc, &
    l1_prec_calc, &
    l1_preeffect, &
    l1_rain, &
    l1_runoffseal, &
    l1_satstw, &
    l1_sealstw, &
    l1_slowrunoff, &
    l1_snow, &
    l1_snowpack, &
    l1_soilmoist, &
    l1_total_runoff, &
    l1_unsatstw, &
    l1_twsaobs, &
    l1_etobs, &
    l1_smobs, &
    l1_neutronsobs, &
    evap_coeff, &
    nsoilhorizons_sm_input, &
    neutron_integral_afast, &
    outputflxstate, &
    timestep_model_outputs, &
    bfi_qbf_sum, &
    bfi_qt_sum
  use mo_init_states, only : variables_default_init, fluxes_states_default_init
  use mo_julian, only : caldat, julday
  use mo_string_utils, only : num2str
  use mo_mhm, only : mhm
  use mo_restart, only : read_restart_states
  use mo_write_fluxes_states, only : mhm_updatedataset, mhm_outputdataset
  use mo_mrm_write_fluxes_states, only : mrm_updatedataset, mrm_outputdataset, gw_outputdataset, gw_updatedataset
  use mo_constants, only : hoursecs
  use mo_common_variables, only : resolutionhydrology
  use mo_mrm_global_variables, only : &
    inflowgauge, &
    l11_c1, &
    l11_c2, &
    l11_l1_id, &
    l11_tsrout, &
    l11_fromn, &
    l11_length, &
    sink_cells, &
    l11_nlinkfracfpimp, &
    l11_noutlets, &
    l11_netperm, &
    l11_qmod, &
    l11_qout, &
    l11_qtin, &
    l11_qtr, &
    l11_slope, &
    l11_ton, &
    l1_l11_id, &
    domain_mrm, &
    level11, &
    mrm_runoff, &
    outputflxstate_mrm, &
    timestep_model_outputs_mrm, &
    gw_coupling, &
    l0_river_head_mon_sum, &
    riv_temp_pcs
  use mo_mpr_global_variables, only : &
    l1_harsamcoeff, &
    l1_prietayalpha, &
    l1_aeroresist, &
    l1_alpha, &
    l1_degday, &
    l1_degdayinc, &
    l1_degdaymax, &
    l1_degdaynopre, &
    l1_fasp, &
    l1_froots, &
    l1_fsealed, &
    l1_jarvis_thresh_c1, &
    l1_kbaseflow, &
    l1_kperco, &
    l1_kslowflow, &
    l1_karstloss, &
    l1_kfastflow, &
    l1_maxinter, &
    l1_petlaicorfactor, &
    l1_sealedthresh, &
    l1_soilmoistexp, &
    l1_soilmoistfc, &
    l1_soilmoistsat, &
    l1_surfresist, &
    l1_tempthresh, &
    l1_unsatthresh, &
    l1_wiltingpoint, &
    l1_no_count, &
    l1_bulkdens, &
    l1_latticewater, &
    l1_cosmicl3, &
    horizondepth_mhm, &
    nsoilhorizons_mhm
  use mo_mrm_init, only : variables_default_init_routing, fluxes_states_default_init_routing
  use mo_mrm_mpr, only : mrm_update_param
  use mo_mrm_restart, only : mrm_read_restart_states
  use mo_mrm_routing, only : mrm_routing
  use mo_utils, only : ge
  use mo_mrm_river_head, only: calc_river_head
  use mo_mpr_eval, only : mpr_eval
 
contains
 
  !> \brief prepare single run of mHM
  subroutine mhm_interface_run_prepare(parameterset, opti_domain_indices, runoff_present, BFI_present)
    implicit none
    !> a set of global parameter (gamma) to run mHM, DIMENSION [no. of global_Parameters]
    real(dp), dimension(:), optional, intent(in) :: parameterset
    !> selected domains for optimization
    integer(i4), dimension(:), optional, intent(in) :: opti_domain_indices
    !> whether runoff is present
    logical, optional, intent(in) :: runoff_present
    !> whether BFI is present
    logical, optional, intent(in) :: BFI_present
 
    integer(i4) :: i, iDomain, domainID
 
    ! run_cfg%output_runoff = .false. by default
    if (present(runoff_present)) run_cfg%output_runoff = runoff_present
    ! run_cfg%output_BFI = .false. by default
    if (present(bfi_present)) run_cfg%output_BFI = bfi_present
 
    ! store current parameter set
    allocate(run_cfg%parameterset(size(global_parameters, dim=1)))
    if (.not. present(parameterset) .and. optimize) then
      call error_message("mhm_interface_run_prepare: Can't optimize without parameter!")
    else if (.not. present(parameterset)) then
      run_cfg%parameterset = global_parameters(:, 3)
    else
      run_cfg%parameterset = parameterset
    end if
 
    ! store domain indices for domain loop
    if (optimize .and. present(opti_domain_indices)) then
      run_cfg%nDomains = size(opti_domain_indices)
      allocate(run_cfg%domain_indices(run_cfg%nDomains))
      run_cfg%domain_indices = opti_domain_indices
    else
      run_cfg%nDomains = domainmeta%nDomains
      allocate(run_cfg%domain_indices(run_cfg%nDomains))
      run_cfg%domain_indices = [(i, i=1, run_cfg%nDomains)]
    end if
 
    !----------------------------------------------------------
    ! Check optionals and initialize
    !----------------------------------------------------------
    if (run_cfg%output_runoff) then
      do i = 1, run_cfg%nDomains
        idomain = run_cfg%get_domain_index(i)
        domainid = domainmeta%indices(idomain)
        if (.not. domainmeta%doRouting(idomain)) then
          call error_message("***ERROR: runoff for domain", trim(num2str(domainid)),&
                        "can not be produced, since routing process is off in Process Matrix")
        end if
      end do
    end if
 
    ! prepare BFI calculation
    if (run_cfg%output_BFI) then
      allocate(bfi_qbf_sum(run_cfg%nDomains))
      allocate(bfi_qt_sum(run_cfg%nDomains))
      bfi_qbf_sum = 0.0_dp
      bfi_qt_sum = 0.0_dp
    end if
 
    if (read_restart) then
      do i = 1, run_cfg%nDomains
        idomain = run_cfg%get_domain_index(i)
        domainid = domainmeta%indices(idomain)
        ! this reads the eff. parameters and optionally the states and fluxes
        call read_restart_states(idomain, domainid, mhmfilerestartin(idomain))
      end do
      if (restart_reset_fluxes_states) then
        call message('    Resetting mHM states and fluxes from restart files ...')
        call fluxes_states_default_init()
      end if
    else
      call variables_default_init()
      call mpr_eval(run_cfg%parameterset)
      if (processmatrix(8, 1) > 0) call variables_default_init_routing()
    end if
 
    allocate(run_cfg%L1_fNotSealed(size(l1_fsealed, 1), size(l1_fsealed, 2), size(l1_fsealed, 3)))
    run_cfg%L1_fNotSealed = 1.0_dp - l1_fsealed
 
  end subroutine mhm_interface_run_prepare
 
  !> \brief get number of domains for looping
  subroutine mhm_interface_run_get_ndomains(ndomains)
    implicit none
    integer(i4), intent(inout) :: ndomains !< number of domains
    ndomains = run_cfg%nDomains
  end subroutine mhm_interface_run_get_ndomains
 
  !> \brief prepare single domain to run mHM on
  subroutine mhm_interface_run_prepare_domain(domain, etOptiSim, twsOptiSim, neutronsOptiSim, smOptiSim)
    implicit none
    !> domain loop counter
    integer(i4), intent(in), optional :: domain
    !> returns soil moisture time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
    type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim
    !> dim1=ncells, dim2=time
    type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim
    !> returns evapotranspiration time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
    type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim
    !> returns tws time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
    type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim
 
    integer(i4) :: iDomain, domainID
 
    ! get domain index
    run_cfg%selected_domain = 1_i4
    if (present(domain)) run_cfg%selected_domain = domain
    idomain = run_cfg%get_domain_index(run_cfg%selected_domain)
    domainid = domainmeta%indices(idomain)
 
    ! evapotranspiration optimization
    if (present(etoptisim)) call etoptisim(idomain)%init(l1_etobs(idomain))
    ! total water storage optimization
    if (present(twsoptisim)) call twsoptisim(idomain)%init(l1_twsaobs(idomain))
    ! neutrons optimization
    if (present(neutronsoptisim)) call neutronsoptisim(idomain)%init(l1_neutronsobs(idomain))
    ! sm optimization
    if (present(smoptisim)) call smoptisim(idomain)%init(l1_smobs(idomain))
 
    ! get Domain information
    run_cfg%nCells = level1(idomain)%nCells
    run_cfg%mask1 => level1(idomain)%mask
    run_cfg%s1 = level1(idomain)%iStart
    run_cfg%e1 = level1(idomain)%iEnd
 
    if (domainmeta%doRouting(idomain)) then
      ! ----------------------------------------
      ! initialize factor between routing resolution and hydrologic model resolution
      ! ----------------------------------------
      run_cfg%tsRoutFactor = 1_i4
      allocate(run_cfg%InflowDischarge(size(inflowgauge%Q, dim = 2)))
      run_cfg%InflowDischarge = 0._dp
 
      ! read states from restart
      if (read_restart) then
        call mrm_read_restart_states(idomain, domainid, mrmfilerestartin(idomain))
        if (restart_reset_fluxes_states) then
          call message('    Resetting mRM states and fluxes from restart files for domain ', num2str(idomain), ' ...')
          call fluxes_states_default_init_routing(idomain)
        end if
      end if
      ! get Domain information at L11 if routing is activated
      run_cfg%s11 = level11(idomain)%iStart
      run_cfg%e11 = level11(idomain)%iEnd
      run_cfg%mask11 => level11(idomain)%mask
 
      ! initialize routing parameters (has to be called for routing option 2)
      if ((processmatrix(8, 1) .eq. 2) .or. (processmatrix(8, 1) .eq. 3)) &
          call mrm_update_param( &
            idomain, run_cfg%parameterset(processmatrix(8, 3) - processmatrix(8, 2) + 1 : processmatrix(8, 3)))
      ! initialize variable for runoff for routing
      allocate(run_cfg%RunToRout(run_cfg%e1 - run_cfg%s1 + 1))
      run_cfg%RunToRout = 0._dp
 
      if ( riv_temp_pcs%active ) then
        ! set indices for current L11 domain
        riv_temp_pcs%s11 = run_cfg%s11
        riv_temp_pcs%e11 = run_cfg%e11
        ! allocate current L1 lateral components
        call riv_temp_pcs%alloc_lateral(run_cfg%nCells)
      end if
    end if
 
    ! init datetime variable
    call run_cfg%domainDateTime%init(idomain)
 
    ! init time-loop
    run_cfg%time_step = 0_i4
 
  end subroutine mhm_interface_run_prepare_domain
 
  !> \brief check if current time loop is finished
  subroutine mhm_interface_run_finished(time_loop_finished)
    implicit none
    logical, intent(inout) :: time_loop_finished !< flag to indicate end of timeloop
    time_loop_finished = run_cfg%time_step == run_cfg%domainDateTime%nTimeSteps
  end subroutine mhm_interface_run_finished
 
  !> \brief do one time-step on current domain
  subroutine mhm_interface_run_do_time_step()
    implicit none
 
    integer(i4) :: iDomain, domainID, tt, jj, s1, e1
 
    ! increment time step count (first input is 0)
    run_cfg%time_step = run_cfg%time_step + 1_i4
    ! current time counter and domain indices
    tt = run_cfg%time_step
    s1 = run_cfg%s1
    e1 = run_cfg%e1
    idomain = run_cfg%get_domain_index(run_cfg%selected_domain)
    domainid = domainmeta%indices(idomain)
 
    call run_cfg%domainDateTime%update_LAI_timestep()
 
    ! update the meteo-handler
    call meteo_handler%update_timestep( &
      tt=tt, &
      time=run_cfg%domainDateTime%newTime - 0.5_dp, &
      idomain=idomain, &
      level1=level1, &
      simper=simper)
 
    ! get corrected pet
    call meteo_handler%get_corrected_pet(pet_calc=l1_pet_calc(s1 : e1), &
      ! pet calculation dependencies
      petlaicorfactorl1 = l1_petlaicorfactor(s1 : e1, run_cfg%domainDateTime%iLAI, run_cfg%domainDateTime%yId), &
      fasp              = l1_fasp(s1 : e1, 1, 1), &
      harsamcoeff       = l1_harsamcoeff(s1 : e1, 1, 1), &
      latitude          = pack(level1(idomain)%y, level1(idomain)%mask), &
      prietayalpha      = l1_prietayalpha(s1 : e1, run_cfg%domainDateTime%iLAI, 1), &
      aeroresist        = l1_aeroresist(s1 : e1, run_cfg%domainDateTime%iLAI, run_cfg%domainDateTime%yId), &
      surfresist        = l1_surfresist(s1 : e1, run_cfg%domainDateTime%iLAI, 1))
 
    ! get temperature and precipitation
    call meteo_handler%get_temp(temp_calc=l1_temp_calc(s1 : e1))
    call meteo_handler%get_prec(prec_calc=l1_prec_calc(s1 : e1))
 
    ! -------------------------------------------------------------------------
    ! ARGUMENT LIST KEY FOR mHM
    ! -------------------------------------------------------------------------
    !  C    CONFIGURATION
    !  F    FORCING DATA L2
    !  Q    INFLOW FROM UPSTREAM AREAS
    !  L0   MORPHOLOGIC DATA L0
    !  L1   MORPHOLOGIC DATA L1
    !  L11  MORPHOLOGIC DATA L11
    !  P    GLOBAL PARAMETERS
    !  E    EFFECTIVE PARAMETER FIELDS (L1, L11 levels)
    !  S    STATE VARIABLES L1
    !  X    FLUXES (L1, L11 levels)
    ! --------------------------------------------------------------------------
    call mhm( &
      read_states            = read_restart, & ! IN C
      tt                     = tt, &
      time                   = run_cfg%domainDateTime%newTime - 0.5_dp, &
      processmatrix          = processmatrix, &
      horizon_depth          = horizondepth_mhm, & ! IN C
      ncells1                = run_cfg%nCells, &
      nhorizons_mhm          = nsoilhorizons_mhm, &
      c2tstu                 = c2tstu,  & ! IN C
      neutron_integral_afast = neutron_integral_afast, & ! IN C
      evap_coeff             = evap_coeff, &
      fsealed1               = l1_fsealed(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT L1
      interc                 = l1_inter(s1 : e1), &
      snowpack               = l1_snowpack(s1 : e1), &
      sealedstorage          = l1_sealstw(s1 : e1), & ! INOUT S
      soilmoisture           = l1_soilmoist(s1 : e1, :), &
      unsatstorage           = l1_unsatstw(s1 : e1), &
      satstorage             = l1_satstw(s1 : e1), & ! INOUT S
      neutrons               = l1_neutrons(s1 : e1), & ! INOUT S
      pet_calc               = l1_pet_calc(s1 : e1), & ! INOUT X
      temp_calc              = l1_temp_calc(s1 : e1), & ! INOUT X
      prec_calc              = l1_prec_calc(s1 : e1), & ! INOUT X
      aet_soil               = l1_aetsoil(s1 : e1, :), &
      aet_canopy             = l1_aetcanopy(s1 : e1), &
      aet_sealed             = l1_aetsealed(s1 : e1), & ! INOUT X
      baseflow               = l1_baseflow(s1 : e1), &
      infiltration           = l1_infilsoil(s1 : e1, :), &
      fast_interflow         = l1_fastrunoff(s1 : e1), & ! INOUT X
      melt                   = l1_melt(s1 : e1), &
      perc                   = l1_percol(s1 : e1), &
      prec_effect            = l1_preeffect(s1 : e1), &
      rain                   = l1_rain(s1 : e1), & ! INOUT X
      runoff_sealed          = l1_runoffseal(s1 : e1), &
      slow_interflow         = l1_slowrunoff(s1 : e1), &
      snow                   = l1_snow(s1 : e1), & ! INOUT X
      throughfall            = l1_throughfall(s1 : e1), &
      total_runoff           = l1_total_runoff(s1 : e1), & ! INOUT X
      ! MPR
      alpha                  = l1_alpha(s1 : e1, 1, run_cfg%domainDateTime%yId), &
      deg_day_incr           = l1_degdayinc(s1 : e1, 1, run_cfg%domainDateTime%yId), &
      deg_day_max            = l1_degdaymax(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
      deg_day_noprec         = l1_degdaynopre(s1 : e1, 1, run_cfg%domainDateTime%yId), &
      deg_day                = l1_degday(s1 : e1, 1, 1), & ! INOUT E1
      frac_roots             = l1_froots(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
      interc_max             = l1_maxinter(s1 : e1, run_cfg%domainDateTime%iLAI, 1), &
      karst_loss             = l1_karstloss(s1 : e1, 1, 1), & ! INOUT E1
      k0                     = l1_kfastflow(s1 : e1, 1, run_cfg%domainDateTime%yId), &
      k1                     = l1_kslowflow(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
      k2                     = l1_kbaseflow(s1 : e1, 1, run_cfg%domainDateTime%yId), &
      kp                     = l1_kperco(s1 : e1, 1, run_cfg%domainDateTime%yId), & ! INOUT E1
      soil_moist_fc          = l1_soilmoistfc(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
      soil_moist_sat         = l1_soilmoistsat(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
      soil_moist_exponen     = l1_soilmoistexp(s1 : e1, :, run_cfg%domainDateTime%yId), &
      jarvis_thresh_c1       = l1_jarvis_thresh_c1(s1 : e1, 1, 1), & ! INOUT E1
      temp_thresh            = l1_tempthresh(s1 : e1, 1, run_cfg%domainDateTime%yId), &
      unsat_thresh           = l1_unsatthresh(s1 : e1, 1, 1), & ! INOUT E1
      water_thresh_sealed    = l1_sealedthresh(s1 : e1, 1, 1), & ! INOUT E1
      wilting_point          = l1_wiltingpoint(s1 : e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
      ! >> neutron count
      no_count               = l1_no_count(s1:e1, 1, 1),  &                     ! INOUT E1
      bulkdens               = l1_bulkdens(s1:e1,     :, run_cfg%domainDateTime%yId), & ! INOUT E1
      latticewater           = l1_latticewater(s1:e1, :, run_cfg%domainDateTime%yId), & ! INOUT E1
      cosmicl3               = l1_cosmicl3(s1:e1,     :, run_cfg%domainDateTime%yId)  & ! INOUT E1
    )
 
    ! call mRM routing
    run_cfg%doRoute = .false.
    if (domainmeta%doRouting(idomain)) then
      ! set discharge timestep
      run_cfg%iDischargeTS = ceiling(real(tt, dp) / real(ntstepday, dp))
      ! set input variables for routing
      if (processmatrix(8, 1) .eq. 1) then
        ! >>>
        ! >>> original Muskingum routing, executed every time
        ! >>>
        run_cfg%doRoute = .true.
        run_cfg%tsRoutFactorIn = 1._dp
        run_cfg%timestep_rout = timestep
        run_cfg%RunToRout = l1_total_runoff(s1 : e1) ! runoff [mm TS-1] mm per timestep
        run_cfg%InflowDischarge = inflowgauge%Q(run_cfg%iDischargeTS, :) ! inflow discharge in [m3 s-1]
        !
      else if ((processmatrix(8, 1) .eq. 2) .or. &
               (processmatrix(8, 1) .eq. 3)) then
        ! >>>
        ! >>> adaptive timestep
        ! >>>
        run_cfg%doRoute = .false.
        ! calculate factor
        run_cfg%tsRoutFactor = l11_tsrout(idomain) / (timestep * hoursecs)
        ! print *, 'routing factor: ', tsRoutFactor
        ! prepare routing call
        if (run_cfg%tsRoutFactor .lt. 1._dp) then
          ! ----------------------------------------------------------------
          ! routing timesteps are shorter than hydrologic time steps
          ! ----------------------------------------------------------------
          ! set all input variables
          run_cfg%tsRoutFactorIn = run_cfg%tsRoutFactor
          run_cfg%RunToRout = l1_total_runoff(s1 : e1) ! runoff [mm TS-1] mm per timestep
          run_cfg%InflowDischarge = inflowgauge%Q(run_cfg%iDischargeTS, :) ! inflow discharge in [m3 s-1]
          run_cfg%timestep_rout = timestep
          run_cfg%doRoute = .true.
        else
          ! ----------------------------------------------------------------
          ! routing timesteps are longer than hydrologic time steps
          ! ----------------------------------------------------------------
          ! set all input variables
          run_cfg%tsRoutFactorIn = run_cfg%tsRoutFactor
          ! Runoff is accumulated in [mm]
          run_cfg%RunToRout = run_cfg%RunToRout + l1_total_runoff(s1 : e1)
          run_cfg%InflowDischarge = run_cfg%InflowDischarge + inflowgauge%Q(run_cfg%iDischargeTS, :)
          ! reset tsRoutFactorIn if last period did not cover full period
          if ((tt == run_cfg%domainDateTime%nTimeSteps) .and. (mod(tt, nint(run_cfg%tsRoutFactorIn)) /= 0_i4)) &
            run_cfg%tsRoutFactorIn = mod(tt, nint(run_cfg%tsRoutFactorIn))
          if ((mod(tt, nint(run_cfg%tsRoutFactorIn)) .eq. 0_i4) .or. (tt .eq. run_cfg%domainDateTime%nTimeSteps)) then
            ! Inflow discharge is given as flow-rate and has to be converted to [m3]
            run_cfg%InflowDischarge = run_cfg%InflowDischarge / run_cfg%tsRoutFactorIn
            run_cfg%timestep_rout = timestep * nint(run_cfg%tsRoutFactorIn, i4)
            run_cfg%doRoute = .true.
          end if
        end if
      end if
      ! prepare temperature routing
      if ( riv_temp_pcs%active ) then
        ! init riv-temp from current air temp
        if ( tt .eq. 1_i4 ) call riv_temp_pcs%init_riv_temp( &
          temp_air     = l1_temp_calc(s1 : e1), &
          efecarea     = level1(idomain)%CellArea * 1.e-6_dp, &
          l1_l11_id    = l1_l11_id(s1 : e1), &
          l11_areacell = level11(idomain)%CellArea * 1.e-6_dp, &
          l11_l1_id    = l11_l1_id(run_cfg%s11 : run_cfg%e11), &
          map_flag     = ge(resolutionrouting(idomain), resolutionhydrology(idomain)) &
        )
        ! get riv-temp specific meteo arrays
        call meteo_handler%get_ssrd(riv_temp_pcs%L1_ssrd_calc)
        call meteo_handler%get_strd(riv_temp_pcs%L1_strd_calc)
        call meteo_handler%get_tann(riv_temp_pcs%L1_tann_calc)
        ! accumulate source Energy at L1 level
        call riv_temp_pcs%acc_source_E( &
          fsealed_area_fraction = l1_fsealed(s1 : e1, 1, run_cfg%domainDateTime%yId), &
          fast_interflow        = l1_fastrunoff(s1 : e1), &
          slow_interflow        = l1_slowrunoff(s1 : e1), &
          baseflow              = l1_baseflow(s1 : e1), &
          direct_runoff         = l1_runoffseal(s1 : e1), &
          temp_air              = l1_temp_calc(s1 : e1) &
        )
        ! if routing should be performed, scale source energy to L11 level
        if ( run_cfg%doRoute ) call riv_temp_pcs%finalize_source_E( &
          efecarea     = level1(idomain)%CellArea * 1.e-6_dp, &
          l1_l11_id    = l1_l11_id(s1 : e1), &
          l11_areacell = level11(idomain)%CellArea * 1.e-6_dp, &
          l11_l1_id    = l11_l1_id(run_cfg%s11 : run_cfg%e11), &
          timestep     = run_cfg%timestep_rout, &
          map_flag     = ge(resolutionrouting(idomain), resolutionhydrology(idomain)) &
        )
      end if
      ! -------------------------------------------------------------------
      ! execute routing
      ! -------------------------------------------------------------------
      if (run_cfg%doRoute) call mrm_routing(&
        ! general INPUT variables
        read_restart, &
        processmatrix(8, 1), & ! parse process Case to be used
        run_cfg%parameterset(processmatrix(8, 3) - processmatrix(8, 2) + 1 : processmatrix(8, 3)), & ! routing par.
        run_cfg%RunToRout, & ! runoff [mm TS-1] mm per timestep old: L1_total_runoff_in(run_cfg%s1:run_cfg%e1, tt), &
        level1(idomain)%CellArea * 1.e-6_dp, &
        l1_l11_id(s1 : e1), &
        level11(idomain)%CellArea * 1.e-6_dp, &
        l11_l1_id(run_cfg%s11 : run_cfg%e11), &
        l11_netperm(run_cfg%s11 : run_cfg%e11), & ! routing order at L11
        l11_fromn(run_cfg%s11 : run_cfg%e11), & ! link source at L11
        l11_ton(run_cfg%s11 : run_cfg%e11), & ! link target at L11
        l11_noutlets(idomain), & ! number of outlets
        run_cfg%timestep_rout, & ! timestep of runoff to rout [h]
        run_cfg%tsRoutFactorIn, & ! simulate timestep in [h]
        level11(idomain)%nCells, & ! number of Nodes
        domain_mrm(idomain)%nInflowGauges, &
        domain_mrm(idomain)%InflowGaugeIndexList(:), &
        domain_mrm(idomain)%InflowGaugeHeadwater(:), &
        domain_mrm(idomain)%InflowGaugeNodeList(:), &
        run_cfg%InflowDischarge, &
        domain_mrm(idomain)%nGauges, &
        domain_mrm(idomain)%gaugeIndexList(:), &
        domain_mrm(idomain)%gaugeNodeList(:), &
        ge(resolutionrouting(idomain), resolutionhydrology(idomain)), &
        ! original routing specific input variables
        l11_length(run_cfg%s11 : run_cfg%e11 - 1), & ! link length
        l11_slope(run_cfg%s11 : run_cfg%e11 - 1), &
        l11_nlinkfracfpimp(run_cfg%s11 : run_cfg%e11, run_cfg%domainDateTime%yId), & ! fraction of impervious layer at L11 scale
        sink_cells(idomain)%ids, &
        ! general INPUT/OUTPUT variables
        l11_c1(run_cfg%s11 : run_cfg%e11), & ! first muskingum parameter
        l11_c2(run_cfg%s11 : run_cfg%e11), & ! second muskigum parameter
        l11_qout(run_cfg%s11 : run_cfg%e11), & ! routed runoff flowing out of L11 cell
        l11_qtin(run_cfg%s11 : run_cfg%e11, :), & ! inflow water into the reach at L11
        l11_qtr(run_cfg%s11 : run_cfg%e11, :), & !
        l11_qmod(run_cfg%s11 : run_cfg%e11), &
        mrm_runoff(tt, :) &
      )
      ! -------------------------------------------------------------------
      ! groundwater coupling
      ! -------------------------------------------------------------------
      if (gw_coupling) then
        call calc_river_head(idomain, l11_qmod, l0_river_head_mon_sum)
      end if
      ! -------------------------------------------------------------------
      ! reset variables
      ! -------------------------------------------------------------------
      if (processmatrix(8, 1) .eq. 1) then
        ! reset Input variables
        run_cfg%InflowDischarge = 0._dp
        run_cfg%RunToRout = 0._dp
      else if ((processmatrix(8, 1) .eq. 2) .or. (processmatrix(8, 1) .eq. 3)) then
        if ((.not. (run_cfg%tsRoutFactorIn .lt. 1._dp)) .and. run_cfg%doRoute) then
          do jj = 1, nint(run_cfg%tsRoutFactorIn) ! BUG: this should start at 2
            mrm_runoff(tt - jj + 1, :) = mrm_runoff(tt, :)
          end do
          ! reset Input variables
          run_cfg%InflowDischarge = 0._dp
          run_cfg%RunToRout = 0._dp
          ! reset lateral fluxes and time-step counter if routing was done
          if ( riv_temp_pcs%active ) call riv_temp_pcs%reset_timestep()
        end if
        ! if routing is done every time-step, reset river-temp time step
        if (run_cfg%tsRoutFactor .lt. 1._dp .and. riv_temp_pcs%active ) call riv_temp_pcs%reset_timestep()
      end if
    end if
 
    ! output only for evaluation period
    run_cfg%domainDateTime%tIndex_out = (tt - warmingdays(idomain) * ntstepday) ! tt if write out of warming period
 
    call run_cfg%domainDateTime%increment()
 
    ! update the year-dependent domainDateTime%yId (land cover id)
    if (run_cfg%domainDateTime%is_new_year .and. tt < run_cfg%domainDateTime%nTimeSteps) then
      run_cfg%domainDateTime%yId = lcyearid(run_cfg%domainDateTime%year, idomain)
    end if
 
    ! calculate BFI releated after warming days if wanted
    if ( run_cfg%output_BFI .and. (run_cfg%domainDateTime%tIndex_out > 0_i4) ) then
      bfi_qbf_sum(idomain) = bfi_qbf_sum(idomain) &
        + sum(l1_baseflow(s1 : e1) * level1(idomain)%CellArea) / level1(idomain)%nCells
      bfi_qt_sum(idomain) = bfi_qt_sum(idomain) &
        + sum(l1_total_runoff(s1 : e1) * level1(idomain)%CellArea) / level1(idomain)%nCells
    end if
 
  end subroutine mhm_interface_run_do_time_step
 
  !> \brief write output after current time-step
  subroutine mhm_interface_run_write_output()
    implicit none
    integer(i4) :: iDomain, tt, s0, e0
 
    ! get time step
    tt = run_cfg%time_step
 
    ! get domain index
    idomain = run_cfg%get_domain_index(run_cfg%selected_domain)
 
    if ( (.not. optimize) .and. (run_cfg%domainDateTime%tIndex_out > 0_i4)) then
      if (any(outputflxstate_mrm) .AND. (domainmeta%doRouting(idomain))) then
 
        if (run_cfg%domainDateTime%tIndex_out == 1) then
          run_cfg%nc_mrm = mrm_outputdataset(idomain, run_cfg%mask11)
        end if
 
        ! update Dataset (riv-temp as optional input)
        if ( riv_temp_pcs%active ) then
          call mrm_updatedataset(run_cfg%nc_mrm, &
            l11_qmod(run_cfg%s11 : run_cfg%e11), riv_temp_pcs%river_temp(riv_temp_pcs%s11 : riv_temp_pcs%e11))
        else
          call mrm_updatedataset(run_cfg%nc_mrm, l11_qmod(run_cfg%s11 : run_cfg%e11))
        end if
 
        ! write data
        if (run_cfg%domainDateTime%writeout(timestep_model_outputs_mrm, tt)) then
          call run_cfg%nc_mrm%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
        end if
 
        if(tt == run_cfg%domainDateTime%nTimeSteps) then
          call run_cfg%nc_mrm%close()
        end if
 
        if ( gw_coupling ) then
          ! create
          if (run_cfg%domainDateTime%tIndex_out == 1) then
            run_cfg%nc_gw = gw_outputdataset(idomain, level0(idomain)%mask)
          end if
          ! add data
          s0 = level0(idomain)%iStart
          e0 = level0(idomain)%iEnd
          call gw_updatedataset(run_cfg%nc_gw, l0_river_head_mon_sum(s0 : e0))
          ! write
          if (run_cfg%domainDateTime%writeout(-2, tt)) then ! -2 for monthly
            call run_cfg%nc_gw%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
          end if
          ! close
          if(tt == run_cfg%domainDateTime%nTimeSteps) then
            call run_cfg%nc_gw%close()
          end if
        end if
 
      end if
 
      if (any(outputflxstate)) then
 
        if (run_cfg%domainDateTime%tIndex_out == 1) then
          run_cfg%nc_mhm = mhm_outputdataset(idomain, run_cfg%mask1)
        end if
 
        call mhm_updatedataset(&
          run_cfg%nc_mhm, &
          run_cfg%s1, run_cfg%e1, &
          l1_fsealed(:, 1, run_cfg%domainDateTime%yId), &
          run_cfg%L1_fNotSealed(:, 1, run_cfg%domainDateTime%yId), &
          l1_inter, &
          l1_snowpack, &
          l1_soilmoist, &
          l1_soilmoistsat(:, :, run_cfg%domainDateTime%yId), &
          l1_sealstw, &
          l1_unsatstw, &
          l1_satstw, &
          l1_neutrons, &
          l1_pet_calc, &
          l1_aetsoil, &
          l1_aetcanopy, &
          l1_aetsealed, &
          l1_total_runoff, &
          l1_runoffseal, &
          l1_fastrunoff, &
          l1_slowrunoff, &
          l1_baseflow, &
          l1_percol, &
          l1_infilsoil, &
          l1_preeffect, &
          l1_melt &
        )
 
        ! write data
        if (run_cfg%domainDateTime%writeout(timestep_model_outputs, tt)) then
          call run_cfg%nc_mhm%writeTimestep(run_cfg%domainDateTime%tIndex_out * timestep)
        end if
 
        if(tt == run_cfg%domainDateTime%nTimeSteps) then
          call run_cfg%nc_mhm%close()
        end if
 
      end if
    end if ! <-- if (.not. optimize)
 
  end subroutine mhm_interface_run_write_output
 
  !> \brief add simulation data to optimization data types
  subroutine mhm_interface_run_update_optisim(etOptiSim, twsOptiSim, neutronsOptiSim, smOptiSim)
    implicit none
    !> returns soil moisture time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
    type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim
    !> dim1=ncells, dim2=time
    type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim
    !> returns evapotranspiration time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
    type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim
    !> returns tws time series for all grid cells (of multiple Domains concatenated),DIMENSION [nCells, nTimeSteps]
    type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim
 
    integer(i4) :: gg, s1, e1, lcId
 
    integer(i4) :: iDomain, tt
 
    ! get time step
    tt = run_cfg%time_step
 
    ! get domain index
    idomain = run_cfg%get_domain_index(run_cfg%selected_domain)
 
    ! slice settings
    s1 = run_cfg%s1
    e1 = run_cfg%e1
    lcid = run_cfg%domainDateTime%yId
 
    !----------------------------------------------------------------------
    ! FOR SOIL MOISTURE
    ! NOTE:: modeled soil moisture is averaged according to input time step
    !        soil moisture (timeStep_sm_input)
    !----------------------------------------------------------------------
    if (present(smoptisim)) then
      ! only for evaluation period - ignore warming days
      if ((tt - warmingdays(idomain) * ntstepday) .GT. 0) then
        ! decide for daily, monthly or yearly aggregation
        call smoptisim(idomain)%average_per_timestep(l1_smobs(idomain)%timeStepInput, &
                     run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
        ! last timestep is already done - write_counter exceeds size(smOptiSim(iDomain)%dataSim, dim=2)
        if (tt /= run_cfg%domainDateTime%nTimeSteps) then
          ! aggregate soil moisture to needed time step for optimization
          call smoptisim(idomain)%average_add(sum(l1_soilmoist(:, 1 : nsoilhorizons_sm_input), dim = 2) / &
                          sum(l1_soilmoistsat(:, 1 : nsoilhorizons_sm_input, lcid), dim = 2))
        end if
      end if
    end if
 
    !----------------------------------------------------------------------
    ! FOR NEUTRONS
    ! NOTE:: modeled neutrons are averaged daily
    !----------------------------------------------------------------------
    if (present(neutronsoptisim)) then
      ! only for evaluation period - ignore warming days
      if ((tt - warmingdays(idomain) * ntstepday) .GT. 0) then
        ! decide for daily, monthly or yearly aggregation
        ! daily
        if (run_cfg%domainDateTime%is_new_day)   then
          call neutronsoptisim(idomain)%average()
        end if
 
        ! last timestep is already done - write_counter exceeds size(sm_opti, dim=2)
        if (tt /= run_cfg%domainDateTime%nTimeSteps) then
          ! aggregate neutrons to needed time step for optimization
          call neutronsoptisim(idomain)%average_add(l1_neutrons(s1 : e1))
        end if
      end if
    end if
 
    !----------------------------------------------------------------------
    ! FOR EVAPOTRANSPIRATION
    ! NOTE:: modeled evapotranspiration is averaged according to input time step
    !        evapotranspiration (timeStep_et_input)
    !----------------------------------------------------------------------
    if (present(etoptisim)) then
      ! only for evaluation period - ignore warming days
      if ((tt - warmingdays(idomain) * ntstepday) .GT. 0) then
        ! decide for daily, monthly or yearly aggregation
        call etoptisim(idomain)%increment_counter(l1_etobs(idomain)%timeStepInput, &
          run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
 
        ! last timestep is already done - write_counter exceeds size(etOptiSim(iDomain)%dataSim, dim=2)
        if (tt /= run_cfg%domainDateTime%nTimeSteps) then
          ! aggregate evapotranspiration to needed time step for optimization
          call etoptisim(idomain)%add(sum(l1_aetsoil(s1 : e1, :), dim = 2) * &
            run_cfg%L1_fNotSealed(s1 : e1, 1, lcid) + &
            l1_aetcanopy(s1 : e1) + &
            l1_aetsealed(s1 : e1) * l1_fsealed(s1 : e1, 1, lcid))
        end if
      end if
    end if
 
    !----------------------------------------------------------------------
    ! FOR TWS
    ! NOTE:: modeled tws is averaged according to input time step
    !        (timeStepInput)
    !----------------------------------------------------------------------
    if (present(twsoptisim)) then
      ! only for evaluation period - ignore warming days
      if ((tt - warmingdays(idomain) * ntstepday) > 0) then
        ! decide for daily, monthly or yearly aggregation
        call twsoptisim(idomain)%average_per_timestep(l1_twsaobs(idomain)%timeStepInput, &
          run_cfg%domainDateTime%is_new_day, run_cfg%domainDateTime%is_new_month, run_cfg%domainDateTime%is_new_year)
 
        ! last timestep is already done - write_counter exceeds size(twsOptiSim(iDomain)%dataSim, dim=2)
        if (tt /= run_cfg%domainDateTime%nTimeSteps) then
          ! aggregate evapotranspiration to needed time step for optimization
          call twsoptisim(idomain)%average_add(l1_inter(s1 : e1) + l1_snowpack(s1 : e1) + l1_sealstw(s1 : e1) + &
            l1_unsatstw(s1 : e1) + l1_satstw(s1 : e1))
          do gg = 1, nsoilhorizons_mhm
            call twsoptisim(idomain)%add(l1_soilmoist(s1 : e1, gg))
          end do
        end if
      end if
    end if
 
    ! TODO-RIV-TEMP: add OptiSim for river temperature
 
  end subroutine mhm_interface_run_update_optisim
 
  !> \brief finalize current domain after running
  subroutine mhm_interface_run_finalize_domain()
    implicit none
 
    integer(i4) :: iDomain
 
    ! get domain index
    idomain = run_cfg%get_domain_index(run_cfg%selected_domain)
 
    if (allocated(run_cfg%InflowDischarge)) deallocate(run_cfg%InflowDischarge)
      if (domainmeta%doRouting(idomain)) then
      ! clean runoff variable
      deallocate(run_cfg%RunToRout)
      if ( riv_temp_pcs%active ) call riv_temp_pcs%dealloc_lateral()
    end if
 
  end subroutine mhm_interface_run_finalize_domain
 
  !> \brief finalize run
  subroutine mhm_interface_run_finalize(runoff, BFI)
    implicit none
    !> returns runoff time series, DIMENSION [nTimeSteps, nGaugesTotal]
    real(dp), dimension(:, :), allocatable, optional, intent(out) :: runoff
    !> baseflow index, dim1=domainID
    real(dp), dimension(:), allocatable, optional, intent(out) :: BFI
 
    if (present(runoff) .and. (processmatrix(8, 1) > 0)) runoff = mrm_runoff
    if (present(bfi)) then
      bfi = bfi_qbf_sum / bfi_qt_sum
      deallocate(bfi_qbf_sum)
      deallocate(bfi_qt_sum)
    end if
 
    if (allocated(run_cfg%parameterset)) deallocate(run_cfg%parameterset)
    if (allocated(run_cfg%L1_fNotSealed)) deallocate(run_cfg%L1_fNotSealed)
    if (allocated(run_cfg%domain_indices)) deallocate(run_cfg%domain_indices)
 
  end subroutine mhm_interface_run_finalize
 
end module mo_mhm_interface_run