mo_common_mHM_mRM_read_config.f90

Go to the documentation of this file.

!> \file mo_common_mHM_mRM_read_config.f90
!> \brief \copybrief mo_common_mhm_mrm_read_config
!> \details \copydetails mo_common_mhm_mrm_read_config
 
!> \brief Reading of main model configurations.
!> \details This routine reads the configurations of common program parts
!> \authors Matthias Zink
!> \date Dec 2012
!> \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_common
MODULE mo_common_mhm_mrm_read_config
 
  use mo_kind, only : i4, dp
  use mo_message, only : message, error_message
 
  IMPLICIT NONE
 
  PRIVATE
 
  PUBLIC :: common_mhm_mrm_read_config, common_check_resolution, check_optimization_settings
 
CONTAINS
 
 
  !> \brief Read main configurations for common parts
  !> \changelog
  !! - Robert Schweppe Dec  2017
  !!   - based on mhm_read_config
  !! - Stephan Thober Jan 2022
  !!   - added nTStepForcingDay
  !> \authors Matthias Zink
  !> \date Dec 2012
  subroutine common_mhm_mrm_read_config(file_namelist, unamelist)
 
    use mo_common_constants, only : maxnodomains, nodata_i4
    use mo_common_mhm_mrm_variables, only : lcyearid, dds_r, mhmfilerestartin, mrmfilerestartin, evalper,&
                                            mcmc_error_params, mcmc_opti, niterations, &
                                            ntstepday, opti_function, opti_method, optimize, optimize_restart, &
                                            read_restart, mrm_read_river_network, resolutionrouting, sa_temp, &
                                            sce_ngs, sce_npg, sce_nps, seed, &
                                            simper, timestep, warmper, warmingdays, read_old_style_restart_bounds, &
                                            restart_reset_fluxes_states
    use mo_common_read_config, only : set_land_cover_scenes_id
    use mo_common_types, only: period
    use mo_common_variables, only : lcfilename, domainmeta, processmatrix
    use mo_julian, only : caldat, julday
    use mo_nml, only : close_nml, open_nml, position_nml
    use mo_string_utils, only : num2str
 
    implicit none
 
    character(*), intent(in) :: file_namelist !< namelist file name
    integer, intent(in) :: unamelist !< unit to open namelist file
 
    integer(i4) :: jday
 
    integer(i4) :: domainid, idomain
 
    integer(i4), dimension(maxNoDomains) :: warming_days
 
    type(period), dimension(maxNoDomains) :: eval_per
 
    real(dp), dimension(maxNoDomains) :: resolution_routing
 
    character(256), dimension(maxNoDomains) :: mhm_file_restartin
    character(256), dimension(maxNoDomains) :: mrm_file_restartin
 
 
    ! namelist spatial & temporal resolution, otmization information
    namelist /mainconfig_mhm_mrm/ timestep, resolution_routing, optimize, &
            optimize_restart, opti_method, opti_function, &
            read_restart, mrm_read_river_network, read_old_style_restart_bounds, restart_reset_fluxes_states, &
            mhm_file_restartin, mrm_file_restartin
    ! namelist for optimization settings
    namelist /optimization/ niterations, seed, dds_r, sa_temp, sce_ngs, &
            sce_npg, sce_nps, mcmc_opti, mcmc_error_params
    ! namelist for time settings
    namelist /time_periods/ warming_days, eval_per
 
    ! set default values for optional arguments
    mrm_read_river_network = .false.
    read_old_style_restart_bounds = .false.
    restart_reset_fluxes_states = .false.
 
    !===============================================================
    !  Read namelist main directories
    !===============================================================
    call open_nml(file_namelist, unamelist, quiet = .true.)
 
    !===============================================================
    !  Read namelist specifying the model configuration
    !===============================================================
    call position_nml('mainconfig_mhm_mrm', unamelist)
    read(unamelist, nml = mainconfig_mhm_mrm)
    ! consistency between read_restart and mrm_read_river_network
    if (read_restart) then
       if (.not. mrm_read_river_network) then
          call message('***WARNING: mrm_read_river_network is set to .true. because read_restart is .true.')
       end if
       mrm_read_river_network = .true.
    end if
 
    allocate(resolutionrouting(domainmeta%nDomains))
    allocate(mhmfilerestartin(domainmeta%nDomains))
    allocate(mrmfilerestartin(domainmeta%nDomains))
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      mhmfilerestartin(idomain) = mhm_file_restartin(domainid)
      mrmfilerestartin(idomain) = mrm_file_restartin(domainid)
      resolutionrouting(idomain) = resolution_routing(domainid)
    end do
 
    ! check for optimize and read restart
    if ((read_restart) .and. (optimize)) then
      call message()
      call error_message('***ERROR: cannot read states from restart file when optimizing')
    end if
 
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      if (processmatrix(8, 1) > 0 .and. domainmeta%optidata(idomain) > 1 .and. optimize) then
        domainmeta%doRouting(idomain) = .false.
        call message('Warning: although defined in namelist, routing is switched off for domain', trim(num2str(domainid)))
        call message('         since the calibration of Q is not possible with the chosen opti input')
      end if
    end do
 
    !===============================================================
    !  INIT !!! (merged from mo_startup and mo_mrm_read_config)
    !===============================================================
    ! transformation of time units & constants
    if (mod(24, timestep) > 0) then
      call error_message('mo_startup: timeStep must be a divisor of 24: ', num2str(timestep))
    end if
    ntstepday = 24_i4 / timestep            ! # of time steps per day
 
    ! allocate time periods
    allocate(simper(domainmeta%nDomains))
    allocate(evalper(domainmeta%nDomains))
    allocate(warmingdays(domainmeta%nDomains))
    allocate(warmper(domainmeta%nDomains))
 
    !===============================================================
    !  read simulation time periods incl. warming days
    !===============================================================
    call position_nml('time_periods', unamelist)
    read(unamelist, nml = time_periods)
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      warmingdays(idomain) = warming_days(domainid)
      ! this will be a procedure subroutine
      ! therefore inout first, in second
      call period_copy_period_data(evalper(idomain), eval_per(domainid))
    end do
    ! evalPer = eval_Per(1 : domainMeta%nDomains)
 
    !===============================================================
    !  determine simulation time period incl. warming days for each
    !  domain
    !===============================================================
    do idomain = 1, domainmeta%nDomains
      ! julian days for evaluation period
      jday = julday(dd = evalper(idomain)%dStart, mm = evalper(idomain)%mStart, yy = evalper(idomain)%yStart)
      evalper(idomain)%julStart = jday
 
      jday = julday(dd = evalper(idomain)%dEnd, mm = evalper(idomain)%mEnd, yy = evalper(idomain)%yEnd)
      evalper(idomain)%julEnd = jday
 
      ! determine warming period
      warmper(idomain)%julStart = evalper(idomain)%julStart - warmingdays(idomain)
      warmper(idomain)%julEnd = evalper(idomain)%julStart - 1
 
      call caldat(warmper(idomain)%julStart, dd = warmper(idomain)%dStart, mm = warmper(idomain)%mStart, &
                  yy = warmper(idomain)%yStart)
      call caldat(warmper(idomain)%julEnd, dd = warmper(idomain)%dEnd, mm = warmper(idomain)%mEnd, &
                  yy = warmper(idomain)%yEnd)
 
      ! simulation Period = warming Period + evaluation Period
      simper(idomain)%dStart = warmper(idomain)%dStart
      simper(idomain)%mStart = warmper(idomain)%mStart
      simper(idomain)%yStart = warmper(idomain)%yStart
      simper(idomain)%julStart = warmper(idomain)%julStart
      simper(idomain)%dEnd = evalper(idomain)%dEnd
      simper(idomain)%mEnd = evalper(idomain)%mEnd
      simper(idomain)%yEnd = evalper(idomain)%yEnd
      simper(idomain)%julEnd = evalper(idomain)%julEnd
    end do
 
    call set_land_cover_scenes_id(simper, lcyearid)
 
    !===============================================================
    ! Settings for Optimization
    !===============================================================
    ! namelist for Optimization settings
    call position_nml('Optimization', unamelist)
    read(unamelist, nml = optimization)
    ! checking of settings and default value initialization moved to new subroutine
    ! because global_parameters need to be set, which is not the case right now
    call close_nml(unamelist)
 
  end subroutine common_mhm_mrm_read_config
 
 
  !> \brief check optimization settings
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine check_optimization_settings
 
    use mo_common_mhm_mrm_variables, only : dds_r, niterations, sce_ngs, sce_npg, sce_nps
    use mo_common_variables, only : global_parameters
 
    implicit none
 
    integer(i4) :: n_true_pars
 
 
    ! check and set default values
    if (niterations .le. 0_i4) then
      call error_message('Number of iterations for Optimization (nIterations) must be greater than zero')
    end if
    if (dds_r .lt. 0.0_dp .or. dds_r .gt. 1.0_dp) then
      call error_message('dds_r must be between 0.0 and 1.0')
    end if
    if (sce_ngs .lt. 1_i4) then
      call error_message('number of complexes in SCE (sce_ngs) must be at least 1')
    end if
    ! number of points in each complex: default = 2n+1
    if (sce_npg .lt. 0_i4) then
      n_true_pars = count(nint(global_parameters(:, 4)) .eq. 1)
      sce_npg = 2 * n_true_pars + 1_i4
    end if
    ! number of points in each sub-complex: default = n+1
    if (sce_nps .lt. 0_i4) then
      n_true_pars = count(nint(global_parameters(:, 4)) .eq. 1)
      sce_nps = n_true_pars + 1_i4
    end if
    if (sce_npg .lt. sce_nps) then
      call error_message('number of points per complex (sce_npg) must be greater or', raise=.false.)
      call error_message('equal number of points per sub-complex (sce_nps)')
    end if
 
  end subroutine check_optimization_settings
 
 
  !> \brief check resolution
  !> \authors Robert Schweppe
  !> \date Jun 2018
  subroutine common_check_resolution(do_message, allow_subgrid_routing)
 
    use mo_common_mhm_mrm_variables, only : resolutionrouting
    use mo_common_variables, only : domainmeta, resolutionhydrology
    use mo_string_utils, only : num2str
 
    implicit none
 
    logical, intent(in) :: do_message !< flag to print messages
    logical, intent(in) :: allow_subgrid_routing !< flag to allow subgrid routing
 
    integer(i4) :: idomain, domainid
 
    ! conversion factor L11 to L1
    real(dp) :: cellfactorrbyh
 
 
    !===============================================================
    ! check matching of resolutions: hydrology, forcing and routing
    !===============================================================
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      cellfactorrbyh = resolutionrouting(idomain) / resolutionhydrology(idomain)
      if (do_message) then
        call message()
        call message('domain ', trim(adjustl(num2str(domainid))), ': ')
        call message('resolution Hydrology (domain ', trim(adjustl(num2str(domainid))), ')     = ', &
                trim(adjustl(num2str(resolutionhydrology(idomain)))))
        call message('resolution Routing (domain ', trim(adjustl(num2str(domainid))), ')       = ', &
                trim(adjustl(num2str(resolutionrouting(idomain)))))
      end if
      !
      if(nint(cellfactorrbyh * 100.0_dp) .eq. 100) then
        if (do_message) then
          call message()
          call message('Resolution of routing and hydrological modeling are equal!')
        end if
 
      else if ((nint(cellfactorrbyh * 100.0_dp) .gt. 100) .and. .not.allow_subgrid_routing) then
        if(nint(mod(cellfactorrbyh, 2.0_dp) * 100.0_dp) .ne. 0) then
          call error_message('***ERROR: Resolution of routing is not a multiple of hydrological model resolution!', raise=.false.)
          call error_message('   FILE: mhm.nml, namelist: mainconfig, variable: resolutionRouting')
        end if
        !
        if (do_message) then
          call message()
          call message('Resolution of routing is bigger than hydrological model resolution by ', &
                  trim(adjustl(num2str(nint(cellfactorrbyh)))), ' times !')
        end if
      end if
      !
    end do
 
  end subroutine common_check_resolution
 
 
  ! ToDo: make this a procedure of period
  !> \brief copy period data
  subroutine period_copy_period_data(toPeriod, fromPeriod)
    use mo_common_types, only: period
    type(period), intent(inout) :: toPeriod !< copy to this period
    type(period), intent(in)    :: fromPeriod !< copy from this period
 
    toperiod%dStart   = fromperiod%dStart    ! first day
    toperiod%mStart   = fromperiod%mStart    ! first month
    toperiod%yStart   = fromperiod%yStart    ! first year
    toperiod%dEnd     = fromperiod%dEnd      ! last  day
    toperiod%mEnd     = fromperiod%mEnd      ! last  month
    toperiod%yEnd     = fromperiod%yEnd      ! last  year
    toperiod%julStart = 0 ! first julian day
    toperiod%julEnd   = 0 ! last  julian day
    toperiod%nObs     = 0 ! total number of observations
 
  end subroutine period_copy_period_data
 
END MODULE mo_common_mhm_mrm_read_config