mo_mrm_riv_temp_class.f90

Go to the documentation of this file.

!> \file    mo_mrm_riv_temp_class.f90
!> \brief   \copybrief mo_mrm_riv_temp_class
!> \details \copydetails mo_mrm_riv_temp_class
 
!> \brief   Class for the river temperature calculations
!> \details River temperature routing on top of mRM.
!> \warning This feature is still experimental! River freezing is still missing.
!> \version 0.1
!> \authors Sebastian Mueller
!> \date    Sep 2020
!> \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_riv_temp_class
 
  use mo_kind, only: dp, i4
 
  implicit none
 
  private
 
  !> \class   riv_temp_type
  !> \brief   This is a container to define the river temperature routing in the current time step
  !> \details This class provides all procedures to rout river tmperature through the river network.
  !> \warning This feature is still experimental!
  !!          This first version doesn't provide ice covering, which means, that the river temperature
  !!          can drop below 0 [deg C] in winter.
  type, public :: riv_temp_type
    ! config settings
    logical :: active = .false. !< state if this process is active
    integer(i4) :: case = 0_i4 !< the selected process-case option
    character(256) :: nml_name = 'config_riv_temp' !< namelist name in mhm namelist
    ! riv geometry
    character(256), dimension(:), allocatable :: dir_riv_widths !< Directory where river widths are stored
    character(256) :: riv_widths_file !< file name for river widths
    character(256) :: riv_widths_name !< variable name for river widths
    real(dp), public, dimension(:), allocatable :: l11_riv_widths !< river widths in L11
    real(dp), public, dimension(:), allocatable :: l11_riv_areas !< river area in L11
    ! PET related vars
    real(dp) :: albedo_water !< albedo of open water
    real(dp) :: pt_a_water   !< priestley taylor alpha parameter for PET on open water
    ! sensible heat flux related
    real(dp) :: emissivity_water   !< emissivity of water
    real(dp) :: turb_heat_ex_coeff !< lateral heat exchange coefficient water <-> air
    !> cutoff value for temperature
    real(dp) :: delta_t = 0.1_dp
    ! controlling variables for iterative solver
    integer(i4) :: max_iter  !< input: maximal number of iterations done
    real(dp) :: delta_iter   !< input: convergence criteria for iterative solver
    real(dp) :: step_iter    !< input: step-size for linear search
    logical :: first_iter    !< whether it is at the first iteration (to determine search direction)
    logical :: up_iter       !< whether the search direction is upwards
    logical :: bisect_iter   !< whether to do the bisection search part (after the interval is found)
    real(dp) :: up_bnd_iter  !< upper bound for the current bisection step
    real(dp) :: low_bnd_iter !< lower bound for the current bisection step
    ! accumulated later fluxes (in current time-step)
    real(dp), dimension(:), allocatable :: l1_runoff_e !< runoff energy at L1 level
    real(dp), dimension(:), allocatable :: l1_acc_ssrd !< accumulated shortwave radiation at L1 level
    real(dp), dimension(:), allocatable :: l1_acc_strd !< accumulated longwave radiation at L1 level
    real(dp), dimension(:), allocatable :: l1_acc_temp !< accumulated air temperature at L1 level
    real(dp), dimension(:), allocatable :: l1_ssrd_calc !< current shortwave radiation at L1 level
    real(dp), dimension(:), allocatable :: l1_strd_calc !< current longwave radiation at L1 level
    real(dp), dimension(:), allocatable :: l1_tann_calc !< current mean air temperature at L1 level
    integer(i4) :: ts_cnt !< sub time-step counter for accumulation of meteo
    ! vars for routing
    integer(i4) :: s11 !< starting index for current L11 domain
    integer(i4) :: e11 !< ending index for current L11 domain
    real(dp), dimension(:,:), allocatable :: netnode_e_in !< Total energy inputs at t-1 and t
    real(dp), dimension(:,:), allocatable :: netnode_e_r  !< energy leaving at t-1 and t
    real(dp), dimension(:), allocatable :: netnode_e_mod  !< Simulated routed energy
    real(dp), dimension(:), allocatable :: netnode_e_out  !< total energy source from cell in L11
    real(dp), dimension(:), allocatable :: l11_srad_net   !< net short wave radiation at L11
    real(dp), dimension(:), allocatable :: l11_lrad_in    !< incoming long wave radiation at L11
    real(dp), dimension(:), allocatable :: l11_air_temp   !< air temp at L11
    ! variable containing river temp for each domain at L11 level
    real(dp), dimension(:), allocatable :: river_temp !< resulting river temp at L11 in [deg C]
  contains
    ! config and inits
    !> \copydoc mo_mrm_riv_temp_class::config
    procedure :: config !< \see mo_mrm_riv_temp_class::config
    !> \copydoc mo_mrm_riv_temp_class::init
    procedure :: init !< \see mo_mrm_riv_temp_class::init
    !> \copydoc mo_mrm_riv_temp_class::init_area
    procedure :: init_area !< \see mo_mrm_riv_temp_class::init_area
    !> \copydoc mo_mrm_riv_temp_class::init_riv_temp
    procedure :: init_riv_temp !< \see mo_mrm_riv_temp_class::init_riv_temp
 
    ! source accumulations
    !> \copydoc mo_mrm_riv_temp_class::acc_source_e
    procedure :: acc_source_e !< \see mo_mrm_riv_temp_class::acc_source_e
    !> \copydoc mo_mrm_riv_temp_class::finalize_source_e
    procedure :: finalize_source_e !< \see mo_mrm_riv_temp_class::finalize_source_e
 
    ! temp-energy routing routines
    !> \copydoc mo_mrm_riv_temp_class::get_lrad_out
    procedure :: get_lrad_out !< \see mo_mrm_riv_temp_class::get_lrad_out
    !> \copydoc mo_mrm_riv_temp_class::get_lat_heat
    procedure :: get_lat_heat !< \see mo_mrm_riv_temp_class::get_lat_heat
    !> \copydoc mo_mrm_riv_temp_class::get_sens_heat
    procedure :: get_sens_heat !< \see mo_mrm_riv_temp_class::get_sens_heat
    !> \copydoc mo_mrm_riv_temp_class::get_e_io
    procedure :: get_e_io !< \see mo_mrm_riv_temp_class::get_e_io
    !> \copydoc mo_mrm_riv_temp_class::l11_routing_e
    procedure :: l11_routing_e !< \see mo_mrm_riv_temp_class::l11_routing_e
 
    ! helper for iterative solver
    !> \copydoc mo_mrm_riv_temp_class::init_iter
    procedure :: init_iter !< \see mo_mrm_riv_temp_class::init_iter
    !> \copydoc mo_mrm_riv_temp_class::next_iter
    procedure :: next_iter !< \see mo_mrm_riv_temp_class::next_iter
 
    ! care taker
    !> \copydoc mo_mrm_riv_temp_class::reset_timestep
    procedure :: reset_timestep !< \see mo_mrm_riv_temp_class::reset_timestep
    !> \copydoc mo_mrm_riv_temp_class::alloc_lateral
    procedure :: alloc_lateral !< \see mo_mrm_riv_temp_class::alloc_lateral
    !> \copydoc mo_mrm_riv_temp_class::dealloc_lateral
    procedure :: dealloc_lateral !< \see mo_mrm_riv_temp_class::dealloc_lateral
    !> \copydoc mo_mrm_riv_temp_class::clean_up
    procedure :: clean_up !< \see mo_mrm_riv_temp_class::clean_up
 
  end type riv_temp_type
 
