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] "integer :: unamelist, unamelist_param"
  !>       \param[in] "character(*) :: file_namelist, file_namelist_param"
  !>       \param[in] "integer :: unamelist, unamelist_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, unamelist, file_namelist_param, unamelist_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, udefoutput
    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_nml, only : close_nml, open_nml, position_nml
    use mo_string_utils, only : num2str
 
    implicit none
 
    character(*), intent(in) :: file_namelist, file_namelist_param
 
    integer, intent(in) :: unamelist, unamelist_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
 
 
    ! namelist spatial & temporal resolution, optmization information
    namelist /mainconfig_mrm/ alma_convention, &
      filenametotalrunoff, varnametotalrunoff, gw_coupling
    ! namelist directories
    namelist /directories_mrm/ dir_gauges, dir_total_runoff, dir_bankfull_runoff
    namelist /evaluation_gauges/ ngaugestotal, nogauges_domain, gauge_id, gauge_filename
    ! namelist for inflow gauges
    namelist /inflow_gauges/ ninflowgaugestotal, noinflowgauges_domain, inflowgauge_id, &
            inflowgauge_filename, inflowgauge_headwater
    ! name list regarding output
    namelist /nloutputresults/ &
            output_deflate_level_mrm, &
            output_double_precision_mrm, &
            output_time_reference_mrm, &
            timestep_model_outputs_mrm, &
            outputflxstate_mrm
 
    !===============================================================
    ! INITIALIZATION
    !===============================================================
    is_start = .true.
    ngaugestotal = nodata_i4
    nogauges_domain = nodata_i4
    gauge_id = nodata_i4
    gauge_filename = num2str(nodata_i4)
 
    ! default arguments
    alma_convention = .false.
    filenametotalrunoff = 'total_runoff'
    varnametotalrunoff = 'total_runoff'
    gw_coupling = .false.
 
    !===============================================================
    !  Read namelist main directories
    !===============================================================
    call open_nml(file_namelist, unamelist, quiet = .true.)
 
    !===============================================================
    !  Read namelist for mainconfig for mRM
    !===============================================================
    call position_nml('mainconfig_mrm', unamelist)
    read(unamelist, nml = mainconfig_mrm)
 
    !===============================================================
    !  Read namelist for mainpaths
    !===============================================================
    call position_nml('directories_mRM', unamelist)
    read(unamelist, nml = directories_mrm)
 
    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 position_nml('evaluation_gauges', unamelist)
    read(unamelist, nml = evaluation_gauges)
 
    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
    !===============================================================
 
    ninflowgaugestotal = 0
    noinflowgauges_domain = 0
    inflowgauge_id = nodata_i4
    inflowgauge_filename = num2str(nodata_i4)
 
    call position_nml('inflow_gauges', unamelist)
    read(unamelist, nml = inflow_gauges)
 
    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.)
 
    call close_nml(unamelist)
 
    !===============================================================
    ! Read namelist global parameters
    !===============================================================
    call read_mrm_routing_params(processmatrix(8, 1), file_namelist_param, unamelist_param)
 
    !===============================================================
    ! Read Output specifications for mRM
    !===============================================================
    output_deflate_level_mrm = 6
    output_time_reference_mrm = 0
    output_double_precision_mrm = .true.
    outputflxstate_mrm = .false.
    timestep_model_outputs_mrm = -2
    inquire(file = file_defoutput, exist = file_exists)
    if (file_exists) then
      ! file exists
      call open_nml(file_defoutput, udefoutput, quiet = .true.)
      call position_nml('NLoutputResults', udefoutput)
      read(udefoutput, nml = nloutputresults)
      call close_nml(udefoutput)
    else
      call message('')
      call message('No file specifying mRM output fluxes exists')
    end if
    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
  !>       \param[in] "integer(i4) :: unamelist_param"      file name id containing parameter namelist
 
  !    HISTORY
  !>       \authors Robert Schweppe
 
  !>       \date Jun 2018
 
  ! Modifications:
 
  subroutine read_mrm_routing_params(processCase, file_namelist_param, unamelist_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_nml, only : close_nml, open_nml, position_nml
 
    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
 
    ! file name id containing parameter namelist
    integer(i4), intent(in) :: unamelist_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
 
    namelist /routing1/ muskingumtraveltime_constant, muskingumtraveltime_riverlength, &
            muskingumtraveltime_riverslope, muskingumtraveltime_impervious, muskingumattenuation_riverslope
    namelist /routing2/ streamflow_celerity
    namelist /routing3/ slope_factor
    !
    call open_nml(file_namelist_param, unamelist_param, quiet = .true.)
 
    if (processcase .eq. 1_i4) then
      call position_nml('routing1', unamelist_param)
      read(unamelist_param, nml = routing1)
    else if (processcase .eq. 2_i4) then
       call position_nml('routing2', unamelist_param)
       read(unamelist_param, nml = routing2)
    else if (processcase .eq. 3_i4) then
       call position_nml('routing3', unamelist_param)
       read(unamelist_param, nml = routing3)
    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
 
    call close_nml(unamelist_param)
 
  end subroutine read_mrm_routing_params
end module mo_mrm_read_config