mo_mrm_read_config.f90

Go to the documentation of this file.

!> \file mo_mrm_read_config.f90
!> \brief \copybrief mo_mrm_read_config
!> \details \copydetails mo_mrm_read_config
 
!> \brief read mRM config
!> \details This module contains all mRM subroutines related to reading the mRM configuration either from file or copy from mHM.
!> \authors 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_read_config
 
  use mo_kind, only : i4, dp
  use mo_message, only: message, error_message
 
  implicit none
 
  public :: mrm_read_config
 
contains
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_read_config
 
  !    PURPOSE
  !>       \brief Read the general config of mRM
 
  !>       \details Depending on the variable mrm_coupling_config, the
  !>       mRM config is either read from mrm.nml and parameters from
  !>       mrm_parameter.nml or copied from mHM.
 
  !    INTENT(IN)
  !>       \param[in] "character(*) :: file_namelist, file_namelist_param"
  !>       \param[in] "character(*) :: file_namelist, file_namelist_param"
  !>       \param[in] "logical :: do_message"                              - flag for writing mHM standard messages
 
  !    HISTORY
  !>       \authors Stephan Thober
 
  !>       \date Aug 2015
 
  ! Modifications:
  ! Stephan Thober  Sep 2015 - removed stop condition when routing resolution is smaller than hydrologic resolution
  ! Stephan Thober  Oct 2015 - added NLoutputResults namelist, fileLatLon to directories_general namelist, and readLatLon flag
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine mrm_read_config(file_namelist, file_namelist_param, do_message)
 
    use mo_common_constants, only : maxnodomains, nodata_i4
    use mo_common_mhm_mrm_read_config, only : common_check_resolution
    use mo_common_variables, only : alma_convention, domainmeta, processmatrix
    use mo_mrm_constants, only : maxnogauges
    use mo_mrm_file, only : file_defoutput
    use mo_mrm_global_variables, only : inflowgauge, domaininfo_mrm, domain_mrm, &
                                        dirgauges, dirtotalrunoff, filenametotalrunoff, dirbankfullrunoff, gauge, is_start, &
                                        ngaugestotal, ngaugeslocal, ninflowgaugestotal, outputflxstate_mrm, &
                                        timestep_model_outputs_mrm, &
                                        varnametotalrunoff, gw_coupling, &
                                        output_deflate_level_mrm, output_double_precision_mrm, output_time_reference_mrm, &
                                        readlatlon
    use mo_string_utils, only : num2str
    use mo_namelists, only : &
      nml_mainconfig_mrm, &
      nml_directories_mrm, &
      nml_evaluation_gauges, &
      nml_inflow_gauges, &
      nml_mrm_outputs
 
    implicit none
 
    character(*), intent(in) :: file_namelist, file_namelist_param
 
    ! - flag for writing mHM standard messages
    logical, intent(in) :: do_message
 
    integer(i4), dimension(maxNoDomains) :: nogauges_domain
    integer(i4), dimension(maxNoDomains, maxNoGauges) :: gauge_id
    character(256), dimension(maxNoDomains, maxNoGauges) :: gauge_filename
 
    integer(i4), dimension(maxNoDomains) :: noinflowgauges_domain
    integer(i4), dimension(maxNoDomains, maxNoGauges) :: inflowgauge_id
    character(256), dimension(maxNoDomains, maxNoGauges) :: inflowgauge_filename
    logical, dimension(maxNoDomains, maxNoGauges) :: inflowgauge_headwater
 
    integer(i4) :: domainid, idomain
 
    integer(i4) :: igauge
 
    integer(i4) :: idx
 
    character(256), dimension(maxNoDomains) :: dir_gauges
    character(256), dimension(maxNoDomains) :: dir_total_runoff
    character(256), dimension(maxNoDomains) :: dir_bankfull_runoff
 
    logical :: file_exists
 
    type(domaininfo_mrm), pointer :: domain_mrm_idomain
 
    !===============================================================
    ! INITIALIZATION
    !===============================================================
    is_start = .true.
 
    !===============================================================
    !  Read namelist for mainconfig for mRM
    !===============================================================
    call nml_mainconfig_mrm%read(file_namelist)
    alma_convention = nml_mainconfig_mrm%ALMA_convention
    filenametotalrunoff = nml_mainconfig_mrm%filenameTotalRunoff
    varnametotalrunoff = nml_mainconfig_mrm%varnameTotalRunoff
    gw_coupling = nml_mainconfig_mrm%gw_coupling
 
    !===============================================================
    !  Read namelist for mainpaths
    !===============================================================
    call nml_directories_mrm%read(file_namelist)
    dir_gauges = nml_directories_mrm%dir_Gauges
    dir_total_runoff = nml_directories_mrm%dir_Total_Runoff
    dir_bankfull_runoff = nml_directories_mrm%dir_Bankfull_Runoff
 
    allocate(dirgauges(domainmeta%nDomains), dirtotalrunoff(domainmeta%nDomains), dirbankfullrunoff(domainmeta%nDomains))
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      dirgauges(idomain)         = dir_gauges(domainid)
      dirtotalrunoff(idomain)    = dir_total_runoff(domainid)
      dirbankfullrunoff(idomain) = dir_bankfull_runoff(domainid)
    end do
 
    !===============================================================
    ! READ EVALUATION GAUGES
    !===============================================================
    call nml_evaluation_gauges%read(file_namelist)
    ngaugestotal = nml_evaluation_gauges%nGaugesTotal
    nogauges_domain = nml_evaluation_gauges%NoGauges_domain
    gauge_id = nml_evaluation_gauges%Gauge_id
    gauge_filename = nml_evaluation_gauges%gauge_filename
 
    if (ngaugestotal .GT. maxnogauges) then
      call error_message('***ERROR: ', trim(file_namelist), ': Total number of evaluation gauges is restricted to', &
              num2str(maxnogauges), raise=.false.)
      call error_message('          Error occured in namlist: evaluation_gauges')
    end if
 
    ! ToDo: check
    ngaugeslocal = 0
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      ngaugeslocal = ngaugeslocal + nogauges_domain(domainid)
    end do
    ! End ToDo
 
    allocate(gauge%gaugeId(ngaugeslocal)) ; gauge%gaugeId = nodata_i4
    allocate(gauge%domainId(ngaugeslocal)) ; gauge%domainId = nodata_i4
    allocate(gauge%fName  (ngaugeslocal))
    if (ngaugeslocal > 0) then
      gauge%fName(1) = num2str(nodata_i4)
    end if
    allocate(domain_mrm(domainmeta%nDomains))
 
    idx = 0
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      domain_mrm_idomain => domain_mrm(idomain)
      ! initialize
      domain_mrm_idomain%nGauges = nodata_i4
      allocate(domain_mrm_idomain%gaugeIdList(maxval(nogauges_domain(:))))
      domain_mrm_idomain%gaugeIdList = nodata_i4
      allocate(domain_mrm_idomain%gaugeIndexList(maxval(nogauges_domain(:))))
      domain_mrm_idomain%gaugeIndexList = nodata_i4
      allocate(domain_mrm_idomain%gaugeNodeList(maxval(nogauges_domain(:))))
      domain_mrm_idomain%gaugeNodeList = nodata_i4
      ! check if NoGauges_domain has a valid value
      if (nogauges_domain(domainid) .EQ. nodata_i4) then
        call error_message('***ERROR: ', trim(file_namelist), ': Number of evaluation gauges for subdomain ', &
                trim(adjustl(num2str(domainid))), ' is not defined!', raise=.false.)
        call error_message('          Error occured in namelist: evaluation_gauges')
      end if
 
      domain_mrm_idomain%nGauges = nogauges_domain(domainid)
 
      do igauge = 1, nogauges_domain(domainid)
        ! check if NoGauges_domain has a valid value
        if (gauge_id(domainid, igauge) .EQ. nodata_i4) then
          call error_message('***ERROR: ', trim(file_namelist), ': ID ', &
                  trim(adjustl(num2str(gauge_id(domainid, igauge)))), ' of evaluation gauge ', &
                  trim(adjustl(num2str(igauge))), ' for subdomain ', &
                  trim(adjustl(num2str(idomain))), ' is not defined!', raise=.false.)
          call error_message('          Error occured in namelist: evaluation_gauges')
        else if (trim(gauge_filename(domainid, igauge)) .EQ. trim(num2str(nodata_i4))) then
          call error_message('***ERROR: ', trim(file_namelist), ': Filename of evaluation gauge ', &
                  trim(adjustl(num2str(igauge))), ' for subdomain ', &
                  trim(adjustl(num2str(idomain))), ' is not defined!', raise=.false.)
          call error_message('          Error occured in namelist: evaluation_gauges')
        end if
        !
        idx = idx + 1
        gauge%domainId(idx) = idomain
        gauge%gaugeId(idx) = gauge_id(domainid, igauge)
        gauge%fname(idx) = trim(dirgauges(idomain)) // trim(gauge_filename(domainid, igauge))
        domain_mrm_idomain%gaugeIdList(igauge) = gauge_id(domainid, igauge)
        domain_mrm_idomain%gaugeIndexList(igauge) = idx
      end do
    end do
 
    if (ngaugeslocal .NE. idx) then
      call error_message('***ERROR: ', trim(file_namelist), ': Total number of evaluation gauges (', &
              trim(adjustl(num2str(ngaugeslocal))), &
              ') different from sum of gauges in subdomains (', trim(adjustl(num2str(idx))), ')!', raise=.false.)
      call error_message('          Error occured in namelist: evaluation_gauges')
    end if
 
    !===============================================================
    ! Read inflow gauge information
    !===============================================================
 
    call nml_inflow_gauges%read(file_namelist)
    ninflowgaugestotal = nml_inflow_gauges%nInflowGaugesTotal
    noinflowgauges_domain = nml_inflow_gauges%NoInflowGauges_domain
    inflowgauge_id = nml_inflow_gauges%InflowGauge_id
    inflowgauge_filename = nml_inflow_gauges%InflowGauge_filename
    inflowgauge_headwater = nml_inflow_gauges%InflowGauge_Headwater
 
    if (ninflowgaugestotal .GT. maxnogauges) then
      call error_message('***ERROR: ', trim(file_namelist), &
              ':read_gauge_lut: Total number of inflow gauges is restricted to', num2str(maxnogauges), raise=.false.)
      call error_message('          Error occured in namlist: inflow_gauges')
    end if
 
    ! allocation - max() to avoid allocation with zero, needed for mhm call
    allocate(inflowgauge%gaugeId (max(1, ninflowgaugestotal)))
    allocate(inflowgauge%domainId (max(1, ninflowgaugestotal)))
    allocate(inflowgauge%fName   (max(1, ninflowgaugestotal)))
    ! dummy initialization
    inflowgauge%gaugeId = nodata_i4
    inflowgauge%domainId = nodata_i4
    inflowgauge%fName = num2str(nodata_i4)
 
    idx = 0
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      domain_mrm_idomain => domain_mrm(idomain)
 
      allocate(domain_mrm_idomain%InflowGaugeIdList    (max(1, maxval(noinflowgauges_domain(:)))))
      allocate(domain_mrm_idomain%InflowGaugeHeadwater (max(1, maxval(noinflowgauges_domain(:)))))
      allocate(domain_mrm_idomain%InflowGaugeIndexList (max(1, maxval(noinflowgauges_domain(:)))))
      allocate(domain_mrm_idomain%InflowGaugeNodeList  (max(1, maxval(noinflowgauges_domain(:)))))
      ! dummy initialization
      domain_mrm_idomain%nInflowGauges = 0
      domain_mrm_idomain%InflowGaugeIdList = nodata_i4
      domain_mrm_idomain%InflowGaugeHeadwater = .false.
      domain_mrm_idomain%InflowGaugeIndexList = nodata_i4
      domain_mrm_idomain%InflowGaugeNodeList = nodata_i4
      ! no inflow gauge for subdomain i
      if (noinflowgauges_domain(domainid) .EQ. nodata_i4) then
        noinflowgauges_domain(domainid) = 0
      end if
 
      domain_mrm_idomain%nInflowGauges = noinflowgauges_domain(domainid)
 
      do igauge = 1, noinflowgauges_domain(domainid)
        ! check if NoInflowGauges_domain has a valid value
        if (inflowgauge_id(domainid, igauge) .EQ. nodata_i4) then
          call error_message('***ERROR: ', trim(file_namelist), ':ID of inflow gauge ', &
                  trim(adjustl(num2str(igauge))), ' for subdomain ', &
                  trim(adjustl(num2str(idomain))), ' is not defined!', raise=.false.)
          call error_message('          Error occured in namlist: inflow_gauges')
        else if (trim(inflowgauge_filename(domainid, igauge)) .EQ. trim(num2str(nodata_i4))) then
          call error_message('***ERROR: ', trim(file_namelist), ':Filename of inflow gauge ', &
                  trim(adjustl(num2str(igauge))), ' for subdomain ', &
                  trim(adjustl(num2str(idomain))), ' is not defined!', raise=.false.)
          call error_message('          Error occured in namlist: inflow_gauges')
        end if
        !
        idx = idx + 1
        inflowgauge%domainId(idx) = idomain
        inflowgauge%gaugeId(idx) = inflowgauge_id(domainid, igauge)
        inflowgauge%fname(idx) = trim(dirgauges(domainid)) // trim(inflowgauge_filename(domainid, igauge))
        domain_mrm_idomain%InflowGaugeIdList(igauge) = inflowgauge_id(domainid, igauge)
        domain_mrm_idomain%InflowGaugeHeadwater(igauge) = inflowgauge_headwater(domainid, igauge)
        domain_mrm_idomain%InflowGaugeIndexList(igauge) = idx
      end do
    end do
 
    if (ninflowgaugestotal .NE. idx) then
      call error_message('***ERROR: ', trim(file_namelist), ': Total number of inflow gauges (', &
              trim(adjustl(num2str(ninflowgaugestotal))), &
              ') different from sum of inflow gauges in subdomains (', trim(adjustl(num2str(idx))), ')!', raise=.false.)
      call error_message('          Error occured in namlist: inflow_gauges')
    end if
 
    call common_check_resolution(do_message, .true.)
 
    !===============================================================
    ! Read namelist global parameters
    !===============================================================
    call read_mrm_routing_params(processmatrix(8, 1), file_namelist_param)
 
    !===============================================================
    ! Read Output specifications for mRM
    !===============================================================
    call nml_mrm_outputs%read(file_defoutput)
    output_deflate_level_mrm = nml_mrm_outputs%output_deflate_level_mrm
    output_double_precision_mrm = nml_mrm_outputs%output_double_precision_mrm
    output_time_reference_mrm = nml_mrm_outputs%output_time_reference_mrm
    timestep_model_outputs_mrm = nml_mrm_outputs%timeStep_model_outputs_mrm
    outputflxstate_mrm = nml_mrm_outputs%outputFlxState_mrm
 
    readlatlon = any(outputflxstate_mrm)
 
    if (any(outputflxstate_mrm)) then
      call message('')
      call message('    Following output will be written:')
      call message('    NetCDF deflate level: ', adjustl(trim(num2str(output_deflate_level_mrm))))
      if ( output_double_precision_mrm ) then
        call message('    NetCDF output precision: double')
      else
        call message('    NetCDF output precision: single')
      end if
      select case(output_time_reference_mrm)
        case(0)
          call message('    NetCDF output time reference point: start of time interval')
        case(1)
          call message('    NetCDF output time reference point: center of time interval')
        case(2)
          call message('    NetCDF output time reference point: end of time interval')
      end select
      call message('    FLUXES:')
      if (outputflxstate_mrm(1)) then
        call message('      routed streamflow      (L11_qMod)                [m3 s-1]')
      end if
      if (outputflxstate_mrm(2)) then
        call message('      river temperature      (RivTemp)                 [deg C]')
      end if
      if (gw_coupling) then
        call message('      river head             (river_head)              [m]')
      end if
    end if
 
    call message('')
    call message('    FINISHED reading config')
    call message('')
 
  end subroutine mrm_read_config
 
  ! ---------------------------------------------------------------------------
  ! SUBROUTINE READ_MRM_ROUTING_PARAMS
  ! ---------------------------------------------------------------------------
  !    NAME
  !        read_mrm_routing_params
 
  !    PURPOSE
  !>       \brief TODO: add description
 
  !>       \details TODO: add description
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4) :: processCase"          it is the default case, should be one
  !>       \param[in] "character(*) :: file_namelist_param" file name containing parameter namelist
 
  !    HISTORY
  !>       \authors Robert Schweppe
 
  !>       \date Jun 2018
 
  ! Modifications:
 
  subroutine read_mrm_routing_params(processCase, file_namelist_param)
 
    use mo_common_constants, only : ncolpars
    use mo_common_functions, only : in_bound
    use mo_common_variables, only : global_parameters, global_parameters_name, processmatrix
    use mo_message, only : message
    use mo_namelists, only : &
      nml_routing1, &
      nml_routing2, &
      nml_routing3
 
    implicit none
 
    ! it is the default case, should be one
    integer(i4), intent(in) :: processCase
 
    ! file name containing parameter namelist
    character(*), intent(in) :: file_namelist_param
 
    ! equals sum of previous parameters
    integer(i4) :: start_index
 
    real(dp), dimension(nColPars) :: muskingumTravelTime_constant
    real(dp), dimension(nColPars) :: muskingumTravelTime_riverLength
    real(dp), dimension(nColPars) :: muskingumTravelTime_riverSlope
    real(dp), dimension(nColPars) :: muskingumTravelTime_impervious
    real(dp), dimension(nColPars) :: muskingumAttenuation_riverSlope
 
    real(dp), dimension(nColPars) :: streamflow_celerity
    real(dp), dimension(nColPars) :: slope_factor
 
    if (processcase .eq. 1_i4) then
      call nml_routing1%read(file_namelist_param)
      muskingumtraveltime_constant = nml_routing1%muskingumTravelTime_constant
      muskingumtraveltime_riverlength = nml_routing1%muskingumTravelTime_riverLength
      muskingumtraveltime_riverslope = nml_routing1%muskingumTravelTime_riverSlope
      muskingumtraveltime_impervious = nml_routing1%muskingumTravelTime_impervious
      muskingumattenuation_riverslope = nml_routing1%muskingumAttenuation_riverSlope
    else if (processcase .eq. 2_i4) then
      call nml_routing2%read(file_namelist_param)
      streamflow_celerity = nml_routing2%streamflow_celerity
    else if (processcase .eq. 3_i4) then
      call nml_routing3%read(file_namelist_param)
      slope_factor = nml_routing3%slope_factor
    end if
 
    ! -------------------------------------------------------------------------
    ! INCLUDE MRM PARAMETERS IN PARAMETERS OF MHM
    ! -------------------------------------------------------------------------
    ! Muskingum routing parameters with MPR
    if (processcase .eq. 1_i4) then
      ! insert parameter values and names at position required by mhm
      processmatrix(8, 1) = processcase
      processmatrix(8, 2) = 5_i4
      processmatrix(8, 3) = sum(processmatrix(1 : 8, 2))
      start_index = processmatrix(8, 3) - processmatrix(8, 2)
      global_parameters(start_index + 1, :) = muskingumtraveltime_constant
      global_parameters(start_index + 2, :) = muskingumtraveltime_riverlength
      global_parameters(start_index + 3, :) = muskingumtraveltime_riverslope
      global_parameters(start_index + 4, :) = muskingumtraveltime_impervious
      global_parameters(start_index + 5, :) = muskingumattenuation_riverslope
 
      global_parameters_name(start_index + 1 : start_index + processmatrix(8, 2)) = (/ &
              'muskingumTravelTime_constant   ', &
                      'muskingumTravelTime_riverLength', &
                      'muskingumTravelTime_riverSlope ', &
                      'muskingumTravelTime_impervious ', &
                      'muskingumAttenuation_riverSlope'/)
      ! adaptive timestep routing
    else if (processcase .eq. 2_i4) then
      processmatrix(8, 1) = processcase
      processmatrix(8, 2) = 1_i4
      processmatrix(8, 3) = sum(processmatrix(1 : 8, 2))
      start_index = processmatrix(8, 3) - processmatrix(8, 2)
      global_parameters(start_index + 1, :) = streamflow_celerity
 
      global_parameters_name(start_index + 1 : start_index + processmatrix(8, 2)) = (/ &
              'streamflow_celerity'/)
     ! adaptive timestep routing - varying celerity
     else if (processcase .eq. 3_i4) then
        ! insert parameter values and names at position required by mhm
        processmatrix(8, 1) = processcase
        processmatrix(8, 2) = 1_i4
        processmatrix(8, 3) = sum(processmatrix(1:8, 2))
        start_index         = processmatrix(8, 3) - processmatrix(8, 2)
        global_parameters(start_index + 1, :) = slope_factor
 
        global_parameters_name(start_index + 1 : start_index + processmatrix(8,2)) = (/ &
             'slope_factor'/)
    end if
 
    ! check if parameter are in range
    if (.not. in_bound(global_parameters)) then
      call error_message('***ERROR: parameter in routing namelist out of bound in ', &
              trim(adjustl(file_namelist_param)))
    end if
 
  end subroutine read_mrm_routing_params
end module mo_mrm_read_config