mo_mrm_init.f90

Go to the documentation of this file.

!> \file mo_mrm_init.f90
!> \brief \copybrief mo_mrm_init
!> \details \copydetails mo_mrm_init
 
!> \brief Wrapper for initializing Routing.
!> \details Calling all routines to initialize all mRM variables
!> \authors Luis Samaniego, Rohini Kumar and Stephan Thober
!> \date Aug 2015
!> \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_mrm
MODULE mo_mrm_init
 
    use mo_common_variables, only : dirout
    use mo_message, only : message, error_message
 
  ! This module sets the river network characteristics and routing order.
 
  ! Written  Luis Samaniego, Mar 2005
 
  IMPLICIT NONE
 
  public :: mrm_init, mrm_configuration
  public :: variables_default_init_routing
  public :: fluxes_states_default_init_routing
 
  private
 
CONTAINS
 
  !> \brief read mRM configuration from namelists
  subroutine mrm_configuration(file_namelist, file_namelist_param)
    use mo_common_mhm_mrm_variables, only : mrm_coupling_mode, optimize
    use mo_common_variables, only : processmatrix
    use mo_mrm_read_config, only : mrm_read_config
    use mo_mrm_global_variables, only: riv_temp_pcs
    use mo_common_read_config, only : common_read_config
    use mo_common_mhm_mrm_read_config, only : check_optimization_settings, common_mhm_mrm_read_config
    use mo_kind, only : i4
    implicit none
 
    character(*), intent(in) :: file_namelist !< namelist file name
    character(*), intent(in) :: file_namelist_param !< parameter namelist file name
 
    if (mrm_coupling_mode .eq. 0_i4) then
      call common_read_config(file_namelist)
      call common_mhm_mrm_read_config(file_namelist)
      !-----------------------------------------------------------
      ! PRINT STARTUP MESSAGE
      !-----------------------------------------------------------
      call print_startup_message(file_namelist, file_namelist_param)
    else
      call message('')
      call message('  Inititalize mRM')
      if ( processmatrix(11, 1) .ne. 0 ) then
        ! processCase(11): river temperature routing
        riv_temp_pcs%active = .true.
        riv_temp_pcs%case = processmatrix(11, 1)
        call message('')
        call message('    Read config: river temperature routing')
        call riv_temp_pcs%config(file_namelist, file_namelist_param)
      end if
    end if
 
    ! read config for mrm, readlatlon is set here depending on whether output is needed
    call mrm_read_config(file_namelist, file_namelist_param, (mrm_coupling_mode .eq. 0_i4))
 
    ! this was moved here, because it depends on global_parameters that are only set in mrm_read_config
    if (mrm_coupling_mode .eq. 0_i4) then
      if (optimize) call check_optimization_settings()
      !-----------------------------------------------------------
      ! CONFIG OUTPUT
      !-----------------------------------------------------------
      call config_output()
    end if
  end subroutine mrm_configuration
 
 
  !> \brief Initialize all mRM variables at all levels (i.e., L0, L1, and L11).
  !> \details Initialize all mRM variables at all levels (i.e., L0, L1, and L11)
  !! either with default values or with values from restart file. The L0 mask (L0_mask),
  !! L0 elevation (L0_elev), and L0 land cover (L0_LCover) can be provided as optional
  !! variables to save memory because these variable will then not be read in again.
  !> \changelog
  !! - Stephan Thober Sep 2015
  !!   - added L0_mask, L0_elev, and L0_LCover
  !! - Stephan Thober May 2016
  !!   - added warning message in case no gauge is found in modelling domain
  !! - Matthias Kelbling Aug 2017
  !!   - added L11_flow_accumulation to Initialize Stream Netwo
  !! - Lennart Schueler May 2018
  !!   - added initialization for groundwater coupling
  !! - Stephan Thober Jun 2018
  !!   - refactored for mpr_extract version
  !! - Stephan Thober May 2019
  !!   - added init of level0 in case of read restart
  !> \authors Stephan Thober
  !> \date Aug 2015
  subroutine mrm_init(file_namelist, file_namelist_param)
 
    use mo_common_constants, only : nodata_dp, nodata_i4
    use mo_common_mhm_mrm_variables, only : mrmfilerestartin, mrm_coupling_mode, mrm_read_river_network, &
                                            resolutionrouting
    use mo_common_restart, only : read_grid_info
    use mo_common_variables, only : global_parameters, l0_l1_remap, level0, level1, domainmeta, &
                                    processmatrix, resolutionhydrology
    use mo_common_grid, only : l0_grid_setup, init_lowres_level, set_domain_indices
    use mo_kind, only : i4
    use mo_mrm_global_variables, only : domain_mrm, &
                                        l0_l11_remap, level11, &
                                        gw_coupling, l0_river_head_mon_sum, &
                                        l11_netperm, l11_fromn, l11_length, l11_noutlets, &
                                        riv_temp_pcs, &
                                        readlatlon, &
                                        sink_cells
    use mo_mrm_net_startup, only : l11_flow_direction, l11_flow_accumulation, l11_fraction_sealed_floodplain, &
                                   l11_link_location, l11_routing_order, l11_set_drain_outlet_gauges, &
                                   l11_set_network_topology, l11_stream_features, l11_l1_mapping
    use mo_mrm_read_data, only : mrm_read_l0_data, mrm_read_discharge, &
                                 mrm_read_total_runoff, mrm_read_bankfull_runoff
    use mo_mrm_restart, only : mrm_read_restart_config
    use mo_read_latlon, only : read_latlon
    use mo_mrm_river_head, only: init_masked_zeros_l0, calc_channel_elevation
    use mo_mrm_mpr, only : mrm_init_param
    use mo_timer, only : timer_get, timer_start, timer_stop, timer_clear
    use mo_string_utils, only : num2str
 
    implicit none
 
    character(*), intent(in) :: file_namelist !< namelist file name
    character(*), intent(in) :: file_namelist_param !< parameter namelist file name
 
    ! start and end index for routing parameters
    integer(i4) :: istart, iend
    ! start and end index at L11
    integer(i4) :: s11, e11
 
    integer(i4) :: domainid, idomain, gauge_counter
 
 
    if (mrm_coupling_mode .eq. 0_i4) then
      allocate(l0_l1_remap(domainmeta%nDomains))
      allocate(level1(domainmeta%nDomains))
    end if
 
    ! ----------------------------------------------------------
    ! READ DATA
    ! ----------------------------------------------------------
    allocate(level11(domainmeta%nDomains))
    allocate(l0_l11_remap(domainmeta%nDomains))
    allocate(sink_cells(domainmeta%nDomains))
 
    if (.not. mrm_read_river_network) then
      ! read all (still) necessary level 0 data
      if (processmatrix(8, 1) .eq. 1_i4) call mrm_read_l0_data(mrm_coupling_mode .eq. 0_i4, readlatlon, .true.)
      if (processmatrix(8, 1) .eq. 2_i4) call mrm_read_l0_data(mrm_coupling_mode .eq. 0_i4, readlatlon, .false.)
      if (processmatrix(8, 1) .eq. 3_i4) call mrm_read_l0_data(mrm_coupling_mode .eq. 0_i4, readlatlon, .false.)
    end if
 
    do idomain = 1, domainmeta%nDomains
      allocate(sink_cells(idomain)%ids(0))
      domainid = domainmeta%indices(idomain)
      if (mrm_read_river_network) then
        ! this reads the domain properties
        if (.not. allocated(level0)) allocate(level0(domainmeta%nDomains))
        ! ToDo: L0_Domain, parallel
        call read_grid_info(mrmfilerestartin(idomain), "0", level0(domainmeta%L0DataFrom(idomain)))
        if (mrm_coupling_mode .eq. 0_i4) then
          call read_grid_info(mrmfilerestartin(idomain), "1", level1(idomain))
        end if
        call read_grid_info(mrmfilerestartin(idomain), "11", level11(idomain))
        call mrm_read_restart_config(idomain, domainid, mrmfilerestartin(idomain))
      else
        if (idomain .eq. 1) then
          call l0_check_input_routing(domainmeta%L0DataFrom(idomain))
          if (mrm_coupling_mode .eq. 0_i4) then
            call l0_grid_setup(level0(domainmeta%L0DataFrom(idomain)))
          end if
        else if ((domainmeta%L0DataFrom(idomain) == idomain)) then
          call l0_check_input_routing(domainmeta%L0DataFrom(idomain))
          if (mrm_coupling_mode .eq. 0_i4) then
            call l0_grid_setup(level0(domainmeta%L0DataFrom(idomain)))
          end if
        end if
 
        if (mrm_coupling_mode .eq. 0_i4) then
          call init_lowres_level(level0(domainmeta%L0DataFrom(idomain)), resolutionhydrology(idomain), &
                  level1(idomain), l0_l1_remap(idomain))
        end if
        call init_lowres_level(level0(domainmeta%L0DataFrom(idomain)), resolutionrouting(idomain), &
                level11(idomain), l0_l11_remap(idomain))
        call l11_l1_mapping(idomain)
 
        if (readlatlon) then
          ! read lat lon coordinates of each domain
          call read_latlon(idomain, "lon", "lat", "level1", level1(idomain))
          call read_latlon(idomain, "lon_l11", "lat_l11", "level11", level11(idomain))
        else
          ! allocate the memory and set to nodata
          allocate(level11(idomain)%x(level11(idomain)%nrows, level11(idomain)%ncols))
          allocate(level11(idomain)%y(level11(idomain)%nrows, level11(idomain)%ncols))
          level11(idomain)%x = nodata_dp
          level11(idomain)%y = nodata_dp
        end if
      end if
    end do
 
    call set_domain_indices(level11)
    call set_domain_indices(level1)
    call set_domain_indices(level0, indices=domainmeta%L0DataFrom)
 
    ! ----------------------------------------------------------
    ! INITIALIZE STATES AND AUXILLIARY VARIABLES
    ! ----------------------------------------------------------
    do idomain = 1, domainmeta%nDomains
      call variables_alloc_routing(idomain)
    end do
 
    ! ----------------------------------------------------------
    ! INITIALIZE STREAM NETWORK
    ! ----------------------------------------------------------
    do idomain = 1, domainmeta%nDomains
       if (.not. mrm_read_river_network) then
        
        call timer_clear(1)
        call timer_start(1)
        call l11_flow_direction(idomain)
        call message(' ')
        call message('    Flow direction upscaled      ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        
        call timer_clear(1)
        call timer_start(1)
        call l11_flow_accumulation(idomain)
        call message('    Flow accumulation upscaled   ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        
        call timer_clear(1)
        call timer_start(1)
        call l11_set_network_topology(idomain)
        call message('    Topology configured          ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        
        call timer_clear(1)
        call timer_start(1)
        call l11_routing_order(idomain)
        call message('    Routing order ready          ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        
        call timer_clear(1)
        call timer_start(1)
        call l11_link_location(idomain)
        call message('    Link location done           ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        
        call timer_clear(1)
        call timer_start(1)
        call l11_set_drain_outlet_gauges(idomain)
        call message('    Gauges assigned with nodes   ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        
        ! stream characteristics
        call timer_clear(1)
        call timer_start(1)
        call l11_stream_features(idomain)
        call message('    Stream features generated    ...')
        call timer_stop(1)
        call message('    in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
        call timer_clear(1)
      end if
    end do
 
    ! ----------------------------------------------------------
    ! INITIALIZE PARAMETERS
    ! ----------------------------------------------------------
    do idomain = 1, domainmeta%nDomains
      istart = processmatrix(8, 3) - processmatrix(8, 2) + 1
      iend = processmatrix(8, 3)
      call mrm_init_param(idomain, global_parameters(istart : iend, 3))
    end do
 
    ! check whether there are gauges within the modelling domain
    if (allocated(domain_mrm)) then
      gauge_counter = 0
      do idomain = 1, domainmeta%nDomains
        if (.not. all(domain_mrm(idomain)%gaugeNodeList .eq. nodata_i4)) then
          gauge_counter = gauge_counter + 1
        end if
      end do
      if (gauge_counter .lt. 1) then
        call message('')
        call message('    WARNING: no gauge found within modelling domain')
      end if
    end if
    ! mpr-like definiton of sealed floodplain fraction
    if ((processmatrix(8, 1) .eq. 1_i4) .and. (.not. mrm_read_river_network)) then
      call l11_fraction_sealed_floodplain(2_i4, .true.)
    else
      ! dummy initialization
      call l11_fraction_sealed_floodplain(2_i4, .false.)
    end if
 
    ! -------------------------------------------------------
    ! READ INPUT DATA AND OBSERVED DISCHARGE DATA
    ! -------------------------------------------------------
    ! read simulated runoff at level 1
    if (mrm_coupling_mode .eq. 0_i4) then
      do idomain = 1, domainmeta%nDomains
        call mrm_read_total_runoff(idomain)
      end do
    end if
    ! discharge data
    call mrm_read_discharge()
 
    ! init groundwater coupling
    if (gw_coupling) then
      do idomain = 1, domainmeta%nDomains
        call init_masked_zeros_l0(idomain, l0_river_head_mon_sum)
        call mrm_read_bankfull_runoff(idomain)
      end do
      call calc_channel_elevation()
    end if
 
    ! init riv temp
    if ( riv_temp_pcs%active ) then
      call message('')
      call message('    Initialization of river temperature routing.')
      do idomain = 1, domainmeta%nDomains
        s11 = level11(idomain)%iStart
        e11 = level11(idomain)%iEnd
        call riv_temp_pcs%init(level11(idomain)%nCells)
        call riv_temp_pcs%init_area( &
          idomain, &
          l11_netperm(s11 : e11), & ! routing order at L11
          l11_fromn(s11 : e11), & ! link source at L11
          l11_length(s11 : e11 - 1), & ! link length
          level11(idomain)%nCells - l11_noutlets(idomain), &
          level11(idomain)%nCells, &
          level11(idomain)%nrows, &
          level11(idomain)%ncols, &
          level11(idomain)%mask &
        )
      end do
    end if
    call message('')
    call message('  Finished Initialization of mRM')
 
  end subroutine mrm_init
 
 
  !> \brief Print mRM startup message
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine print_startup_message(file_namelist, file_namelist_param)
 
    use mo_kind, only : i4
    use mo_mrm_file, only : file_defoutput, file_main, version, version_date
    use mo_string_utils, only : num2str, separator
 
    implicit none
 
    character(*), intent(in) :: file_namelist !< namelist file name
    character(*), intent(in) :: file_namelist_param !< parameter namelist file name
 
    ! Date and time
    integer(i4), dimension(8) :: datetime
 
    CHARACTER(len=1024) :: message_text = ''
 
    call message(separator)
    call message('              mRM-UFZ')
    call message()
    call message('    MULTISCALE ROUTING MODEL')
    call message('           Version ', trim(version))
    call message('           ', trim(version_date))
    call message()
    call message('Made available by S. Thober & M. Cuntz')
    call message()
    call message('Based on mHM-UFZ by L. Samaniego & R. Kumar')
 
    call message(separator)
 
    call message()
    call date_and_time(values = datetime)
    message_text = trim(num2str(datetime(3), '(I2.2)')) // "." // trim(num2str(datetime(2), '(I2.2)')) &
            // "." // trim(num2str(datetime(1), '(I4.4)')) // " " // trim(num2str(datetime(5), '(I2.2)')) &
            // ":" // trim(num2str(datetime(6), '(I2.2)')) // ":" // trim(num2str(datetime(7), '(I2.2)'))
    call message('Start at ', trim(message_text), '.')
    call message('Using main file ', trim(file_main), ' and namelists: ')
    call message('     ', trim(file_namelist))
    call message('     ', trim(file_namelist_param))
    call message('     ', trim(file_defoutput), ' (if it is given)')
    call message()
 
  end subroutine print_startup_message
 
 
  !> \brief print mRM configuration
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine config_output
 
    use mo_common_variables, only : dirlcover, dirmorpho, dirout, domainmeta
    use mo_kind, only : i4
    use mo_mrm_file, only : file_defoutput, file_namelist_mrm, file_namelist_param_mrm
    use mo_mrm_global_variables, only : domain_mrm, &
                                        dirgauges
    use mo_string_utils, only : num2str
 
    implicit none
 
    integer(i4) :: domainID, iDomain
 
    integer(i4) :: jj
 
 
    !
    call message()
    call message('Read namelist file: ', trim(file_namelist_mrm))
    call message('Read namelist file: ', trim(file_namelist_param_mrm))
    call message('Read namelist file: ', trim(file_defoutput), ' (if it is given)')
 
    call message()
    call message('  # of domains:         ', trim(num2str(domainmeta%nDomains)))
    call message()
    call message('  Input data directories:')
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      call message('  --------------')
      call message('      DOMAIN                   ', num2str(domainid, '(I3)'))
      call message('  --------------')
      call message('    Morphological directory:    ', trim(dirmorpho(idomain)))
      call message('    Land cover directory:       ', trim(dirlcover(idomain)))
      call message('    Discharge directory:        ', trim(dirgauges(idomain)))
      call message('    Output directory:           ', trim(dirout(idomain)))
      call message('    Evaluation gauge            ', 'ID')
      do jj = 1, domain_mrm(idomain)%nGauges
        call message('    ', trim(adjustl(num2str(jj))), '                           ', &
                trim(adjustl(num2str(domain_mrm(idomain)%gaugeIdList(jj)))))
      end do
      if (domain_mrm(idomain)%nInflowGauges .GT. 0) then
        call message('    Inflow gauge              ', 'ID')
        do jj = 1, domain_mrm(idomain)%nInflowGauges
          call message('    ', trim(adjustl(num2str(jj))), '                         ', &
                  trim(adjustl(num2str(domain_mrm(idomain)%InflowGaugeIdList(jj)))))
        end do
      end if
    end do
  end subroutine config_output
 
 
  !> \brief Default initalization mRM related L11 variables
  !> \details Default initalization of mHM related L11 variables (e.g., states,
  !! fluxes, and parameters) as per given constant values given in mo_mhm_constants.
  !! Variables initalized here is defined in the mo_global_variables.f90 file.
  !! Only Variables that are defined in the variables_alloc subroutine are
  !! intialized here.
  !! If a variable is added or removed here, then it also has to be added or removed
  !! in the subroutine state_variables_set in the module mo_restart and in the
  !! subroutine set_state in the module mo_set_netcdf_restart.
  !> \authors  Stephan Thober, Rohini Kumar, and Juliane Mai
  !> \date    Aug 2015
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine variables_default_init_routing
 
    use mo_common_constants, only : p1_initstatefluxes
    use mo_mrm_global_variables, only : l11_c1, l11_c2, l11_k, l11_xi
 
    implicit none
 
    !-------------------------------------------
    ! L11 ROUTING STATE VARIABLES, FLUXES AND
    !             PARAMETERS
    !-------------------------------------------
 
    ! fluxes and states
    call fluxes_states_default_init_routing()
 
    ! kappa: Muskingum travel time parameter.
    l11_k = p1_initstatefluxes
    ! xi:    Muskingum diffusion parameter
    l11_xi = p1_initstatefluxes
    ! Routing parameter C1=f(K,xi, DT) (Chow, 25-41)
    l11_c1 = p1_initstatefluxes
    ! Routing parameter C2 =f(K,xi, DT) (Chow, 25-41)
    l11_c2 = p1_initstatefluxes
 
  end subroutine variables_default_init_routing
 
  !> \brief initialize fluxes and states with default values for mRM
  subroutine fluxes_states_default_init_routing(iDomain)
 
    use mo_kind, only: i4
    use mo_mrm_global_variables, only : level11
    use mo_common_constants, only : p1_initstatefluxes
    use mo_mrm_global_variables, only : l11_qmod, l11_qout, l11_qtin, l11_qtr
 
    implicit none
 
    !> number of Domain (if not present, set for all)
    integer(i4), intent(in), optional :: idomain
 
    integer(i4) :: s11, e11
 
    !-------------------------------------------
    ! L11 ROUTING STATE VARIABLES, FLUXES AND
    !             PARAMETERS
    !-------------------------------------------
 
    if (present(idomain)) then
      s11 = level11(idomain)%iStart
      e11 = level11(idomain)%iEnd
      ! simulated discharge at each node
      l11_qmod(s11 : e11) = p1_initstatefluxes
      ! Total outflow from cells L11 at time tt
      l11_qout(s11 : e11) = p1_initstatefluxes
      ! Total discharge inputs at t-1 and t
      l11_qtin(s11 : e11, :) = p1_initstatefluxes
      !  Routed outflow leaving a node
      l11_qtr(s11 : e11, :) = p1_initstatefluxes
    else
      ! simulated discharge at each node
      l11_qmod = p1_initstatefluxes
      ! Total outflow from cells L11 at time tt
      l11_qout = p1_initstatefluxes
      ! Total discharge inputs at t-1 and t
      l11_qtin = p1_initstatefluxes
      !  Routed outflow leaving a node
      l11_qtr = p1_initstatefluxes
    end if
 
  end subroutine fluxes_states_default_init_routing
 
 
  !> \brief check routing input on level-0
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine l0_check_input_routing(L0Domain_iDomain)
 
    use mo_common_constants, only : nodata_i4
    use mo_common_variables, only : level0
    use mo_kind, only : i4
    use mo_mrm_global_variables, only : l0_facc, l0_fdir
    use mo_string_utils, only : num2str
 
    implicit none
 
    integer(i4), intent(in) :: L0Domain_iDomain !< domain index for associated level-0 data
 
    integer(i4) :: k
 
    CHARACTER(len=1024) :: message_text = ''
 
    do k = level0(l0domain_idomain)%iStart, level0(l0domain_idomain)%iEnd
      ! flow direction [-]
      if (l0_fdir(k) .eq. nodata_i4) then
        message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(l0domain_idomain, '(I5)'))
        call error_message(' Error: flow direction has missing value within the valid masked area at cell in domain ', &
                trim(message_text))
      end if
      ! flow accumulation [-]
      if (l0_facc(k) .eq. nodata_i4) then
        message_text = trim(num2str(k, '(I5)')) // ',' // trim(num2str(l0domain_idomain, '(I5)'))
        call error_message(' Error: flow accumulation has missing values within the valid masked area at cell in domain ', &
                trim(message_text))
      end if
    end do
 
  end subroutine l0_check_input_routing
 
 
  !> \brief allocated routing related variables
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine variables_alloc_routing(iDomain)
 
    use mo_append, only : append
    use mo_kind, only : dp, i4
    use mo_mrm_constants, only : nroutingstates
    use mo_common_variables, only : level0, domainmeta
    use mo_mrm_global_variables, only : l11_c1, l11_c2, l11_k, &
         l11_qmod, l11_qout, l11_qtin, l11_qtr, l11_xi, &
         level11, l11_celerity, l0_celerity
 
    implicit none
 
    integer(i4), intent(in) :: iDomain !< domain index
 
    real(dp), dimension(:), allocatable :: dummy_Vector11
 
    real(dp), dimension(:, :), allocatable :: dummy_Matrix11_IT
 
 
    ! dummy vector and matrix
    allocate(dummy_vector11(level11(idomain)%nCells))
    allocate(dummy_matrix11_it(level11(idomain)%nCells, nroutingstates))
 
    ! simulated discharge at each node
    dummy_vector11(:) = 0.0_dp
    call append(l11_qmod, dummy_vector11)
 
    ! Total outflow from cells L11 at time tt
    dummy_vector11(:) = 0.0_dp
    call append(l11_qout, dummy_vector11)
 
    ! Total discharge inputs at t-1 and t
    dummy_matrix11_it(:, :) = 0.0_dp
    call append(l11_qtin, dummy_matrix11_it)
 
    !  Routed outflow leaving a node
    dummy_matrix11_it(:, :) = 0.0_dp
    call append(l11_qtr, dummy_matrix11_it)
 
    ! kappa: Muskingum travel time parameter.
    dummy_vector11(:) = 0.0_dp
    call append(l11_k, dummy_vector11)
 
    ! xi:    Muskingum diffusion parameter
    dummy_vector11(:) = 0.0_dp
    call append(l11_xi, dummy_vector11)
 
    ! Routing parameter C1=f(K,xi, DT) (Chow, 25-41)
    dummy_vector11(:) = 0.0_dp
    call append(l11_c1, dummy_vector11)
 
    ! Routing parameter C2 =f(K,xi, DT) (Chow, 25-41)
    dummy_vector11(:) = 0.0_dp
    call append(l11_c2, dummy_vector11)
 
    ! Celerity at each link
    dummy_vector11(:) = 0.0_dp
    call append(l11_celerity, dummy_vector11)
 
    ! celerity at level 0
    if (allocated(dummy_vector11)) deallocate(dummy_vector11)
    allocate(dummy_vector11(level0(domainmeta%L0DataFrom(idomain))%ncells))
    dummy_vector11(:) = 0.0_dp
    call append(l0_celerity, dummy_vector11)
 
    ! free space
    if (allocated(dummy_vector11)) deallocate(dummy_vector11)
    if (allocated(dummy_matrix11_it)) deallocate(dummy_matrix11_it)
 
  end subroutine variables_alloc_routing
 
END MODULE mo_mrm_init