contains
 
 
  !> \brief clean up
  subroutine clean_up( &
    self &
  )
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
 
    if ( allocated(self%L1_runoff_E) ) deallocate(self%L1_runoff_E)
    if ( allocated(self%L1_acc_strd) ) deallocate(self%L1_acc_strd)
    if ( allocated(self%L1_acc_ssrd) ) deallocate(self%L1_acc_ssrd)
    if ( allocated(self%L1_acc_temp) ) deallocate(self%L1_acc_temp)
    if ( allocated(self%dir_riv_widths) ) deallocate(self%dir_riv_widths)
    if ( allocated(self%L11_riv_widths) ) deallocate(self%L11_riv_widths)
    if ( allocated(self%L11_riv_areas) ) deallocate(self%L11_riv_areas)
    if ( allocated(self%netNode_E_IN) ) deallocate(self%netNode_E_IN)
    if ( allocated(self%netNode_E_R) ) deallocate(self%netNode_E_R)
    if ( allocated(self%netNode_E_mod) ) deallocate(self%netNode_E_mod)
    if ( allocated(self%netNode_E_out) ) deallocate(self%netNode_E_out)
    if ( allocated(self%L11_srad_net) ) deallocate(self%L11_srad_net)
    if ( allocated(self%L11_lrad_in) ) deallocate(self%L11_lrad_in)
    if ( allocated(self%L11_air_temp) ) deallocate(self%L11_air_temp)
    if ( allocated(self%river_temp) ) deallocate(self%river_temp)
    ! meteo arrays
    if ( allocated(self%L1_ssrd_calc) ) deallocate(self%L1_ssrd_calc)
    if ( allocated(self%L1_strd_calc) ) deallocate(self%L1_strd_calc)
    if ( allocated(self%L1_tann_calc) ) deallocate(self%L1_tann_calc)
 
  end subroutine clean_up
 
  !> \brief configure the \ref riv_temp_type class from the mhm namelist
  subroutine config( &
    self, &
    file_namelist, &
    file_namelist_param &
  )
 
    use mo_common_constants, only : maxnodomains, nodata_i4
    use mo_common_variables, only : domainmeta
    use mo_check, only : check_dir
    USE mo_string_utils, ONLY : num2str
    use mo_namelists, only : nml_config_riv_temp
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    character(*), intent(in) :: file_namelist !< mhm namelist file
    character(*), intent(in) :: file_namelist_param !< mhm parameter namelist file
 
    integer(i4) :: iDomain, domainID
 
    ! allocate the directory arrays
    allocate(self%dir_riv_widths(domainmeta%nDomains))
 
    call nml_config_riv_temp%read(file_namelist)
    self%albedo_water = nml_config_riv_temp%albedo_water
    self%pt_a_water = nml_config_riv_temp%pt_a_water
    self%emissivity_water = nml_config_riv_temp%emissivity_water
    self%turb_heat_ex_coeff = nml_config_riv_temp%turb_heat_ex_coeff
    self%delta_iter = nml_config_riv_temp%delta_iter
    self%max_iter = nml_config_riv_temp%max_iter
    self%step_iter = nml_config_riv_temp%step_iter
    self%riv_widths_file = nml_config_riv_temp%riv_widths_file
    self%riv_widths_name = nml_config_riv_temp%riv_widths_name
    self%dir_riv_widths = nml_config_riv_temp%dir_riv_widths(1:domainmeta%nDomains)
 
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      call check_dir( &
        path=self%dir_riv_widths(idomain), &
        text="(domain "//trim(num2str(domainid,'(I3)'))//") River widths directory:", &
        raise=.true., &
        tab=4, &
        text_length=40 &
      )
    end do
 
  end subroutine config
 
  !> \brief initalize the \ref riv_temp_type class for the current domain
  subroutine init( &
    self, &
    nCells &
  )
    use mo_append, only : append
    use mo_mrm_constants, only : nroutingstates
    ! use mo_common_variables, only : level0, domainMeta
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    integer(i4), intent(in) :: nCells !< number of level-11 cells for the current domain
 
    real(dp), dimension(:), allocatable :: dummy_Vector11
    real(dp), dimension(:, :), allocatable :: dummy_Matrix11_IT
 
    ! dummy vector and matrix
    allocate(dummy_vector11(ncells))
    allocate(dummy_matrix11_it(ncells, nroutingstates))
    dummy_vector11(:) = 0.0_dp
    dummy_matrix11_it(:, :) = 0.0_dp
 
    ! simulated energy flux at each node
    call append(self%netNode_E_mod, dummy_vector11)
    ! simulated river temperature at each node
    call append(self%river_temp, dummy_vector11)
    ! Total outflow from cells L11 at time tt
    call append(self%netNode_E_out, dummy_vector11)
    ! Total discharge inputs at t-1 and t
    call append(self%netNode_E_IN, dummy_matrix11_it)
    !  Routed outflow leaving a node
    call append(self%netNode_E_R, dummy_matrix11_it)
 
    ! free space
    deallocate(dummy_vector11)
    deallocate(dummy_matrix11_it)
 
  end subroutine init
 
  !> \brief initialize the river area of \ref riv_temp_type class for the current domain
  subroutine init_area( &
    self, &
    iDomain, &
    L11_netPerm, &
    L11_fromN, &
    L11_length, &
    nLinks, &
    nCells, &
    nrows, &
    ncols, &
    L11_mask &
  )
    use mo_append, only : append
    use mo_read_nc, only: read_const_nc
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    integer(i4), intent(in) :: iDomain !< Domain ID
    integer(i4), dimension(:), intent(in) :: L11_netPerm !< L11 routing order
    integer(i4), dimension(:), intent(in) :: L11_fromN !< L11 source grid cell order
    real(dp), dimension(:), intent(in) :: L11_length !< L11 link length
    integer(i4), intent(in) :: nLinks !< number of L11 links in the current domain
    integer(i4), intent(in) :: nCells !< number of L11 cells of the current domain
    integer(i4), intent(in) :: ncols  !< Number of columns
    integer(i4), intent(in) :: nrows  !< Number of rows
    logical, dimension(:, :), intent(in) :: L11_mask !< the mask for valid cells in the original grid (nrows, ncols)
 
    real(dp), dimension(:,:), allocatable :: L11_data ! read data from file
    real(dp), dimension(:), allocatable :: L11_riv_widths, L11_riv_areas
 
    integer(i4) :: i, k, iNode
 
    call read_const_nc(&
      trim(self%dir_riv_widths(idomain)), &
      nrows, &
      ncols, &
      self%riv_widths_name, &
      l11_data, &
      self%riv_widths_file &
    )
 
    allocate(l11_riv_widths(ncells))
    l11_riv_widths(:) = pack(l11_data(:,:), mask=l11_mask)
    call append(self%L11_riv_widths, l11_riv_widths)
 
    allocate(l11_riv_areas(ncells))
    ! at area at Outlet is 0
    l11_riv_areas = 0.0_dp
    do k = 1, nlinks
      ! get LINK routing order -> i
      i = l11_netperm(k)
      inode = l11_fromn(i)
      l11_riv_areas(inode) = l11_length(i) * l11_riv_widths(inode)
    end do
    call append(self%L11_riv_areas, l11_riv_areas)
 
    deallocate(l11_data)
    deallocate(l11_riv_widths)
    deallocate(l11_riv_areas)
 
  end subroutine init_area
 
  !> \brief initialize the river temperature of \ref riv_temp_type class for the current domain
  subroutine init_riv_temp( &
    self, &
    temp_air, &
    efecarea, &
    L1_L11_Id, &
    L11_areacell, &
    L11_L1_Id, &
    map_flag &
  )
 
    use mo_mrm_pre_routing, only : l11_meteo_acc
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    real(dp), dimension(:), intent(in) :: temp_air !< air temperature [degC] for current timestep
    real(dp), intent(in), dimension(:) :: efecarea !< effective area in [km2] at Level 1
    integer(i4), intent(in), dimension(:) :: L1_L11_Id !< L11 Ids mapped on L1
    real(dp), intent(in), dimension(:) :: L11_areacell !< effective area in [km2] at Level 11
    integer(i4), intent(in), dimension(:) :: L11_L1_Id !< L1 Ids mapped on L11
    logical, intent(in) :: map_flag !< Flag indicating whether routing resolution is higher than hydrologic one
 
    ! map temperature from L1 to L11
    call l11_meteo_acc( &
      temp_air, efecarea, l1_l11_id, l11_areacell, l11_l1_id, map_flag, self%river_temp(self%s11 : self%e11))
 
    ! assure positive temperature
    self%river_temp(self%s11 : self%e11) = max(self%delta_T, self%river_temp(self%s11 : self%e11))
 
  end subroutine init_riv_temp
 
  !> \brief reset \ref riv_temp_type class for next timestep
  subroutine reset_timestep(self)
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
 
    ! set sub time-step counter to 0
    self%ts_cnt = 0_i4
 
    self%L1_runoff_E = 0.0_dp
    self%L1_acc_strd = 0.0_dp
    self%L1_acc_ssrd = 0.0_dp
    self%L1_acc_temp = 0.0_dp
 
  end subroutine reset_timestep
 
  !> \brief allocate lateral temp components of \ref riv_temp_type class for current domain
  subroutine alloc_lateral( &
    self, &
    nCells &
  )
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    integer(i4), intent(in) :: nCells !< number of level-1 cells for the current domain
 
    allocate(self%L1_runoff_E(ncells))
    allocate(self%L1_acc_strd(ncells))
    allocate(self%L1_acc_ssrd(ncells))
    allocate(self%L1_acc_temp(ncells))
    ! meteo arrays
    allocate(self%L1_ssrd_calc(ncells))
    allocate(self%L1_strd_calc(ncells))
    allocate(self%L1_tann_calc(ncells))
    ! init these arrays to 0
    call self%reset_timestep()
 
  end subroutine alloc_lateral
 
  !> \brief deallocate lateral temp components of \ref riv_temp_type
  subroutine dealloc_lateral( &
    self &
  )
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
 
    deallocate(self%L1_runoff_E)
    deallocate(self%L1_acc_strd)
    deallocate(self%L1_acc_ssrd)
    deallocate(self%L1_acc_temp)
    deallocate(self%L1_ssrd_calc)
    deallocate(self%L1_strd_calc)
    deallocate(self%L1_tann_calc)
 
  end subroutine dealloc_lateral
 
  !> \brief accumulate energy sources of \ref riv_temp_type
  subroutine acc_source_e( &
    self, &
    fSealed_area_fraction, &
    fast_interflow, &
    slow_interflow, &
    baseflow, &
    direct_runoff, &
    temp_air &
  )
 
    use mo_mrm_pre_routing, only : calc_l1_runoff_e
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    real(dp), dimension(:), intent(in) :: fSealed_area_fraction !< sealed area fraction [1]
    real(dp), dimension(:), intent(in) :: fast_interflow !< \f$ q_0 \f$ Fast runoff component [mm TS-1]
    real(dp), dimension(:), intent(in) :: slow_interflow !< \f$ q_1 \f$ Slow runoff component [mm TS-1]
    real(dp), dimension(:), intent(in) :: baseflow !< \f$ q_2 \f$ Baseflow [mm TS-1]
    real(dp), dimension(:), intent(in) :: direct_runoff !< \f$ q_D \f$ Direct runoff from impervious areas  [mm TS-1]
    real(dp), dimension(:), intent(in) :: temp_air !< air temperature [K]
 
    ! increase the sub time-step counter
    self%ts_cnt = self%ts_cnt + 1_i4
 
    ! caclucate the temperature energy of the runoffs at L1 in [K mm]
    ! automatically accumulate them
    call calc_l1_runoff_e( &
      fsealed_area_fraction, &
      fast_interflow, slow_interflow, baseflow, direct_runoff, &
      temp_air, self%L1_tann_calc, &
      self%L1_runoff_E & ! will be added here
    )
    ! accumulate meteo forcings (will be averaged with sub time-step counter later)
    self%L1_acc_ssrd = self%L1_acc_ssrd + self%L1_ssrd_calc
    self%L1_acc_strd = self%L1_acc_strd + self%L1_strd_calc
    self%L1_acc_temp = self%L1_acc_temp + temp_air
 
  end subroutine acc_source_e
 
  !> \brief finalize energy sources of \ref riv_temp_type
  subroutine finalize_source_e( &
    self, &
    efecarea, &
    L1_L11_Id, &
    L11_areacell, &
    L11_L1_Id, &
    timestep, &
    map_flag &
  )
 
    use mo_mrm_pre_routing, only : l11_meteo_acc, l11_runoff_acc
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    !> effective area in [km2] at Level 1
    real(dp), intent(in), dimension(:) :: efecarea
    !> L11 Ids mapped on L1
    integer(i4), intent(in), dimension(:) :: L1_L11_Id
    !> effective area in [km2] at Level 11
    real(dp), intent(in), dimension(:) :: L11_areacell
    !> L1 Ids mapped on L11
    integer(i4), intent(in), dimension(:) :: L11_L1_Id
    !> simulation timestep in [h]
    integer(i4), intent(in) :: timestep
    !> Flag indicating whether routing resolution is higher than hydrologic one
    logical, intent(in) :: map_flag
 
    ! prepare temporal variables for temp-routing
    if ( allocated(self%L11_srad_net) ) deallocate(self%L11_srad_net)
    if ( allocated(self%L11_lrad_in) ) deallocate(self%L11_lrad_in)
    if ( allocated(self%L11_air_temp) ) deallocate(self%L11_air_temp)
    allocate(self%L11_srad_net(size(l11_l1_id, dim = 1)))
    allocate(self%L11_lrad_in(size(l11_l1_id, dim = 1)))
    allocate(self%L11_air_temp(size(l11_l1_id, dim = 1)))
 
    ! convert the L1 runoff temperature energy from [K mm] to [K m3 s-1] on L11
    call l11_runoff_acc( &
      self%L1_runoff_E, &
      efecarea, &
      l1_l11_id, &
      l11_areacell, &
      l11_l1_id, &
      timestep, &
      map_flag, &
      self%netNode_E_out(self%s11 : self%e11) &
    )
    ! all meteo forcings are accumulated over sub time-steps and need to be divided by the sub time-step counter
    ! short wave radiation on L11
    self%L1_acc_ssrd = self%L1_acc_ssrd / real(self%ts_cnt, dp)
    call l11_meteo_acc(self%L1_acc_ssrd, efecarea, l1_l11_id, l11_areacell, l11_l1_id, map_flag, self%L11_srad_net)
    ! calculate the net-radiation from incoming/outging short radiation
    self%L11_srad_net = self%L11_srad_net * (1._dp - self%albedo_water)
    ! long wave radiation on L11
    self%L1_acc_strd = self%L1_acc_strd / real(self%ts_cnt, dp)
    call l11_meteo_acc(self%L1_acc_strd, efecarea, l1_l11_id, l11_areacell, l11_l1_id, map_flag, self%L11_lrad_in)
    ! air temperature on L11
    self%L1_acc_temp = self%L1_acc_temp / real(self%ts_cnt, dp)
    call l11_meteo_acc(self%L1_acc_temp, efecarea, l1_l11_id, l11_areacell, l11_l1_id, map_flag, self%L11_air_temp)
 
  end subroutine finalize_source_e
 
  !> \brief get outgoing longwave radiation of \ref riv_temp_type
  !> \return outgoing longwave radiation
  real(dp) function get_lrad_out(self, riv_temp) result(lrad_out)
 
    use mo_constants, only: sigma_dp
 
    implicit none
 
    class(riv_temp_type), intent(in) :: self
    !> river temperature in K
    real(dp), intent(in) :: riv_temp
 
    ! outgoing longwave radiation from Boltzmann equation
    lrad_out = self%emissivity_water * sigma_dp * riv_temp ** 4_i4
 
  end function get_lrad_out
 
  !> \brief latent heat flux of \ref riv_temp_type
  !> \return latent heat flux
  real(dp) function get_lat_heat(self, air_temp, netrad) result(lat_heat)
 
    use mo_constants, only : psychro_dp
    use mo_pet, only : slope_satpressure
 
    implicit none
 
    class(riv_temp_type), intent(in) :: self
    !> air temperature in deg C
    real(dp), intent(in) :: air_temp
    !> net radiation in W * m-2
    real(dp), intent(in) :: netrad
 
    ! save slope of saturation vapor pressure curve
    real(dp) :: delta
 
    delta = slope_satpressure(air_temp) ! slope of saturation vapor pressure curve
    ! latent heat for priestley taylor PET formula
    lat_heat = self%pt_a_water * delta / (psychro_dp + delta) * max(netrad, 0._dp)
 
  end function get_lat_heat
 
  !> \brief sensible heat flux of \ref riv_temp_type
  !> \return sensible heat flux
  real(dp) function get_sens_heat(self, air_temp, riv_temp) result(sens_heat)
 
    implicit none
 
    class(riv_temp_type), intent(in) :: self
    real(dp), intent(in) :: air_temp !< air temperature in [deg C]
    real(dp), intent(in) :: riv_temp !< river temperature in [deg C]
 
    ! sensible heat flux resulting for temp. diff.: river <-> air
    sens_heat = self%turb_heat_ex_coeff * (riv_temp - air_temp)
 
  end function get_sens_heat
 
  !> \brief get complete energy source of \ref riv_temp_type at given cell
  !> \return energy IO
  real(dp) function get_e_io(self, riv_temp, cell) result(e_io)
 
    use mo_constants, only : t0_dp, cp_w_dp
    use mo_mhm_constants, only : h2odens
 
    implicit none
 
    class(riv_temp_type), intent(in) :: self
    !> given river temperature in K to calculate heat fluxes
    real(dp), intent(in) :: riv_temp
    !> cell index in the current domain
    integer(i4), intent(in) :: cell
    !> net radiation calc from short and longwave radiation in/out
    real(dp) :: netrad, sens_heat, lat_heat
 
    ! net radiation
    netrad = self%L11_srad_net(cell) + self%L11_lrad_in(cell) - self%get_lrad_out(riv_temp)
    ! sensible heat flux
    sens_heat = self%get_sens_heat(self%L11_air_temp(cell), riv_temp - t0_dp)
    ! latent heat flux
    lat_heat = self%get_lat_heat(self%L11_air_temp(cell), netrad)
 
    ! convert energy flux [W m-2] to [K m s-1]
    ! needs to be multiplied with river-area to result in temp-energy flux [K m3 s-1]
    e_io = (netrad - sens_heat - lat_heat) / (h2odens * cp_w_dp)
 
  end function get_e_io
 
  !> \brief execute the temperature routing of \ref riv_temp_type
  subroutine l11_routing_e( &
    self, &
    nLinks, &
    netPerm, &
    netLink_fromN, &
    netLink_toN, &
    netLink_C1, &
    netLink_C2, &
    nInflowGauges, &
    InflowHeadwater, &
    InflowNodeList, &
    sink_cells, &
    L11_qTR, &
    L11_Qmod &
  )
    use mo_constants, only : t0_dp
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    !> number of stream segment (reaches)
    integer(i4), intent(in) :: nLinks
    !> routing order of a given domain (permutation)
    integer(i4), dimension(:), intent(in) :: netPerm
    !> from node
    integer(i4), dimension(:), intent(in) :: netLink_fromN
    !> to node
    integer(i4), dimension(:), intent(in) :: netLink_toN
    !> routing parameter  C1 (\cite CMM1988 p. 25-41)
    real(dp), dimension(:), intent(in) :: netLink_C1
    !> routing parameters C2 (id)
    real(dp), dimension(:), intent(in) :: netLink_C2
    !> [-]      number of inflow points
    integer(i4), intent(in) :: nInflowGauges
    !> [-]      if to consider headwater cells of inflow gauge
    logical, dimension(:), intent(in) :: InflowHeadwater
    !> [-]      L11 ID of inflow points
    integer(i4), dimension(:), intent(in) :: InflowNodeList
    ! [-]      sink nodes
    integer(i4), dimension(:), intent(in) :: sink_cells
    !> [m3 s-1] Transformed outflow leaving node I at current timestep(Muskingum)
    real(dp), intent(in), dimension(:) :: L11_qTR
    !> [m3 s-1] Simulated routed discharge
    real(dp), intent(in), dimension(:) :: L11_Qmod
 
    integer(i4) :: i, k, m, iNode, tNode, L11in, L11to
    real(dp) :: E_IO_in, E_IO, T_est, T_rout
 
    ! initialize total input at point time (t+1) in all nodes
    self%netNode_E_IN(self%s11 : self%e11, 2) = 0.0_dp
 
    riv_loop: do k = 1, nlinks
      ! get LINK routing order -> i
      i = netperm(k)
      inode = netlink_fromn(i)
      tnode = netlink_ton(i)
      ! indices for concaternated domains
      l11in = inode + self%s11 - 1_i4
      l11to = tnode + self%s11 - 1_i4
 
      ! accumulate all inputs in iNode
      self%netNode_E_IN(l11in, 2) = self%netNode_E_IN(l11in, 2) + self%netNode_E_out(l11in)
      ! calculate the river temp at the upstream node
      self%river_temp(l11in) = self%netNode_E_IN(l11in, 2) / l11_qmod(inode) - t0_dp
 
      ! first guess for the routed temperature is the upstream temp
      t_est = self%river_temp(l11in) + t0_dp
      ! calculate the meteo energy source at the upstream node
      e_io_in = self%get_E_IO(t_est, inode)
 
      ! initialize iteration
      call self%init_iter()
      ! perform iteration
      iterloop: do m=1, self%max_iter
        ! meteo energy flux depending on the estimated routed temperature (avg. IN- and OUT-node)
        e_io = (e_io_in + self%get_E_IO(t_est, tnode)) * 0.5_dp * self%L11_riv_areas(l11in)
        ! routing iNode
        self%netNode_E_R(l11in, 2) = self%netNode_E_R(l11in, 1) &
          + netlink_c1(i) * (self%netNode_E_IN(l11in, 1) - self%netNode_E_R(l11in, 1)) &
          + netlink_c2(i) * (self%netNode_E_IN(l11in, 2) + e_io - self%netNode_E_IN(l11in, 1))
        ! calculate the routed temperature
        t_rout = self%netNode_E_R(l11in, 2) / l11_qtr(inode)
        ! check for convergence
        if ( abs(t_rout - t_est) < self%delta_iter ) exit iterloop
        ! get new estimate for next iteration
        call self%next_iter(t_est, t_rout)
      end do iterloop
 
      ! add the meteo energy flux to the accumulated incoming energy at the IN-node
      self%netNode_E_IN(l11in, 2) = self%netNode_E_IN(l11in, 2) + e_io
 
      ! check if the inflow from upstream cells should be deactivated
      if (ninflowgauges .GT. 0) then
        do i = 1, ninflowgauges
          ! check if downstream Node (tNode) is inflow gauge and headwaters should be ignored
          if ((tnode == inflownodelist(i)) .AND. (.NOT. inflowheadwater(i))) &
            self%netNode_E_R(l11in, 2) = 0.0_dp
        end do
      end if
 
      ! add routed temp energy to downstream node
      self%netNode_E_IN(l11to, 2) = self%netNode_E_IN(l11to, 2) + self%netNode_E_R(l11in, 2)
 
    end do riv_loop
 
    ! Accumulate all inputs at sinks
    ! calculate riv temp at sinks
    do k = 1, size(sink_cells)
      tnode = sink_cells(k)
      l11to = tnode + self%s11 - 1_i4
      self%netNode_E_IN(l11to, 2) = self%netNode_E_IN(l11to, 2) + self%netNode_E_out(l11to)
      self%river_temp(l11to) = self%netNode_E_IN(l11to, 2) / l11_qmod(tnode) - t0_dp
    end do
 
    ! backflow t-> t-1
    self%netNode_E_R(self%s11 : self%e11, 1) = self%netNode_E_R(self%s11 : self%e11, 2)
    self%netNode_E_IN(self%s11 : self%e11, 1) = self%netNode_E_IN(self%s11 : self%e11, 2)
 
  end subroutine l11_routing_e
 
  !> \brief initialize iterative solver of \ref riv_temp_type
  subroutine init_iter(self)
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
 
    ! first iteration is used to determine the search direction for the interval
    self%first_iter = .true.
    ! after the interval is found, we will perform a biseciton to nail down the temperature
    self%bisect_iter = .false.
 
  end subroutine init_iter
 
  !> \brief execute next iteration with iterative solver of \ref riv_temp_type
  subroutine next_iter(self, T_est, T_rout)
 
    implicit none
 
    class(riv_temp_type), intent(inout) :: self
    real(dp), intent(inout) :: T_est !< estimated river temperature
    real(dp), intent(in) :: T_rout !< calculated (routed) river temperature
 
    ! before performing a bisection we need to search for the interval (with given step-size)
    if ( .not. self%bisect_iter ) then
      ! to determine the direction of interval search, we have to wait for the first iteration
      if ( self%first_iter ) then
        self%first_iter = .false.
        ! determine search direction
        self%up_iter = ( t_rout > t_est )
      end if
      ! check if we need to take another step
      if ( self%up_iter .eqv. ( t_rout > t_est ) ) then
        if ( self%up_iter ) then
          t_est = t_est + self%step_iter
        else
          t_est = t_est - self%step_iter
        end if
      ! if direction changes, we have found the interval for bisection
      else
        ! start bisection in next iteration
        self%bisect_iter = .true.
        ! set interval bounds for bisection
        if ( self%up_iter ) then
          self%up_bnd_iter = t_est
          self%low_bnd_iter = t_est - self%step_iter
        else
          self%up_bnd_iter = t_est + self%step_iter
          self%low_bnd_iter = t_est
        end if
        ! set estimation to interval center
        t_est = (self%up_bnd_iter + self%low_bnd_iter) / 2.0_dp
      end if
    ! perform bisection
    else
      ! if routed temp is lower than estimated, use lower interval
      if ( t_rout < t_est ) then
        self%up_bnd_iter = t_est
      ! if routed temp is higher than estimated, use upper interval
      else
        self%low_bnd_iter = t_est
      end if
      ! set estimation to interval center
      t_est = (self%up_bnd_iter + self%low_bnd_iter) / 2.0_dp
    end if
 
  end subroutine next_iter
 
end module mo_mrm_riv_temp_class