mo_mrm_write.f90

Go to the documentation of this file.

!> \file mo_mrm_write.f90
!> \brief   \copybrief mo_mrm_write
!> \details \copydetails mo_mrm_write
 
!> \brief write of discharge and restart files
!> \details This module contains the subroutines for writing the discharge files and optionally the restart files.
!> \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_write
 
  use mo_kind, only : i4, dp
  use mo_mrm_write_fluxes_states, only : outputdataset
  use mo_mrm_global_variables, only : output_time_reference_mrm
  use mo_message, only : message, error_message
 
  implicit none
 
  public :: mrm_write
  public :: mrm_write_optinamelist
  public :: mrm_write_optifile
 
  private
 
contains
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_write
 
  !    PURPOSE
  !>       \brief write discharge and restart files
 
  !>       \details First, this subroutine calls the writing or restart files that only
  !>       succeeds if it happens after the write of mHM restart files because
  !>       mHM restart files must exist. Second, simulated discharge is aggregated to the daily
  !>       scale and then written to file jointly with observed discharge
 
  !    HISTORY
  !>       \authors Juliane Mai, Rohini Kumar & Stephan Thober
 
  !>       \date Aug 2015
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
  ! Pallav Shrestha, Husain Najafi Mar 2022 - refactoring for measured timestep output
 
  subroutine mrm_write
 
    use mo_common_constants, only : nodata_dp
    use mo_common_mhm_mrm_variables, only : evalper, mrm_coupling_mode, ntstepday, simper, warmingdays
    use mo_common_variables, only : mrmfilerestartout, domainmeta, write_restart
    use mo_mrm_global_variables, only : domain_mrm, &
                                        gauge, mrm_runoff, ngaugestotal, nmeasperday
    use mo_mrm_restart, only : mrm_write_restart
 
    implicit none
 
    integer(i4) :: domainid, idomain
 
    integer(i4) :: iday, isubday, is, ie
 
    integer(i4) :: maxdailytimesteps, maxmeastimesteps
 
    integer(i4) :: tt
 
    integer(i4) :: gg
 
    integer(i4) :: ntimesteps
 
    real(dp), dimension(:, :), allocatable :: d_qmod, subd_qmod ! Sim discharge at daily and subdaily time step
    real(dp), dimension(:, :), allocatable :: d_qobs, subd_qobs ! Obs discharge at daily and subdaily time step
 
    ! between simulated and measured time scale
    integer(i4) :: factor
 
    ! simulated Timesteps per Day
    integer(i4) :: tpd_sim
 
    ! observed Timesteps per Day
    integer(i4) :: tpd_obs
 
    ! --------------------------------------------------------------------------
    ! WRITE CONFIG FILE
    ! --------------------------------------------------------------------------
    if (mrm_coupling_mode .eq. 0_i4) call write_configfile()
 
    ! --------------------------------------------------------------------------
    ! WRITE RESTART
    ! --------------------------------------------------------------------------
    if (write_restart) then
      do idomain = 1, domainmeta%nDomains
        domainid = domainmeta%indices(idomain)
        if (domainmeta%doRouting(idomain)) call mrm_write_restart(idomain, domainid, mrmfilerestartout)
      end do
    end if
 
    ! --------------------------------------------------------------------------
    ! STORE DAILY DISCHARGE TIMESERIES OF EACH GAUGING STATION
    ! FOR SIMULATIONS DURING THE EVALUATION PERIOD
    !
    !  **** AT DAILY TIME STEPS ****
    ! Note:: Observed Q are stored only for the evaluation period and not for
    !        the warming days
    ! --------------------------------------------------------------------------
 
    ! copy time resolution to local variables
    tpd_sim = ntstepday
    tpd_obs = nmeasperday
 
    ! check if modelled timestep is an integer multiple of measured timesteps
    if (modulo(tpd_sim, tpd_obs) .eq. 0) then
      factor = tpd_sim / tpd_obs
    else
      call error_message(' Error: Number of modelled datapoints is no multiple of measured datapoints per day')
    end if
 
 
    maxdailytimesteps = maxval(evalper(1 : domainmeta%nDomains)%julEnd - evalper(1 : domainmeta%nDomains)%julStart + 1)
    maxmeastimesteps  = maxval(evalper(1 : domainmeta%nDomains)%julEnd - evalper(1 : domainmeta%nDomains)%julStart + 1) * tpd_obs
    allocate(d_qmod(maxdailytimesteps, ngaugestotal))
    allocate(d_qobs(maxdailytimesteps, ngaugestotal))
    allocate(subd_qmod( maxmeastimesteps, ngaugestotal))
    allocate(subd_qobs( maxmeastimesteps, ngaugestotal))
    d_qmod     = nodata_dp
    d_qobs     = nodata_dp
    subd_qmod  = nodata_dp
    subd_qobs  = nodata_dp
 
 
    ! loop over domains
    do idomain = 1, domainmeta%nDomains
      if (domainmeta%doRouting(idomain)) then
        domainid = domainmeta%indices(idomain)
 
        ! Convert simulated values to daily
        ntimesteps = (simper(idomain)%julEnd - simper(idomain)%julStart + 1) * ntstepday
        iday  = 0
        do tt = warmingdays(idomain) * ntstepday + 1, ntimesteps, ntstepday
          is = tt
          ie = tt + ntstepday - 1
          iday = iday + 1
          ! over gauges
          do gg = 1, domain_mrm(idomain)%nGauges
            ! simulation
            d_qmod(iday, domain_mrm(idomain)%gaugeIndexList(gg)) = &
                    sum(mrm_runoff(is : ie, domain_mrm(idomain)%gaugeIndexList(gg))) / real(ntstepday, dp)
          end do
        end do
 
        dailycheck: if (nmeasperday > 1) then
          ! Convert observed values to daily
          ntimesteps = (simper(idomain)%julEnd - simper(idomain)%julStart + 1) * nmeasperday
          iday  = 0
          do tt = 1, ntimesteps, nmeasperday
            is = tt
            ie = tt + nmeasperday - 1
            iday = iday + 1
            ! over gauges
            do gg = 1, domain_mrm(idomain)%nGauges
              ! when -9999 value/s are present in current day, daily value remain -9999.
              if (.not.(any(gauge%Q(is : ie, domain_mrm(idomain)%gaugeIndexList(gg)) == nodata_dp))) then
                ! observation
                d_qobs(iday, domain_mrm(idomain)%gaugeIndexList(gg)) = &
                        sum( gauge%Q(is : ie, domain_mrm(idomain)%gaugeIndexList(gg))) / real(nmeasperday, dp)
              end if
            end do
          end do
        else
          ! observed values are already at daily (nMeasPerDay = 1) and stored for evalper
          ! over gauges
          do gg = 1, domain_mrm(idomain)%nGauges
            ! observation
            d_qobs(:, domain_mrm(idomain)%gaugeIndexList(gg)) = gauge%Q(:, domain_mrm(idomain)%gaugeIndexList(gg))
          end do
        end if dailycheck
 
 
        subdailycheck: if (nmeasperday > 1) then
 
          ! Convert simulated values to subdaily
          ntimesteps = (simper(idomain)%julEnd - simper(idomain)%julStart + 1) * ntstepday
          isubday = 0
          do tt = warmingdays(idomain) * ntstepday + 1, ntimesteps, factor
            is = tt
            ie = tt + factor - 1
            isubday = isubday + 1
            ! over gauges
            do gg = 1, domain_mrm(idomain)%nGauges
              ! simulation
              subd_qmod(isubday, domain_mrm(idomain)%gaugeIndexList(gg)) = &
                      sum(mrm_runoff(is : ie, domain_mrm(idomain)%gaugeIndexList(gg))) / real(factor, dp)
            end do
          end do
 
          ! Convert observed values to subdaily
 
          ! observed values are already at subdaily (nMeasPerDay) and stored for evalper
          ! over gauges
          do gg = 1, domain_mrm(idomain)%nGauges
            ! observation
            subd_qobs(:, domain_mrm(idomain)%gaugeIndexList(gg)) = gauge%Q(:, domain_mrm(idomain)%gaugeIndexList(gg))
          end do
 
        end if subdailycheck
 
      end if
    end do
 
 
    ! write in an ASCII file          ! OBS[nModeling_days X nGauges_total] , SIM[nModeling_days X nGauges_total]
    ! ToDo: is this if statement reasonable
    if (allocated(gauge%Q)) call write_daily_obs_sim_discharge(d_qobs(:, :), d_qmod(:, :))
 
    ! write in an ASCII file          ! OBS[nMeasPerDay X nGauges_total] , SIM[nMeasPerDay X nGauges_total]
    if (nmeasperday > 1 .and. allocated(gauge%Q)) call write_subdaily_obs_sim_discharge(subd_qobs(:, :), subd_qmod(:, :), factor)
    ! The subdaily routine is only called if subdaily Q data is provided
 
    ! free space
    deallocate(d_qmod, d_qobs, subd_qmod, subd_qobs)
  end subroutine mrm_write
  ! ------------------------------------------------------------------
 
  !    NAME
  !        write_configfile
 
  !    PURPOSE
  !>       \brief This modules writes the results of the configuration into an ASCII-file
 
  !>       \details TODO: add description
 
  !    HISTORY
  !>       \authors Christoph Schneider
 
  !>       \date May 2013
 
  ! Modifications:
  ! Juliane Mai    May 2013 - module version and documentation
  ! Stephan Thober Jun 2014 - bug fix in L11 config print out
  ! Stephan Thober Jun 2014 - updated read_restart
  ! Rohini, Luis   Jul 2015 - updated version, L1 level prints
  ! Stephan Thober Sep 2015 - updated write of stream network
  ! Stephan Thober Nov 2016 - adapted write to selected case for routing process
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  Subroutine write_configfile
 
    use mo_common_constants, only : nodata_dp
    use mo_common_file, only : file_config, uconfig
    use mo_common_mhm_mrm_variables, only : lcyearid, simper, evalper, mrm_coupling_mode, read_restart, &
                                            resolutionrouting, timestep, warmper
    use mo_common_variables, only : lc_year_end, lc_year_start, lcfilename, &
                                    dirconfigout, dirlcover, dirmorpho, dirout, mrmfilerestartout, global_parameters, &
                                    global_parameters_name, level0, level1, domainmeta, nlcoverscene, processmatrix, &
                                    resolutionhydrology, write_restart
    use mo_kind, only : dp, i4
    use mo_mrm_file, only : version
    use mo_mrm_global_variables, only : inflowgauge, l11_l1_id, l11_fromn, l11_label, &
                                        l11_length, l11_netperm, l11_rorder, l11_slope, l11_ton, l1_l11_id, domain_mrm, &
                                        dirgauges, dirtotalrunoff, gauge, level11, ngaugestotal, ninflowgaugestotal
    use mo_string_utils, only : num2str
    use mo_utils, only : ge
    use mo_os, only : check_path_isdir
 
    implicit none
 
    character(256) :: fName
 
    integer(i4) :: i, iDomain, domainID, j
 
    integer(i4) :: err
 
 
    fname = trim(adjustl(dirconfigout)) // trim(adjustl(file_config))
    call message()
    call message('  Log-file written to ', trim(fname))
    !checking whether the directory exists where the file shall be created or opened
    call check_path_isdir(trim(adjustl(dirconfigout)), raise=.true.)
    open(uconfig, file = fname, status = 'unknown', action = 'write', iostat = err)
    if (err .ne. 0) then
      call error_message('  Problems while creating File. Error-Code ', num2str(err))
    end if
    write(uconfig, 200)
    write(uconfig, 100) 'mRM-UFZ v-' // trim(version)
    write(uconfig, 100) 'S. Thober, L. Samaniego & R. Kumar, UFZ'
    write(uconfig, 200)
    write(uconfig, 100)
    write(uconfig, 201) '         M A I N  mRM  C O N F I G U R A T I O N  I N F O R M A T I O N         '
    write(uconfig, 100)
    write(uconfig, 103) 'Number of domains            ', domainmeta%nDomains
    write(uconfig, 103) 'Total No. of gauges         ', ngaugestotal
    write(uconfig, 103)    'Time Step [h]               ', timestep
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      write(uconfig, 103) 'Total No. of nodes          ', level11(idomain)%nCells
      write(uconfig, 103) 'No. of cells L0             ', level0(domainmeta%L0DataFrom(idomain))%nCells
      write(uconfig, 103) 'No. of cells L1             ', level1(idomain)%nCells
      write(uconfig, 103) 'No. of cells L11            ', level11(idomain)%nCells
 
      !    select case (iFlag_cordinate_sys)
      !    case (0)
      write(uconfig, 301)      'domain  ', domainid, '   Hydrology Resolution [m]      ', resolutionhydrology(idomain)
      write(uconfig, 301)   'domain  ', domainid, '   Routing Resolution [m]        ', resolutionrouting(idomain)
      !    case(1)
      !      write(uconfig, 302)       'domain  ',domainID, '   Hydrology Resolution [o]      ', resolutionHydrology(iDomain)
      !      write(uconfig, 302)   'domain  ',domainID, '   Routing Resolution [o]        ', resolutionRouting(iDomain)
      !    end select
    end do
    write(uconfig, 126)    'Flag READ  restart            ', read_restart
    write(uconfig, 126)    'Flag WRITE restart            ', write_restart
    !
    !******************
    ! Model Run period
    !******************
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      write(uconfig, 115) '                      Model Run Periods for domain ', num2str(domainid)
      write(uconfig, 116) &
              'From                To', &
              '   Day Month  Year   Day Month  Year'
      write(uconfig, 117)  &
              'Warming Period (1)            ', &
              warmper(idomain)%dStart, warmper(idomain)%mStart, warmper(idomain)%yStart, &
              warmper(idomain)%dEnd, warmper(idomain)%mEnd, warmper(idomain)%yEnd
      write(uconfig, 117)  &
              'Evaluation Period (2)         ', &
              evalper(idomain)%dStart, evalper(idomain)%mStart, evalper(idomain)%yStart, &
              evalper(idomain)%dEnd, evalper(idomain)%mEnd, evalper(idomain)%yEnd
      write(uconfig, 117)  &
              'Simulation Period (1)+(2)     ', &
              simper(idomain)%dStart, simper(idomain)%mStart, simper(idomain)%yStart, &
              simper(idomain)%dEnd, simper(idomain)%mEnd, simper(idomain)%yEnd
    end do
 
    !*********************************
    ! Model Land Cover Observations
    !*********************************
    if (processmatrix(8, 1) .eq. 1) then
      do idomain = 1, domainmeta%nDomains
        domainid = domainmeta%indices(idomain)
        write(uconfig, 118) '       Land Cover Observations for domain ', num2str(domainid)
        write(uconfig, 119) ' Start Year', ' End Year', '    Land cover scene', 'Land Cover File'
        do i = 1, nlcoverscene
          write(uconfig, 120) lc_year_start(i), lc_year_end(i), &
                  lcyearid(max(evalper(idomain)%yStart, lc_year_start(i)), idomain), trim(lcfilename(i))
        end do
      end do
    end if
    !*********************************
    ! Initial Parameter Ranges
    !*********************************
    write(uconfig, 121) '  Initial Transfer Function Parameter Ranges (gammas)  '
    !
    ! Transfer functions
    write(uconfig, 122)      &
            '         i', '            min', '            max', '        current', &
            '                               name'
    do i = 1, size(global_parameters, 1)
      write(uconfig, 123) &
              i, global_parameters(i, 1), global_parameters(i, 2), global_parameters(i, 3), &
              trim(adjustl(global_parameters_name(i)))
    end do
    ! domain runoff data
    write(uconfig, 202) '                domain Runoff Data                '
    write(uconfig, 107) ' Gauge No.', '  domain Id', '     Qmax[m3/s]', '     Qmin[m3/s]'
    do i = 1, ngaugestotal
      if(any(gauge%Q(:, i) > nodata_dp)) then
        write(uconfig, 108) i, gauge%domainId(i), maxval(gauge%Q(:, i), gauge%Q(:, i) > nodata_dp), &
                minval(gauge%Q(:, i), gauge%Q(:, i) > nodata_dp)
      else
        write(uconfig, 108) i, gauge%domainId(i), nodata_dp, nodata_dp
      end if
    end do
    ! inflow gauge data
    if (ninflowgaugestotal .GT. 0) then
      write(uconfig, 202) '                domain Inflow Data                 '
      write(uconfig, 107) ' Gauge No.', '  domain Id', '     Qmax[m3/s]', '     Qmin[m3/s]'
      do i = 1, ninflowgaugestotal
        if(all(inflowgauge%Q(:, i) > nodata_dp)) then
          write(uconfig, 108) i, inflowgauge%domainId(i), maxval(inflowgauge%Q(:, i), inflowgauge%Q(:, i) > nodata_dp), &
                  minval(inflowgauge%Q(:, i), inflowgauge%Q(:, i) > nodata_dp)
        else
          write(uconfig, 108) i, inflowgauge%domainId(i), nodata_dp, nodata_dp
        end if
      end do
    end if
    ! domain config
    write(uconfig, 218) 'domain-wise Configuration'
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      write(uconfig, 103) 'domain No.                   ', domainid, &
              'No. of gauges               ', domain_mrm(idomain)%nGauges
 
      write(uconfig, 222)   'Directory list'
 
      write(uconfig, 224) 'Directory to morphological input         ', dirmorpho(idomain)
      write(uconfig, 224) 'Directory to land cover input            ', dirlcover(idomain)
      write(uconfig, 224) 'Directory to gauging station input       ', dirgauges(idomain)
      if (mrm_coupling_mode .eq. 0) then
        write(uconfig, 224) 'Directory to simulated runoff input      ', dirtotalrunoff(idomain)
      end if
      write(uconfig, 224) 'Directory to write output by default     ', dirout(idomain)
      write(uconfig, 224) 'File to write mRM output when restarted  ', mrmfilerestartout(idomain)
 
      write(uconfig, 102) 'River Network  (Routing level)'
      write(uconfig, 100) 'Label 0 = intermediate draining cell '
      write(uconfig, 100) 'Label 1 = headwater cell             '
      write(uconfig, 100) 'Label 2 = sink cell                  '
 
      if (processmatrix(8, 1) .eq. 1_i4) then
        write(uconfig, 104) '   Overall', &
                '      From', &
                '        To', &
                '   Routing', &
                '     Label', &
                '    Length', &
                '      Mean', &
                '      Link', &
                '   Routing', &
                '   Routing', &
                '  Sequence', &
                '          ', &
                '          ', &
                '     Slope'
        !
        write(uconfig, 105) '        Id', &
                '      Node', &
                '      Node', &
                '', &
                '', &
                '      [km]', &
                '    [o/oo]'
        !
        do j = level11(idomain)%iStart, level11(idomain)%iEnd - 1
          i = l11_netperm(j) + level11(idomain)%iStart - 1 ! adjust permutation for multi-domain option
          write(uconfig, 106) i, l11_fromn(i), l11_ton(i), l11_rorder(i), l11_label(i), &
                  l11_length(i) / 1000.0_dp, l11_slope(i) * 1.0e3_dp
        end do
 
      else if (processmatrix(8, 1) .eq. 2_i4) then
        write(uconfig, 134) '   Overall', &
                '      From', &
                '        To', &
                '   Routing', &
                '     Label', &
                '      Link', &
                '   Routing', &
                '   Routing', &
                '  Sequence', &
                '          '
        !
        write(uconfig, 135) '        Id', &
                '      Node', &
                '      Node', &
                '', &
                ''
        !
        do j = level11(idomain)%iStart, level11(idomain)%iEnd - 1
          i = l11_netperm(j) + level11(idomain)%iStart - 1 ! adjust permutation for multi-domain option
          write(uconfig, 136) i, l11_fromn(i), l11_ton(i), l11_rorder(i), l11_label(i)
        end do
      end if
      ! draining node at L11
      write(uconfig, 109)  '   Overall', '     domain', &
              '      Cell', '   Routing', &
              '        Id', '   Node Id'
      do i = level11(idomain)%Id(1), level11(idomain)%Id(level11(idomain)%nCells)
        write(uconfig, 110) i + level11(idomain)%iStart - 1, i
      end do
 
      ! L1 level information
      write(uconfig, 111)  '  Modeling', '   Routing', ' Effective', &
              '      Cell', '   Cell Id', '      Area', &
              '        Id', '       [-]', '     [km2]'
      if (ge(resolutionrouting(idomain), resolutionhydrology(idomain))) then
        do i = level1(idomain)%Id(1), level1(idomain)%Id(level1(idomain)%nCells)
          write(uconfig, 113) i + level1(idomain)%iStart - 1, l1_l11_id(i + level1(idomain)%iStart - 1), &
            level1(idomain)%CellArea(i) * 1.e-6_dp
      domainid = domainmeta%indices(idomain)
        end do
      else
        do i = level11(idomain)%Id(1), level11(idomain)%Id(level11(idomain)%nCells)
          write(uconfig, 110) i + level11(idomain)%iStart - 1, l11_l1_id(i + level11(idomain)%iStart - 1)
        end do
      end if
      write(uconfig, 114)  ' Total[km2]', sum(level1(idomain)%CellArea) * 1.e-6_dp
    end do
 
    write(uconfig, *)
    close(uconfig)
 
    !! Formats
    100 format (a80)
    102 format (/ 30('-') / a30 / 30('-'))
    103 format (a20, 10x, i10)
    104 format (/ 75('-') / 5a10, 5x, 2a10 / 5a10, 5x, 2a10)
    105 format (5a10, 5x, 2a10 / 75('-'))
    106 format (5i10, 5x, 2f10.3)
    107 format (2a10, 2a15)
    108 format (2i10, 2f15.3)
    !
    109 format (/ 20('-') / 2a10 / 2a10 / 2a10 / 20('-'))
    110 format (2i10)
    !
    111 format (/ 30('-') / 3a10 / 3a10 / 3a10 /  30('-'))
    113 format (            2i10,   1f10.3         )
    114 format (30('-') / a15, 5x, 1f10.3 /)
    !
    115 format (/61('-')/ a50, a10 /61('-'))
    116 format (39x, a22 / 25x, a36)
    117 format (3(a25, 6(i6)))
    !
    118 format (/50('-')/ a40, a10  /50('-'))
    119 format (a10, a10, a20, a20/)
    120 format (i10, i10, 10x, i10, a20)
    !
    121 format (/55('-')/ a55 /55('-'))
    122 format (a10, 3a15, a35)
    123 format (i10, 3f15.3, a35)
    !
    126 format (a30, 9x, l1)
    !
    134 format (/ 50('-') / 5a10 / 5a10)
    135 format (5a10 / 50('-'))
    136 format (5i10)
    !
    200 format (80('-'))
    201 format (a80)
    202 format (/50('-')/ a50 /50('-'))
    !
    218 format (/ 80('-')/ 26x, a24, 26x, /80('-'))
    222 format (/80('-')/ 26x, a21 /80('-'))
    224 format (a40, 5x, a256)
 
    301 format (a7, i2, a32, f15.0)
    ! 302 format (a7, i2, a32,es20.8)
  end Subroutine write_configfile
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        write_daily_obs_sim_discharge
 
  !    PURPOSE
  !>       \brief Write a file for the daily observed and simulated discharge timeseries
  !>       during the evaluation period for each gauging station
 
  !>       \details Write a file for the daily observed and simulated discharge timeseries
  !>       during the evaluation period for each gauging station
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(:, :) :: Qobs" daily time series of observed dischargedims = (nModeling_days
  !>       , nGauges_total)
  !>       \param[in] "real(dp), dimension(:, :) :: Qsim" daily time series of modeled dischargedims = (nModeling_days ,
  !>       nGauges_total)
 
  !    HISTORY
  !>       \authors Rohini Kumar
 
  !>       \date August 2013
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine write_daily_obs_sim_discharge(Qobs, Qsim)
 
    use mo_common_constants, only : nodata_dp
    use mo_common_mhm_mrm_variables, only : evalper
    use mo_common_variables, only : dirout, domainmeta
    use mo_errormeasures, only : kge, nse
    use mo_julian, only : dec2date
    use mo_mrm_file, only : file_daily_discharge, ncfile_discharge, udaily_discharge
    use mo_mrm_global_variables, only : domain_mrm, gauge
    use mo_string_utils, only : num2str
    use mo_utils, only : ge
    use mo_netcdf, only : ncdataset, ncdimension, ncvariable
 
    implicit none
 
    ! daily time series of observed dischargedims = (nModeling_days , nGauges_total)
    real(dp), dimension(:, :), intent(in) :: Qobs
 
    ! daily time series of modeled dischargedims = (nModeling_days , nGauges_total)
    real(dp), dimension(:, :), intent(in) :: Qsim
 
    character(256) :: fName, formHeader, formData, dummy
 
    integer(i4) :: domainID, iDomain, gg, tt, err
 
    integer(i4) :: igauge_start, igauge_end
 
    integer(i4) :: day, month, year
 
    integer(i4) :: tlength
 
    ! time axis
    integer(i4), allocatable, dimension(:) :: taxis
 
    real(dp) :: newTime
 
    ! nc related variables
    type(ncdataset) :: nc_out
    type(ncdimension) :: dim, dim_bnd
    type(ncvariable) :: var
 
    ! initalize igauge_start
    igauge_start = 1
 
    ! domain loop
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      if(domain_mrm(idomain)%nGauges .lt. 1) cycle
 
      ! estimate igauge_end
      igauge_end = igauge_start + domain_mrm(idomain)%nGauges - 1
 
      ! check the existance of file
      fname = trim(adjustl(dirout(idomain))) // trim(adjustl(file_daily_discharge))
      open(udaily_discharge, file = trim(fname), status = 'unknown', action = 'write', iostat = err)
      if(err .ne. 0) then
        call error_message('  IOError while openening "', trim(fname), '". Error-Code ', num2str(err))
      end if
 
      ! header
      write(formheader, *) '( 4a8, ', domain_mrm(idomain)%nGauges, '(2X, a5, i10.10, 2X, a5, i10.10) )'
      write(udaily_discharge, formheader) 'No', 'Day', 'Mon', 'Year', &
              ('Qobs_', gauge%gaugeId(gg), &
                      'Qsim_', gauge%gaugeId(gg), gg = igauge_start, igauge_end)
 
      ! form data
      write(formdata, *) '( 4I8, ', domain_mrm(idomain)%nGauges, '(2X,   f15.7 , 2X,  f15.7  ) )'
 
      ! write data
      newtime = real(evalper(idomain)%julStart, dp) - 0.5_dp
 
      do tt = 1, (evalper(idomain)%julEnd - evalper(idomain)%julStart + 1)
        call dec2date(newtime, yy = year, mm = month, dd = day)
        write(udaily_discharge, formdata) tt, day, month, year, (qobs(tt, gg), qsim(tt, gg), gg = igauge_start, igauge_end)
        newtime = newtime + 1.0_dp
      end do
 
      ! close file
      close(udaily_discharge)
 
      ! ======================================================================
      ! write netcdf file
      ! ======================================================================
      fname = trim(adjustl(dirout(idomain))) // trim(adjustl(ncfile_discharge))
      nc_out = ncdataset(trim(fname), "w")
      tlength = evalper(idomain)%julEnd - evalper(idomain)%julStart + 1
      ! write time
      allocate(taxis(tlength))
 
      ! tt is dependent on the unit of the time axis and is set to hours in mRM
      select case( output_time_reference_mrm)
        case(0)
          forall(tt = 1 : tlength) taxis(tt) = (tt-1) * 24
        case(1)
          forall(tt = 1 : tlength) taxis(tt) = tt * 24 - 12
        case(2)
          forall(tt = 1 : tlength) taxis(tt) = tt * 24
      end select
 
      call dec2date(real(evalper(idomain)%julStart, dp) - 0.5_dp, yy = year, mm = month, dd = day)
      dim = nc_out%setDimension("time", tlength)
      var = nc_out%setVariable("time", "i32", [dim])
      call var%setData(taxis)
      call var%setAttribute( &
        "units", &
        'hours since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' 00:00:00' &
      )
      call var%setAttribute("long_name", "time in hours")
      call var%setAttribute("bounds", "time_bnds")
      call var%setAttribute("axis", "T")
      dim_bnd = nc_out%setDimension("bnds", 2)
      var = nc_out%setVariable("time_bnds", "i32", [dim_bnd, dim])
      do tt = 1, tlength
        call var%setData((tt - 1) * 24, (/1, tt/))
        call var%setData(tt * 24, (/2, tt/))
      end do
      deallocate(taxis)
      ! write gauges
      do gg = igauge_start, igauge_end
        var = nc_out%setVariable('Qsim_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
        call var%setFillValue(nodata_dp)
        call var%setData(qsim(1 : tlength, gg))
        call var%setAttribute("units", "m3 s-1")
        call var%setAttribute("long_name", 'simulated discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
        call var%setAttribute("missing_value", nodata_dp)
        ! write observed discharge at that gauge
        var = nc_out%setVariable('Qobs_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
        call var%setFillValue(nodata_dp)
        call var%setData(qobs(1 : tlength, gg))
        call var%setAttribute("units", "m3 s-1")
        call var%setAttribute("long_name", 'observed discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
        call var%setAttribute("missing_value", nodata_dp)
      end do
      call nc_out%close()
 
      ! ======================================================================
      ! screen output
      ! ======================================================================
 
      ! if ( nMeasPerDay == 1_i4 ) then ! only print daily stats for daily Qobs
 
      call message()
      write(dummy, '(I3)') domainid
      call message('  OUTPUT: saved daily discharge file for domain ', trim(adjustl(dummy)))
      call message('    to ', trim(fname))
      do gg = igauge_start, igauge_end
        if (count(ge(qobs(:, gg), 0.0_dp)) > 1)  then
          call message('    KGE of daily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
                  trim(adjustl(num2str(kge(qobs(:, gg), qsim(:, gg), mask = (ge(qobs(:, gg), 0.0_dp)))))))
          call message('    NSE of daily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
                  trim(adjustl(num2str(nse(qobs(:, gg), qsim(:, gg), mask = (ge(qobs(:, gg), 0.0_dp)))))))
        end if
      end do
 
      ! end if
 
      ! update igauge_start
      igauge_start = igauge_end + 1
      !
    end do
    !
  end subroutine write_daily_obs_sim_discharge
 
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        write_subdaily_obs_sim_discharge
 
  !    PURPOSE
  !>       \brief Write a file for the simulated discharge timeseries
  !>       during the evaluation period for each gauging station
 
  !>       \details Write a file for the simulated discharge timeseries
  !>       during the evaluation period for each gauging station
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(:, :) :: Qobs"  time series of observed discharge dims = (nMeasTimeSteps ,
  !>       nGauges_total)
  !>       \param[in] "real(dp), dimension(:, :) :: Qsim"  time series of modeled discharge dims = (nMeasTimeSteps ,
  !>       nGauges_total)
  !>       \param[in] "integer(i4),              :: factor" ratio of modelled time steps per day to observation time steps per day
 
  !    HISTORY
  !>       \authors Rohini Kumar
 
  !>       \date August 2013
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
  ! Pallav Shrestha Jul 2021 - ported for printing out simulations at hourly time step
  ! Pallav Shrestha, Husain Najafi Mar 2022 - refactoring for subdaily timestep output
 
  subroutine write_subdaily_obs_sim_discharge(Qobs, Qsim, factor)
 
    use mo_common_constants, only : nodata_dp
    use mo_common_mhm_mrm_variables, only : evalper
    use mo_common_variables, only : dirout, domainmeta
    use mo_errormeasures, only : kge, nse
    use mo_julian, only : dec2date
    use mo_mrm_file, only : ncfile_subdaily_discharge, file_subdaily_discharge, &
                            usubdaily_discharge
    use mo_mrm_global_variables, only : domain_mrm, gauge, nmeasperday
    use mo_string_utils, only : num2str
    use mo_utils, only : ge
    use mo_netcdf, only : ncdataset, ncdimension, ncvariable
 
    implicit none
 
    ! time series of observed discharge. dims = (nMeasTimeSteps , nGauges_total)
    real(dp), dimension(:, :), intent(in) :: Qobs
 
    ! time series of modeled discharge. dims = (nMeasTimeSteps , nGauges_total)
    real(dp), dimension(:, :), intent(in) :: Qsim
 
    ! ratio of modelled time steps per day to observation time steps per day
    integer(i4),               intent(in) :: factor
 
    character(256) :: fName, formHeader, formData, dummy
 
    integer(i4) :: domainID, iDomain, gg, tt, err
 
    integer(i4) :: igauge_start, igauge_end
 
    integer(i4) :: hour, day, month, year
 
    integer(i4) :: tlength
 
    ! time axis
    integer(i4), allocatable, dimension(:) :: taxis
 
    real(dp) :: newTime
 
    ! nc related variables
    type(ncdataset) :: nc_out
    type(ncdimension) :: dim, dim_bnd
    type(ncvariable) :: var
 
    ! use minutes if needed
    logical :: use_minutes
    integer(i4) :: time_unit_factor
 
    ! initalize igauge_start
    igauge_start = 1
 
    ! domain loop
    do idomain = 1, domainmeta%nDomains
      domainid = domainmeta%indices(idomain)
      if(domain_mrm(idomain)%nGauges .lt. 1) cycle
 
      ! estimate igauge_end
      igauge_end = igauge_start + domain_mrm(idomain)%nGauges - 1
 
 
      ! ======================================================================
      ! write text file
      ! ======================================================================
 
      ! check the existance of file
      fname = trim(adjustl(dirout(idomain))) // trim(adjustl(file_subdaily_discharge))
      open(usubdaily_discharge, file = trim(fname), status = 'unknown', action = 'write', iostat = err)
      if(err .ne. 0) then
        call error_message('  IOError while openening "', trim(fname), '". Error-Code ', num2str(err))
      end if
 
      ! header
      write(formheader, *) '( 5a8, ', domain_mrm(idomain)%nGauges, '(2X, a5, i10.10, 2X, a5, i10.10) )'
      write(usubdaily_discharge, formheader) 'No', 'Hour', 'Day', 'Mon', 'Year', &
              ('Qobs_', gauge%gaugeId(gg), &
                      'Qsim_', gauge%gaugeId(gg), gg = igauge_start, igauge_end)
 
      ! form data
      write(formdata, *) '( 5I8, ', domain_mrm(idomain)%nGauges, '(2X,   f15.7 , 2X,  f15.7  ) )'
 
      ! write data
      newtime = real(evalper(idomain)%julStart, dp) - 0.5_dp
 
      do tt = 1, (evalper(idomain)%julEnd - evalper(idomain)%julStart + 1) * nmeasperday
        call dec2date(newtime, yy = year, mm = month, dd = day, hh = hour)
        write(usubdaily_discharge, formdata) tt, hour, day, month, year, (qobs(tt, gg), qsim(tt, gg), gg = igauge_start, igauge_end)
        newtime = newtime + 1.0_dp / real(nmeasperday, dp)
      end do
 
      ! close file
      close(usubdaily_discharge)
 
      ! ======================================================================
      ! write netcdf file
      ! ======================================================================
      fname = trim(adjustl(dirout(idomain))) // trim(adjustl(ncfile_subdaily_discharge))
      nc_out = ncdataset(trim(fname), "w")
      tlength = (evalper(idomain)%julEnd - evalper(idomain)%julStart + 1) * nmeasperday
      ! write time
      allocate(taxis(tlength))
 
      use_minutes = .false.
      time_unit_factor = 1
      if ( mod(factor, 2) == 1 ) then
        use_minutes = .true.
        time_unit_factor = 60
      end if
 
      ! tt is dependent on the unit of the time axis and is set to hours in mRM
      select case( output_time_reference_mrm)
        case(0)
          forall(tt = 1 : tlength) taxis(tt) = (tt-1) * factor * time_unit_factor
        case(1)
          forall(tt = 1 : tlength) taxis(tt) = tt * factor - factor * time_unit_factor / 2
        case(2)
          forall(tt = 1 : tlength) taxis(tt) = tt * factor * time_unit_factor
      end select
 
      call dec2date(real(evalper(idomain)%julStart, dp) - 0.5_dp, yy = year, mm = month, dd = day, hh = hour)
      dim = nc_out%setDimension("time", tlength)
      var = nc_out%setVariable("time", "i32", [dim])
      call var%setData(taxis)
      if (use_minutes) then
        call var%setAttribute( &
          "units", &
          'minutes since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' '// &
          trim(num2str(hour, '(i2.2)'))//':00:00' &
        )
        call var%setAttribute("long_name", "time in minutes")
      else
        call var%setAttribute( &
          "units", &
          'hours since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' '// &
          trim(num2str(hour, '(i2.2)'))//':00:00' &
        )
        call var%setAttribute("long_name", "time in hours")
      end if
      call var%setAttribute("bounds", "time_bnds")
      call var%setAttribute("axis", "T")
      dim_bnd = nc_out%setDimension("bnds", 2)
      var = nc_out%setVariable("time_bnds", "i32", [dim_bnd, dim])
      do tt = 1, tlength
        call var%setData((tt - 1) * factor * time_unit_factor, (/1, tt/))
        call var%setData(tt * factor * time_unit_factor, (/2, tt/))
      end do
      deallocate(taxis)
      ! write gauges
      do gg = igauge_start, igauge_end
        var = nc_out%setVariable('Qsim_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
        call var%setFillValue(nodata_dp)
        call var%setData(qsim(1 : tlength, gg))
        call var%setAttribute("units", "m3 s-1")
        call var%setAttribute("long_name", 'simulated discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
        call var%setAttribute("missing_value", nodata_dp)
        ! write observed discharge at that gauge
        var = nc_out%setVariable('Qobs_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
        call var%setFillValue(nodata_dp)
        call var%setData(qobs(1 : tlength, gg))
        call var%setAttribute("units", "m3 s-1")
        call var%setAttribute("long_name", 'observed discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
        call var%setAttribute("missing_value", nodata_dp)
      end do
      call nc_out%close()
 
      ! ======================================================================
      ! screen output
      ! ======================================================================
      call message()
      write(dummy, '(I3)') domainid
      call message('  OUTPUT: saved subdaily discharge file for domain ', trim(adjustl(dummy)))
      call message('    to ', trim(fname))
      do gg = igauge_start, igauge_end
        if (count(ge(qobs(:, gg), 0.0_dp)) > 1)  then
          call message('    KGE of subdaily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
                  trim(adjustl(num2str(kge(qobs(:, gg), qsim(:, gg), mask = (ge(qobs(:, gg), 0.0_dp)))))))
          call message('    NSE of subdaily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
                  trim(adjustl(num2str(nse(qobs(:, gg), qsim(:, gg), mask = (ge(qobs(:, gg), 0.0_dp)))))))
        end if
      end do
 
      ! update igauge_start
      igauge_start = igauge_end + 1
      !
    end do
    !
  end subroutine write_subdaily_obs_sim_discharge
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_write_optifile
 
  !    PURPOSE
  !>       \brief Write briefly final optimization results.
 
  !>       \details Write overall best objective function and the best optimized parameter set to a file_opti.
 
  !    INTENT(IN)
  !>       \param[in] "real(dp) :: best_OF"                             best objective function value as returnedby the
  !>       optimization routine
  !>       \param[in] "real(dp), dimension(:) :: best_paramSet"         best associated global parameter setCalled only
  !>       when optimize is .TRUE.
  !>       \param[in] "character(len = *), dimension(:) :: param_names"
 
  !    HISTORY
  !>       \authors David Schaefer
 
  !>       \date July 2013
 
  ! Modifications:
  ! Rohini Kumar   Aug 2013 - change in structure of the code including call statements
  ! Juliane Mai    Oct 2013 - clear parameter names added
  !                         - double precision written
  ! Stephan Thober Oct 2015 - ported to mRM
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine mrm_write_optifile(best_OF, best_paramSet, param_names)
 
    use mo_common_mhm_mrm_file, only : file_opti, uopti
    use mo_common_variables, only : dirconfigout
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! best objective function value as returnedby the optimization routine
    real(dp), intent(in) :: best_of
 
    ! best associated global parameter setCalled only when optimize is .TRUE.
    real(dp), dimension(:), intent(in) :: best_paramset
 
    character(len = *), dimension(:), intent(in) :: param_names
 
    character(256) :: fname, formheader, formparams
 
    integer(i4) :: ii, err, n_params
 
 
    ! number of parameters
    n_params = size(best_paramset)
 
    ! open file
    fname = trim(adjustl(dirconfigout)) // trim(adjustl(file_opti))
    open(uopti, file = fname, status = 'unknown', action = 'write', iostat = err, recl = (n_params + 1) * 40)
    if(err .ne. 0) then
      call error_message('  IOError while openening "', trim(fname), '" Error-Code ', num2str(err))
    end if
 
    ! header
    write(formheader, *) '(a40,', n_params, 'a40)'
    ! len(param_names(1))=256 but only 39 characters taken here
    ! write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii))), ii=1, n_params)
    write(uopti, formheader) 'OF', (trim(adjustl(param_names(ii)(1 : 39))), ii = 1, n_params)
 
    ! output
    write(formparams, *) '( es40.14, ', n_params, '(es40.14) )'
    write(uopti, formparams) best_of, (best_paramset(ii), ii = 1, n_params)
 
    ! close file
    close(uopti)
 
    ! screen output
    call message()
    call message(' Optimized parameters written to ', trim(fname))
 
  end subroutine mrm_write_optifile
 
  ! ------------------------------------------------------------------
 
  !    NAME
  !        mrm_write_optinamelist
 
  !    PURPOSE
  !>       \brief Write final, optimized parameter set in a namelist format.
 
  !>       \details Write final, optimized parameter set in a namelist format.
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(:, :) :: parameters"                               (min, max, opti)
  !>       \param[in] "logical, dimension(size(parameters, 1)) :: maskpara"                   .true. if parameter was
  !>       calibrated
  !>       \param[in] "character(len = *), dimension(size(parameters, 1)) :: parameters_name" clear names of parameters
 
  !    HISTORY
  !>       \authors Juliane Mai
 
  !>       \date Dec 2013
 
  ! Modifications:
  ! Stephan Thober Oct 2015 - adapted to mRM
  ! Stephan Thober Nov 2016 - adapt header to routing process
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine mrm_write_optinamelist(parameters, maskpara, parameters_name)
 
    use mo_common_mhm_mrm_file, only : file_opti_nml, uopti_nml
    use mo_common_variables, only : dirconfigout, processmatrix
    use mo_string_utils, only : num2str
 
    implicit none
 
    ! (min, max, opti)
    real(dp), dimension(:, :), intent(in) :: parameters
 
    ! .true. if parameter was calibrated
    logical, dimension(size(parameters, 1)), intent(in) :: maskpara
 
    ! clear names of parameters
    character(len = *), dimension(size(parameters, 1)), intent(in) :: parameters_name
 
    character(256) :: fname
 
    character(3) :: flag
 
    integer(i4) :: err
 
    integer(i4) :: ipar
 
 
    ! open file
    fname = trim(adjustl(dirconfigout)) // trim(adjustl(file_opti_nml))
    open(uopti_nml, file = fname, status = 'unknown', action = 'write', iostat = err)
    if(err .ne. 0) then
      call message ('  IOError while openening "', trim(fname), '" Error-Code ', num2str(err))
    end if
 
    write(uopti_nml, *) '!global_parameters'
    write(uopti_nml, *) '!PARAMETER                       lower_bound  upper_bound          value   FLAG  SCALING'
 
    write(uopti_nml, *) '! ', trim(adjustl('routing'))
 
    if (processmatrix(8, 1) .eq. 1_i4) write(uopti_nml, *) '&routing1'
    if (processmatrix(8, 1) .eq. 2_i4) write(uopti_nml, *) '&routing2'
    if (processmatrix(8, 1) .eq. 3_i4) write(uopti_nml, *) '&routing3'
 
    do ipar = 1, size(parameters, 1)
      if (maskpara(ipar)) then
        flag = ' 1 '
      else
        flag = ' 0 '
      end if
      write(uopti_nml, *) trim(adjustl(parameters_name(ipar))), ' = ', &
              parameters(ipar, 1), ' , ', &
              parameters(ipar, 2), ' , ', &
              parameters(ipar, 3), ' , ', &
              flag, ', 1 '
    end do
 
    write(uopti_nml, *) '/'
    write(uopti_nml, *) ' '
 
    ! close file
    close(uopti_nml)
 
    ! screen output
    call message()
    call message(' Optimized parameters written in namelist format to ', trim(fname))
 
  end subroutine mrm_write_optinamelist
 
 
end module mo_mrm_write