LCOV - code coverage report
Current view: top level - mRM - mo_mrm_write.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 109 405 26.9 %
Date: 2023-08-24 16:49:59 Functions: 2 6 33.3 %

          Line data    Source code
       1             : !> \file mo_mrm_write.f90
       2             : !> \brief   \copybrief mo_mrm_write
       3             : !> \details \copydetails mo_mrm_write
       4             : 
       5             : !> \brief write of discharge and restart files
       6             : !> \details This module contains the subroutines for writing the discharge files and optionally the restart files.
       7             : !> \authors Stephan Thober
       8             : !> \date Aug 2015
       9             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      10             : !! mHM is released under the LGPLv3+ license \license_note
      11             : !> \ingroup f_mrm
      12             : module mo_mrm_write
      13             : 
      14             :   use mo_kind, only : i4, dp
      15             :   use mo_mrm_write_fluxes_states, only : OutputDataset
      16             :   use mo_mrm_global_variables, only : output_time_reference_mrm
      17             :   use mo_message, only : message, error_message
      18             : 
      19             :   implicit none
      20             : 
      21             :   public :: mrm_write
      22             :   public :: mrm_write_optinamelist
      23             :   public :: mrm_write_optifile
      24             : 
      25             :   private
      26             : 
      27             : contains
      28             : 
      29             :   ! ------------------------------------------------------------------
      30             : 
      31             :   !    NAME
      32             :   !        mrm_write
      33             : 
      34             :   !    PURPOSE
      35             :   !>       \brief write discharge and restart files
      36             : 
      37             :   !>       \details First, this subroutine calls the writing or restart files that only
      38             :   !>       succeeds if it happens after the write of mHM restart files because
      39             :   !>       mHM restart files must exist. Second, simulated discharge is aggregated to the daily
      40             :   !>       scale and then written to file jointly with observed discharge
      41             : 
      42             :   !    HISTORY
      43             :   !>       \authors Juliane Mai, Rohini Kumar & Stephan Thober
      44             : 
      45             :   !>       \date Aug 2015
      46             : 
      47             :   ! Modifications:
      48             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
      49             :   ! Pallav Shrestha, Husain Najafi Mar 2022 - refactoring for measured timestep output
      50             : 
      51          13 :   subroutine mrm_write
      52             : 
      53             :     use mo_common_constants, only : nodata_dp
      54             :     use mo_common_mhm_mrm_variables, only : evalPer, mrm_coupling_mode, nTstepDay, simPer, warmingDays
      55             :     use mo_common_variables, only : mrmFileRestartOut, domainMeta, write_restart
      56             :     use mo_mrm_global_variables, only : domain_mrm, &
      57             :                                         gauge, mRM_runoff, nGaugesTotal, nMeasPerDay
      58             :     use mo_mrm_restart, only : mrm_write_restart
      59             : 
      60             :     implicit none
      61             : 
      62             :     integer(i4) :: domainID, iDomain
      63             : 
      64             :     integer(i4) :: iDay, iSubDay, iS, iE
      65             : 
      66             :     integer(i4) :: maxDailyTimeSteps, maxMeasTimeSteps
      67             : 
      68             :     integer(i4) :: tt
      69             : 
      70             :     integer(i4) :: gg
      71             : 
      72             :     integer(i4) :: nTimeSteps
      73             : 
      74          13 :     real(dp), dimension(:, :), allocatable :: d_Qmod, subd_Qmod ! Sim discharge at daily and subdaily time step
      75          13 :     real(dp), dimension(:, :), allocatable :: d_Qobs, subd_Qobs ! Obs discharge at daily and subdaily time step
      76             : 
      77             :     ! between simulated and measured time scale
      78             :     integer(i4) :: factor
      79             : 
      80             :     ! simulated Timesteps per Day
      81             :     integer(i4) :: TPD_sim
      82             : 
      83             :     ! observed Timesteps per Day
      84             :     integer(i4) :: TPD_obs
      85             : 
      86             :     ! --------------------------------------------------------------------------
      87             :     ! WRITE CONFIG FILE
      88             :     ! --------------------------------------------------------------------------
      89           0 :     if (mrm_coupling_mode .eq. 0_i4) call write_configfile()
      90             : 
      91             :     ! --------------------------------------------------------------------------
      92             :     ! WRITE RESTART
      93             :     ! --------------------------------------------------------------------------
      94          13 :     if (write_restart) then
      95          38 :       do iDomain = 1, domainMeta%nDomains
      96          25 :         domainID = domainMeta%indices(iDomain)
      97          38 :         if (domainMeta%doRouting(iDomain)) call mrm_write_restart(iDomain, domainID, mrmFileRestartOut)
      98             :       end do
      99             :     end if
     100             : 
     101             :     ! --------------------------------------------------------------------------
     102             :     ! STORE DAILY DISCHARGE TIMESERIES OF EACH GAUGING STATION
     103             :     ! FOR SIMULATIONS DURING THE EVALUATION PERIOD
     104             :     !
     105             :     !  **** AT DAILY TIME STEPS ****
     106             :     ! Note:: Observed Q are stored only for the evaluation period and not for
     107             :     !        the warming days
     108             :     ! --------------------------------------------------------------------------
     109             : 
     110             :     ! copy time resolution to local variables
     111          13 :     TPD_sim = nTstepDay
     112          13 :     TPD_obs = nMeasPerDay
     113             : 
     114             :     ! check if modelled timestep is an integer multiple of measured timesteps
     115          13 :     if (modulo(TPD_sim, TPD_obs) .eq. 0) then
     116          13 :       factor = TPD_sim / TPD_obs
     117             :     else
     118           0 :       call error_message(' Error: Number of modelled datapoints is no multiple of measured datapoints per day')
     119             :     end if
     120             : 
     121             : 
     122          51 :     maxDailyTimeSteps = maxval(evalPer(1 : domainMeta%nDomains)%julEnd - evalPer(1 : domainMeta%nDomains)%julStart + 1)
     123          38 :     maxMeasTimeSteps  = maxval(evalPer(1 : domainMeta%nDomains)%julEnd - evalPer(1 : domainMeta%nDomains)%julStart + 1) * TPD_obs
     124          52 :     allocate(d_Qmod     (maxDailyTimeSteps, nGaugesTotal))
     125          39 :     allocate(d_Qobs     (maxDailyTimeSteps, nGaugesTotal))
     126          52 :     allocate(subd_Qmod  ( maxMeasTimeSteps, nGaugesTotal))
     127          39 :     allocate(subd_Qobs  ( maxMeasTimeSteps, nGaugesTotal))
     128       12101 :     d_Qmod     = nodata_dp
     129       12101 :     d_Qobs     = nodata_dp
     130       12101 :     subd_Qmod  = nodata_dp
     131       12101 :     subd_Qobs  = nodata_dp
     132             : 
     133             : 
     134             :     ! loop over domains
     135          38 :     do iDomain = 1, domainMeta%nDomains
     136          38 :       if (domainMeta%doRouting(iDomain)) then
     137          20 :         domainID = domainMeta%indices(iDomain)
     138             : 
     139             :         ! Convert simulated values to daily
     140          20 :         nTimeSteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nTstepDay
     141          20 :         iDay  = 0
     142        8786 :         do tt = warmingDays(iDomain) * nTstepDay + 1, nTimeSteps, nTstepDay
     143        8766 :           iS = tt
     144        8766 :           iE = tt + nTstepDay - 1
     145        8766 :           iDay = iDay + 1
     146             :           ! over gauges
     147       18282 :           do gg = 1, domain_mrm(iDomain)%nGauges
     148             :             ! simulation
     149       18992 :             d_Qmod(iDay, domain_mrm(iDomain)%gaugeIndexList(gg)) = &
     150      265158 :                     sum(mRM_runoff(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg))) / real(nTstepDay, dp)
     151             :           end do
     152             :         end do
     153             : 
     154          20 :         dailycheck: if (nMeasPerDay > 1) then
     155             :           ! Convert observed values to daily
     156           0 :           nTimeSteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nMeasPerDay
     157           0 :           iDay  = 0
     158           0 :           do tt = 1, nTimeSteps, nMeasPerDay
     159           0 :             iS = tt
     160           0 :             iE = tt + nMeasPerDay - 1
     161           0 :             iDay = iDay + 1
     162             :             ! over gauges
     163           0 :             do gg = 1, domain_mrm(iDomain)%nGauges
     164             :               ! when -9999 value/s are present in current day, daily value remain -9999.
     165           0 :               if (.not.(any(gauge%Q(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg)) == nodata_dp))) then
     166             :                 ! observation
     167           0 :                 d_Qobs(iDay, domain_mrm(iDomain)%gaugeIndexList(gg)) = &
     168           0 :                         sum( gauge%Q(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg))) / real(nMeasPerDay, dp)
     169             :               end if
     170             :             end do
     171             :           end do
     172             :         else
     173             :           ! observed values are already at daily (nMeasPerDay = 1) and stored for evalper
     174             :           ! over gauges
     175          42 :           do gg = 1, domain_mrm(iDomain)%nGauges
     176             :             ! observation
     177       11363 :             d_Qobs(:, domain_mrm(iDomain)%gaugeIndexList(gg)) = gauge%Q(:, domain_mrm(iDomain)%gaugeIndexList(gg))
     178             :           end do
     179             :         end if dailycheck
     180             : 
     181             : 
     182          20 :         subdailycheck: if (nMeasPerDay > 1) then
     183             : 
     184             :           ! Convert simulated values to subdaily
     185           0 :           nTimeSteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nTstepDay
     186           0 :           iSubDay = 0
     187           0 :           do tt = warmingDays(iDomain) * nTstepDay + 1, nTimeSteps, factor
     188           0 :             iS = tt
     189           0 :             iE = tt + factor - 1
     190           0 :             iSubDay = iSubDay + 1
     191             :             ! over gauges
     192           0 :             do gg = 1, domain_mrm(iDomain)%nGauges
     193             :               ! simulation
     194           0 :               subd_Qmod(iSubDay, domain_mrm(iDomain)%gaugeIndexList(gg)) = &
     195           0 :                       sum(mRM_runoff(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg))) / real(factor, dp)
     196             :             end do
     197             :           end do
     198             : 
     199             :           ! Convert observed values to subdaily
     200             : 
     201             :           ! observed values are already at subdaily (nMeasPerDay) and stored for evalper
     202             :           ! over gauges
     203           0 :           do gg = 1, domain_mrm(iDomain)%nGauges
     204             :             ! observation
     205           0 :             subd_Qobs(:, domain_mrm(iDomain)%gaugeIndexList(gg)) = gauge%Q(:, domain_mrm(iDomain)%gaugeIndexList(gg))
     206             :           end do
     207             : 
     208             :         end if subdailycheck
     209             : 
     210             :       end if
     211             :     end do
     212             : 
     213             : 
     214             :     ! write in an ASCII file          ! OBS[nModeling_days X nGauges_total] , SIM[nModeling_days X nGauges_total]
     215             :     ! ToDo: is this if statement reasonable
     216          13 :     if (allocated(gauge%Q)) call write_daily_obs_sim_discharge(d_Qobs(:, :), d_Qmod(:, :))
     217             : 
     218             :     ! write in an ASCII file          ! OBS[nMeasPerDay X nGauges_total] , SIM[nMeasPerDay X nGauges_total]
     219          13 :     if (nMeasPerDay > 1 .and. allocated(gauge%Q)) call write_subdaily_obs_sim_discharge(subd_Qobs(:, :), subd_Qmod(:, :), factor)
     220             :     ! The subdaily routine is only called if subdaily Q data is provided
     221             : 
     222             :     ! free space
     223          13 :     deallocate(d_Qmod, d_Qobs, subd_Qmod, subd_Qobs)
     224          13 :   end subroutine mrm_write
     225             :   ! ------------------------------------------------------------------
     226             : 
     227             :   !    NAME
     228             :   !        write_configfile
     229             : 
     230             :   !    PURPOSE
     231             :   !>       \brief This modules writes the results of the configuration into an ASCII-file
     232             : 
     233             :   !>       \details TODO: add description
     234             : 
     235             :   !    HISTORY
     236             :   !>       \authors Christoph Schneider
     237             : 
     238             :   !>       \date May 2013
     239             : 
     240             :   ! Modifications:
     241             :   ! Juliane Mai    May 2013 - module version and documentation
     242             :   ! Stephan Thober Jun 2014 - bug fix in L11 config print out
     243             :   ! Stephan Thober Jun 2014 - updated read_restart
     244             :   ! Rohini, Luis   Jul 2015 - updated version, L1 level prints
     245             :   ! Stephan Thober Sep 2015 - updated write of stream network
     246             :   ! Stephan Thober Nov 2016 - adapted write to selected case for routing process
     247             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     248             : 
     249           0 :   Subroutine write_configfile
     250             : 
     251          13 :     use mo_common_constants, only : nodata_dp
     252             :     use mo_common_file, only : file_config, uconfig
     253             :     use mo_common_mHM_mRM_variables, only : LCyearId, SimPer, evalPer, mrm_coupling_mode, read_restart, &
     254             :                                             resolutionRouting, timeStep, warmPer
     255             :     use mo_common_variables, only : LC_year_end, LC_year_start, LCfilename, &
     256             :                                     dirConfigOut, dirLCover, dirMorpho, dirOut, mrmFileRestartOut, global_parameters, &
     257             :                                     global_parameters_name, level0, level1, domainMeta, nLCoverScene, processMatrix, &
     258             :                                     resolutionHydrology, write_restart
     259             :     use mo_kind, only : dp, i4
     260             :     use mo_mrm_file, only : version
     261             :     use mo_mrm_global_variables, only : InflowGauge, L11_L1_Id, L11_fromN, L11_label, &
     262             :                                         L11_length, L11_netPerm, L11_rOrder, L11_slope, L11_toN, L1_L11_Id, domain_mrm, &
     263             :                                         dirGauges, dirTotalRunoff, gauge, level11, nGaugesTotal, nInflowGaugesTotal
     264             :     use mo_string_utils, only : num2str
     265             :     use mo_utils, only : ge
     266             :     use mo_os, only : check_path_isdir
     267             : 
     268             :     implicit none
     269             : 
     270             :     character(256) :: fName
     271             : 
     272             :     integer(i4) :: i, iDomain, domainID, j
     273             : 
     274             :     integer(i4) :: err
     275             : 
     276             : 
     277           0 :     fName = trim(adjustl(dirConfigOut)) // trim(adjustl(file_config))
     278           0 :     call message()
     279           0 :     call message('  Log-file written to ', trim(fName))
     280             :     !checking whether the directory exists where the file shall be created or opened
     281           0 :     call check_path_isdir(trim(adjustl(dirConfigOut)), raise=.true.)
     282           0 :     open(uconfig, file = fName, status = 'unknown', action = 'write', iostat = err)
     283           0 :     if (err .ne. 0) then
     284           0 :       call error_message('  Problems while creating File. Error-Code ', num2str(err))
     285             :     end if
     286           0 :     write(uconfig, 200)
     287           0 :     write(uconfig, 100) 'mRM-UFZ v-' // trim(version)
     288           0 :     write(uconfig, 100) 'S. Thober, L. Samaniego & R. Kumar, UFZ'
     289           0 :     write(uconfig, 200)
     290           0 :     write(uconfig, 100)
     291           0 :     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         '
     292           0 :     write(uconfig, 100)
     293           0 :     write(uconfig, 103) 'Number of domains            ', domainMeta%nDomains
     294           0 :     write(uconfig, 103) 'Total No. of gauges         ', nGaugesTotal
     295           0 :     write(uconfig, 103)    'Time Step [h]               ', timeStep
     296           0 :     do iDomain = 1, domainMeta%nDomains
     297           0 :       domainID = domainMeta%indices(iDomain)
     298           0 :       write(uconfig, 103) 'Total No. of nodes          ', level11(iDomain)%nCells
     299           0 :       write(uconfig, 103) 'No. of cells L0             ', level0(domainMeta%L0DataFrom(iDomain))%nCells
     300           0 :       write(uconfig, 103) 'No. of cells L1             ', level1(iDomain)%nCells
     301           0 :       write(uconfig, 103) 'No. of cells L11            ', level11(iDomain)%nCells
     302             : 
     303             :       !    select case (iFlag_cordinate_sys)
     304             :       !    case (0)
     305           0 :       write(uconfig, 301)      'domain  ', domainID, '   Hydrology Resolution [m]      ', resolutionHydrology(iDomain)
     306           0 :       write(uconfig, 301)   'domain  ', domainID, '   Routing Resolution [m]        ', resolutionRouting(iDomain)
     307             :       !    case(1)
     308             :       !      write(uconfig, 302)       'domain  ',domainID, '   Hydrology Resolution [o]      ', resolutionHydrology(iDomain)
     309             :       !      write(uconfig, 302)   'domain  ',domainID, '   Routing Resolution [o]        ', resolutionRouting(iDomain)
     310             :       !    end select
     311             :     end do
     312           0 :     write(uconfig, 126)    'Flag READ  restart            ', read_restart
     313           0 :     write(uconfig, 126)    'Flag WRITE restart            ', write_restart
     314             :     !
     315             :     !******************
     316             :     ! Model Run period
     317             :     !******************
     318           0 :     do iDomain = 1, domainMeta%nDomains
     319           0 :       domainID = domainMeta%indices(iDomain)
     320           0 :       write(uconfig, 115) '                      Model Run Periods for domain ', num2str(domainID)
     321             :       write(uconfig, 116) &
     322           0 :               'From                To', &
     323           0 :               '   Day Month  Year   Day Month  Year'
     324             :       write(uconfig, 117)  &
     325           0 :               'Warming Period (1)            ', &
     326           0 :               warmPer(iDomain)%dStart, warmPer(iDomain)%mStart, warmPer(iDomain)%yStart, &
     327           0 :               warmPer(iDomain)%dEnd, warmPer(iDomain)%mEnd, warmPer(iDomain)%yEnd
     328             :       write(uconfig, 117)  &
     329           0 :               'Evaluation Period (2)         ', &
     330           0 :               evalPer(iDomain)%dStart, evalPer(iDomain)%mStart, evalPer(iDomain)%yStart, &
     331           0 :               evalPer(iDomain)%dEnd, evalPer(iDomain)%mEnd, evalPer(iDomain)%yEnd
     332             :       write(uconfig, 117)  &
     333           0 :               'Simulation Period (1)+(2)     ', &
     334           0 :               SimPer(iDomain)%dStart, SimPer(iDomain)%mStart, SimPer(iDomain)%yStart, &
     335           0 :               SimPer(iDomain)%dEnd, SimPer(iDomain)%mEnd, SimPer(iDomain)%yEnd
     336             :     end do
     337             : 
     338             :     !*********************************
     339             :     ! Model Land Cover Observations
     340             :     !*********************************
     341           0 :     if (processMatrix(8, 1) .eq. 1) then
     342           0 :       do iDomain = 1, domainMeta%nDomains
     343           0 :         domainID = domainMeta%indices(iDomain)
     344           0 :         write(uconfig, 118) '       Land Cover Observations for domain ', num2str(domainID)
     345           0 :         write(uconfig, 119) ' Start Year', ' End Year', '    Land cover scene', 'Land Cover File'
     346           0 :         do i = 1, nLCoverScene
     347           0 :           write(uconfig, 120) LC_year_start(i), LC_year_end(i), &
     348           0 :                   LCyearId(max(evalPer(iDomain)%yStart, LC_year_start(i)), iDomain), trim(LCfilename(i))
     349             :         end do
     350             :       end do
     351             :     end if
     352             :     !*********************************
     353             :     ! Initial Parameter Ranges
     354             :     !*********************************
     355           0 :     write(uconfig, 121) '  Initial Transfer Function Parameter Ranges (gammas)  '
     356             :     !
     357             :     ! Transfer functions
     358             :     write(uconfig, 122)      &
     359           0 :             '         i', '            min', '            max', '        current', &
     360           0 :             '                               name'
     361           0 :     do i = 1, size(global_parameters, 1)
     362             :       write(uconfig, 123) &
     363           0 :               i, global_parameters(i, 1), global_parameters(i, 2), global_parameters(i, 3), &
     364           0 :               trim(adjustl(global_parameters_name(i)))
     365             :     end do
     366             :     ! domain runoff data
     367           0 :     write(uconfig, 202) '                domain Runoff Data                '
     368           0 :     write(uconfig, 107) ' Gauge No.', '  domain Id', '     Qmax[m3/s]', '     Qmin[m3/s]'
     369           0 :     do i = 1, nGaugesTotal
     370           0 :       if(any(gauge%Q(:, i) > nodata_dp)) then
     371           0 :         write(uconfig, 108) i, gauge%domainId(i), maxval(gauge%Q(:, i), gauge%Q(:, i) > nodata_dp), &
     372           0 :                 minval(gauge%Q(:, i), gauge%Q(:, i) > nodata_dp)
     373             :       else
     374           0 :         write(uconfig, 108) i, gauge%domainId(i), nodata_dp, nodata_dp
     375             :       end if
     376             :     end do
     377             :     ! inflow gauge data
     378           0 :     if (nInflowGaugesTotal .GT. 0) then
     379           0 :       write(uconfig, 202) '                domain Inflow Data                 '
     380           0 :       write(uconfig, 107) ' Gauge No.', '  domain Id', '     Qmax[m3/s]', '     Qmin[m3/s]'
     381           0 :       do i = 1, nInflowGaugesTotal
     382           0 :         if(all(InflowGauge%Q(:, i) > nodata_dp)) then
     383           0 :           write(uconfig, 108) i, InflowGauge%domainId(i), maxval(InflowGauge%Q(:, i), InflowGauge%Q(:, i) > nodata_dp), &
     384           0 :                   minval(InflowGauge%Q(:, i), InflowGauge%Q(:, i) > nodata_dp)
     385             :         else
     386           0 :           write(uconfig, 108) i, InflowGauge%domainId(i), nodata_dp, nodata_dp
     387             :         end if
     388             :       end do
     389             :     end if
     390             :     ! domain config
     391           0 :     write(uconfig, 218) 'domain-wise Configuration'
     392           0 :     do iDomain = 1, domainMeta%nDomains
     393           0 :       domainID = domainMeta%indices(iDomain)
     394           0 :       write(uconfig, 103) 'domain No.                   ', domainID, &
     395           0 :               'No. of gauges               ', domain_mrm(iDomain)%nGauges
     396             : 
     397           0 :       write(uconfig, 222)   'Directory list'
     398             : 
     399           0 :       write(uconfig, 224) 'Directory to morphological input         ', dirMorpho(iDomain)
     400           0 :       write(uconfig, 224) 'Directory to land cover input            ', dirLCover(iDomain)
     401           0 :       write(uconfig, 224) 'Directory to gauging station input       ', dirGauges(iDomain)
     402           0 :       if (mrm_coupling_mode .eq. 0) then
     403           0 :         write(uconfig, 224) 'Directory to simulated runoff input      ', dirTotalRunoff(iDomain)
     404             :       end if
     405           0 :       write(uconfig, 224) 'Directory to write output by default     ', dirOut(iDomain)
     406           0 :       write(uconfig, 224) 'File to write mRM output when restarted  ', mrmFileRestartOut(iDomain)
     407             : 
     408           0 :       write(uconfig, 102) 'River Network  (Routing level)'
     409           0 :       write(uconfig, 100) 'Label 0 = intermediate draining cell '
     410           0 :       write(uconfig, 100) 'Label 1 = headwater cell             '
     411           0 :       write(uconfig, 100) 'Label 2 = sink cell                  '
     412             : 
     413           0 :       if (processMatrix(8, 1) .eq. 1_i4) then
     414           0 :         write(uconfig, 104) '   Overall', &
     415           0 :                 '      From', &
     416           0 :                 '        To', &
     417           0 :                 '   Routing', &
     418           0 :                 '     Label', &
     419           0 :                 '    Length', &
     420           0 :                 '      Mean', &
     421           0 :                 '      Link', &
     422           0 :                 '   Routing', &
     423           0 :                 '   Routing', &
     424           0 :                 '  Sequence', &
     425           0 :                 '          ', &
     426           0 :                 '          ', &
     427           0 :                 '     Slope'
     428             :         !
     429           0 :         write(uconfig, 105) '        Id', &
     430           0 :                 '      Node', &
     431           0 :                 '      Node', &
     432           0 :                 '', &
     433           0 :                 '', &
     434           0 :                 '      [km]', &
     435           0 :                 '    [o/oo]'
     436             :         !
     437           0 :         do j = level11(iDomain)%iStart, level11(iDomain)%iEnd - 1
     438           0 :           i = L11_netPerm(j) + level11(iDomain)%iStart - 1 ! adjust permutation for multi-domain option
     439           0 :           write(uconfig, 106) i, L11_fromN(i), L11_toN(i), L11_rOrder(i), L11_label(i), &
     440           0 :                   L11_length(i) / 1000.0_dp, L11_slope(i) * 1.0e3_dp
     441             :         end do
     442             : 
     443           0 :       else if (processMatrix(8, 1) .eq. 2_i4) then
     444           0 :         write(uconfig, 134) '   Overall', &
     445           0 :                 '      From', &
     446           0 :                 '        To', &
     447           0 :                 '   Routing', &
     448           0 :                 '     Label', &
     449           0 :                 '      Link', &
     450           0 :                 '   Routing', &
     451           0 :                 '   Routing', &
     452           0 :                 '  Sequence', &
     453           0 :                 '          '
     454             :         !
     455           0 :         write(uconfig, 135) '        Id', &
     456           0 :                 '      Node', &
     457           0 :                 '      Node', &
     458           0 :                 '', &
     459           0 :                 ''
     460             :         !
     461           0 :         do j = level11(iDomain)%iStart, level11(iDomain)%iEnd - 1
     462           0 :           i = L11_netPerm(j) + level11(iDomain)%iStart - 1 ! adjust permutation for multi-domain option
     463           0 :           write(uconfig, 136) i, L11_fromN(i), L11_toN(i), L11_rOrder(i), L11_label(i)
     464             :         end do
     465             :       end if
     466             :       ! draining node at L11
     467           0 :       write(uconfig, 109)  '   Overall', '     domain', &
     468           0 :               '      Cell', '   Routing', &
     469           0 :               '        Id', '   Node Id'
     470           0 :       do i = level11(iDomain)%Id(1), level11(iDomain)%Id(level11(iDomain)%nCells)
     471           0 :         write(uconfig, 110) i + level11(iDomain)%iStart - 1, i
     472             :       end do
     473             : 
     474             :       ! L1 level information
     475           0 :       write(uconfig, 111)  '  Modeling', '   Routing', ' Effective', &
     476           0 :               '      Cell', '   Cell Id', '      Area', &
     477           0 :               '        Id', '       [-]', '     [km2]'
     478           0 :       if (ge(resolutionRouting(iDomain), resolutionHydrology(iDomain))) then
     479           0 :         do i = level1(iDomain)%Id(1), level1(iDomain)%Id(level1(iDomain)%nCells)
     480           0 :           write(uconfig, 113) i + level1(iDomain)%iStart - 1, L1_L11_Id (i + level1(iDomain)%iStart - 1), &
     481           0 :             level1(iDomain)%CellArea(i) * 1.E-6_dp
     482           0 :       domainID = domainMeta%indices(iDomain)
     483             :         end do
     484             :       else
     485           0 :         do i = level11(iDomain)%Id(1), level11(iDomain)%Id(level11(iDomain)%nCells)
     486           0 :           write(uconfig, 110) i + level11(iDomain)%iStart - 1, L11_L1_Id (i + level11(iDomain)%iStart - 1)
     487             :         end do
     488             :       end if
     489           0 :       write(uconfig, 114)  ' Total[km2]', sum(level1(iDomain)%CellArea) * 1.E-6_dp
     490             :     end do
     491             : 
     492           0 :     write(uconfig, *)
     493           0 :     close(uconfig)
     494             : 
     495             :     !! Formats
     496             :     100 format (a80)
     497             :     102 format (/ 30('-') / a30 / 30('-'))
     498             :     103 format (a20, 10x, i10)
     499             :     104 format (/ 75('-') / 5a10, 5x, 2a10 / 5a10, 5x, 2a10)
     500             :     105 format (5a10, 5x, 2a10 / 75('-'))
     501             :     106 format (5i10, 5x, 2f10.3)
     502             :     107 format (2a10, 2a15)
     503             :     108 format (2i10, 2f15.3)
     504             :     !
     505             :     109 format (/ 20('-') / 2a10 / 2a10 / 2a10 / 20('-'))
     506             :     110 format (2i10)
     507             :     !
     508             :     111 format (/ 30('-') / 3a10 / 3a10 / 3a10 /  30('-'))
     509             :     113 format (            2i10,   1f10.3         )
     510             :     114 format (30('-') / a15, 5x, 1f10.3 /)
     511             :     !
     512             :     115 format (/61('-')/ a50, a10 /61('-'))
     513             :     116 format (39x, a22 / 25x, a36)
     514             :     117 format (3(a25, 6(i6)))
     515             :     !
     516             :     118 format (/50('-')/ a40, a10  /50('-'))
     517             :     119 format (a10, a10, a20, a20/)
     518             :     120 format (i10, i10, 10x, i10, a20)
     519             :     !
     520             :     121 format (/55('-')/ a55 /55('-'))
     521             :     122 format (a10, 3a15, a35)
     522             :     123 format (i10, 3f15.3, a35)
     523             :     !
     524             :     126 format (a30, 9x, L1)
     525             :     !
     526             :     134 format (/ 50('-') / 5a10 / 5a10)
     527             :     135 format (5a10 / 50('-'))
     528             :     136 format (5i10)
     529             :     !
     530             :     200 format (80('-'))
     531             :     201 format (a80)
     532             :     202 format (/50('-')/ a50 /50('-'))
     533             :     !
     534             :     218 format (/ 80('-')/ 26x, a24, 26x, /80('-'))
     535             :     222 format (/80('-')/ 26x, a21 /80('-'))
     536             :     224 format (a40, 5x, a256)
     537             : 
     538             :     301 format (a7, i2, a32, f15.0)
     539             :     ! 302 format (a7, i2, a32,es20.8)
     540           0 :   end Subroutine write_configfile
     541             : 
     542             :   ! ------------------------------------------------------------------
     543             : 
     544             :   !    NAME
     545             :   !        write_daily_obs_sim_discharge
     546             : 
     547             :   !    PURPOSE
     548             :   !>       \brief Write a file for the daily observed and simulated discharge timeseries
     549             :   !>       during the evaluation period for each gauging station
     550             : 
     551             :   !>       \details Write a file for the daily observed and simulated discharge timeseries
     552             :   !>       during the evaluation period for each gauging station
     553             : 
     554             :   !    INTENT(IN)
     555             :   !>       \param[in] "real(dp), dimension(:, :) :: Qobs" daily time series of observed dischargedims = (nModeling_days
     556             :   !>       , nGauges_total)
     557             :   !>       \param[in] "real(dp), dimension(:, :) :: Qsim" daily time series of modeled dischargedims = (nModeling_days ,
     558             :   !>       nGauges_total)
     559             : 
     560             :   !    HISTORY
     561             :   !>       \authors Rohini Kumar
     562             : 
     563             :   !>       \date August 2013
     564             : 
     565             :   ! Modifications:
     566             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     567             : 
     568          26 :   subroutine write_daily_obs_sim_discharge(Qobs, Qsim)
     569             : 
     570           0 :     use mo_common_constants, only : nodata_dp
     571             :     use mo_common_mhm_mrm_variables, only : evalPer
     572             :     use mo_common_variables, only : dirOut, domainMeta
     573             :     use mo_errormeasures, only : kge, nse
     574             :     use mo_julian, only : dec2date
     575             :     use mo_mrm_file, only : file_daily_discharge, ncfile_discharge, udaily_discharge
     576             :     use mo_mrm_global_variables, only : domain_mrm, gauge
     577             :     use mo_string_utils, only : num2str
     578             :     use mo_utils, only : ge
     579             :     use mo_netcdf, only : NcDataset, NcDimension, NcVariable
     580             : 
     581             :     implicit none
     582             : 
     583             :     ! daily time series of observed dischargedims = (nModeling_days , nGauges_total)
     584             :     real(dp), dimension(:, :), intent(in) :: Qobs
     585             : 
     586             :     ! daily time series of modeled dischargedims = (nModeling_days , nGauges_total)
     587             :     real(dp), dimension(:, :), intent(in) :: Qsim
     588             : 
     589             :     character(256) :: fName, formHeader, formData, dummy
     590             : 
     591             :     integer(i4) :: domainID, iDomain, gg, tt, err
     592             : 
     593             :     integer(i4) :: igauge_start, igauge_end
     594             : 
     595             :     integer(i4) :: day, month, year
     596             : 
     597             :     integer(i4) :: tlength
     598             : 
     599             :     ! time axis
     600          13 :     integer(i4), allocatable, dimension(:) :: taxis
     601             : 
     602          13 :     real(dp) :: newTime
     603             : 
     604             :     ! nc related variables
     605             :     type(NcDataset) :: nc_out
     606             :     type(NcDimension) :: dim, dim_bnd
     607             :     type(NcVariable) :: var
     608             : 
     609             :     ! initalize igauge_start
     610          13 :     igauge_start = 1
     611             : 
     612             :     ! domain loop
     613          38 :     do iDomain = 1, domainMeta%nDomains
     614          25 :       domainID = domainMeta%indices(iDomain)
     615          25 :       if(domain_mrm(iDomain)%nGauges .lt. 1) cycle
     616             : 
     617             :       ! estimate igauge_end
     618          19 :       igauge_end = igauge_start + domain_mrm(iDomain)%nGauges - 1
     619             : 
     620             :       ! check the existance of file
     621          19 :       fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(file_daily_discharge))
     622          19 :       open(udaily_discharge, file = trim(fName), status = 'unknown', action = 'write', iostat = err)
     623          19 :       if(err .ne. 0) then
     624           0 :         call error_message('  IOError while openening "', trim(fName), '". Error-Code ', num2str(err))
     625             :       end if
     626             : 
     627             :       ! header
     628          19 :       write(formHeader, *) '( 4a8, ', domain_mrm(iDomain)%nGauges, '(2X, a5, i10.10, 2X, a5, i10.10) )'
     629          19 :       write(udaily_discharge, formHeader) 'No', 'Day', 'Mon', 'Year', &
     630          22 :               ('Qobs_', gauge%gaugeId(gg), &
     631          60 :                       'Qsim_', gauge%gaugeId(gg), gg = igauge_start, igauge_end)
     632             : 
     633             :       ! form data
     634          19 :       write(formData, *) '( 4I8, ', domain_mrm(iDomain)%nGauges, '(2X,   f15.7 , 2X,  f15.7  ) )'
     635             : 
     636             :       ! write data
     637          19 :       newTime = real(evalPer(iDomain)%julStart, dp) - 0.5_dp
     638             : 
     639        8420 :       do tt = 1, (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1)
     640        8401 :         call dec2date(newTime, yy = year, mm = month, dd = day)
     641       17897 :         write(udaily_discharge, formData) tt, day, month, year, (Qobs(tt, gg), Qsim(tt, gg), gg = igauge_start, igauge_end)
     642        8420 :         newTime = newTime + 1.0_dp
     643             :       end do
     644             : 
     645             :       ! close file
     646          19 :       close(udaily_discharge)
     647             : 
     648             :       ! ======================================================================
     649             :       ! write netcdf file
     650             :       ! ======================================================================
     651          19 :       fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(ncfile_discharge))
     652          19 :       nc_out = NcDataset(trim(fName), "w")
     653          19 :       tlength = evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1
     654             :       ! write time
     655          57 :       allocate(taxis(tlength))
     656             : 
     657             :       ! tt is dependent on the unit of the time axis and is set to hours in mRM
     658             :       select case( output_time_reference_mrm)
     659             :         case(0)
     660        8420 :           forall(tt = 1 : tlength) taxis(tt) = (tt-1) * 24
     661             :         case(1)
     662           0 :           forall(tt = 1 : tlength) taxis(tt) = tt * 24 - 12
     663             :         case(2)
     664          19 :           forall(tt = 1 : tlength) taxis(tt) = tt * 24
     665             :       end select
     666             : 
     667          19 :       call dec2date(real(evalPer(iDomain)%julStart, dp) - 0.5_dp, yy = year, mm = month, dd = day)
     668          19 :       dim = nc_out%setDimension("time", tlength)
     669          38 :       var = nc_out%setVariable("time", "i32", [dim])
     670          19 :       call var%setData(taxis)
     671             :       call var%setAttribute( &
     672             :         "units", &
     673             :         'hours since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' 00:00:00' &
     674          19 :       )
     675          19 :       call var%setAttribute("long_name", "time in hours")
     676          19 :       call var%setAttribute("bounds", "time_bnds")
     677          19 :       call var%setAttribute("axis", "T")
     678          19 :       dim_bnd = nc_out%setDimension("bnds", 2)
     679          57 :       var = nc_out%setVariable("time_bnds", "i32", [dim_bnd, dim])
     680        8420 :       do tt = 1, tlength
     681       25203 :         call var%setData((tt - 1) * 24, (/1, tt/))
     682       25222 :         call var%setData(tt * 24, (/2, tt/))
     683             :       end do
     684          19 :       deallocate(taxis)
     685             :       ! write gauges
     686          41 :       do gg = igauge_start, igauge_end
     687          44 :         var = nc_out%setVariable('Qsim_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
     688          22 :         call var%setFillValue(nodata_dp)
     689          22 :         call var%setData(Qsim(1 : tlength, gg))
     690          22 :         call var%setAttribute("units", "m3 s-1")
     691          22 :         call var%setAttribute("long_name", 'simulated discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
     692          22 :         call var%setAttribute("missing_value", nodata_dp)
     693             :         ! write observed discharge at that gauge
     694          44 :         var = nc_out%setVariable('Qobs_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
     695          22 :         call var%setFillValue(nodata_dp)
     696          22 :         call var%setData(Qobs(1 : tlength, gg))
     697          22 :         call var%setAttribute("units", "m3 s-1")
     698          22 :         call var%setAttribute("long_name", 'observed discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
     699          41 :         call var%setAttribute("missing_value", nodata_dp)
     700             :       end do
     701          19 :       call nc_out%close()
     702             : 
     703             :       ! ======================================================================
     704             :       ! screen output
     705             :       ! ======================================================================
     706             : 
     707             :       ! if ( nMeasPerDay == 1_i4 ) then ! only print daily stats for daily Qobs
     708             : 
     709          19 :       call message()
     710          19 :       write(dummy, '(I3)') domainID
     711          19 :       call message('  OUTPUT: saved daily discharge file for domain ', trim(adjustl(dummy)))
     712          19 :       call message('    to ', trim(fname))
     713          41 :       do gg = igauge_start, igauge_end
     714       11362 :         if (count(ge(Qobs(:, gg), 0.0_dp)) > 1)  then
     715             :           call message('    KGE of daily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
     716       10611 :                   trim(adjustl(num2str(kge(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
     717             :           call message('    NSE of daily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
     718       10611 :                   trim(adjustl(num2str(nse(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
     719             :         end if
     720             :       end do
     721             : 
     722             :       ! end if
     723             : 
     724             :       ! update igauge_start
     725          19 :       igauge_start = igauge_end + 1
     726             :       !
     727             :     end do
     728             :     !
     729          13 :   end subroutine write_daily_obs_sim_discharge
     730             : 
     731             : 
     732             :   ! ------------------------------------------------------------------
     733             : 
     734             :   !    NAME
     735             :   !        write_subdaily_obs_sim_discharge
     736             : 
     737             :   !    PURPOSE
     738             :   !>       \brief Write a file for the simulated discharge timeseries
     739             :   !>       during the evaluation period for each gauging station
     740             : 
     741             :   !>       \details Write a file for the simulated discharge timeseries
     742             :   !>       during the evaluation period for each gauging station
     743             : 
     744             :   !    INTENT(IN)
     745             :   !>       \param[in] "real(dp), dimension(:, :) :: Qobs"  time series of observed discharge dims = (nMeasTimeSteps ,
     746             :   !>       nGauges_total)
     747             :   !>       \param[in] "real(dp), dimension(:, :) :: Qsim"  time series of modeled discharge dims = (nMeasTimeSteps ,
     748             :   !>       nGauges_total)
     749             :   !>       \param[in] "integer(i4),              :: factor" ratio of modelled time steps per day to observation time steps per day
     750             : 
     751             :   !    HISTORY
     752             :   !>       \authors Rohini Kumar
     753             : 
     754             :   !>       \date August 2013
     755             : 
     756             :   ! Modifications:
     757             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     758             :   ! Pallav Shrestha Jul 2021 - ported for printing out simulations at hourly time step
     759             :   ! Pallav Shrestha, Husain Najafi Mar 2022 - refactoring for subdaily timestep output
     760             : 
     761           0 :   subroutine write_subdaily_obs_sim_discharge(Qobs, Qsim, factor)
     762             : 
     763          13 :     use mo_common_constants, only : nodata_dp
     764             :     use mo_common_mhm_mrm_variables, only : evalPer
     765             :     use mo_common_variables, only : dirOut, domainMeta
     766             :     use mo_errormeasures, only : kge, nse
     767             :     use mo_julian, only : dec2date
     768             :     use mo_mrm_file, only : ncfile_subdaily_discharge, file_subdaily_discharge, &
     769             :                             usubdaily_discharge
     770             :     use mo_mrm_global_variables, only : domain_mrm, gauge, nMeasPerDay
     771             :     use mo_string_utils, only : num2str
     772             :     use mo_utils, only : ge
     773             :     use mo_netcdf, only : NcDataset, NcDimension, NcVariable
     774             : 
     775             :     implicit none
     776             : 
     777             :     ! time series of observed discharge. dims = (nMeasTimeSteps , nGauges_total)
     778             :     real(dp), dimension(:, :), intent(in) :: Qobs
     779             : 
     780             :     ! time series of modeled discharge. dims = (nMeasTimeSteps , nGauges_total)
     781             :     real(dp), dimension(:, :), intent(in) :: Qsim
     782             : 
     783             :     ! ratio of modelled time steps per day to observation time steps per day
     784             :     integer(i4),               intent(in) :: factor
     785             : 
     786             :     character(256) :: fName, formHeader, formData, dummy
     787             : 
     788             :     integer(i4) :: domainID, iDomain, gg, tt, err
     789             : 
     790             :     integer(i4) :: igauge_start, igauge_end
     791             : 
     792             :     integer(i4) :: hour, day, month, year
     793             : 
     794             :     integer(i4) :: tlength
     795             : 
     796             :     ! time axis
     797           0 :     integer(i4), allocatable, dimension(:) :: taxis
     798             : 
     799           0 :     real(dp) :: newTime
     800             : 
     801             :     ! nc related variables
     802             :     type(NcDataset) :: nc_out
     803             :     type(NcDimension) :: dim, dim_bnd
     804             :     type(NcVariable) :: var
     805             : 
     806             :     ! use minutes if needed
     807             :     logical :: use_minutes
     808             :     integer(i4) :: time_unit_factor
     809             : 
     810             :     ! initalize igauge_start
     811           0 :     igauge_start = 1
     812             : 
     813             :     ! domain loop
     814           0 :     do iDomain = 1, domainMeta%nDomains
     815           0 :       domainID = domainMeta%indices(iDomain)
     816           0 :       if(domain_mrm(iDomain)%nGauges .lt. 1) cycle
     817             : 
     818             :       ! estimate igauge_end
     819           0 :       igauge_end = igauge_start + domain_mrm(iDomain)%nGauges - 1
     820             : 
     821             : 
     822             :       ! ======================================================================
     823             :       ! write text file
     824             :       ! ======================================================================
     825             : 
     826             :       ! check the existance of file
     827           0 :       fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(file_subdaily_discharge))
     828           0 :       open(usubdaily_discharge, file = trim(fName), status = 'unknown', action = 'write', iostat = err)
     829           0 :       if(err .ne. 0) then
     830           0 :         call error_message('  IOError while openening "', trim(fName), '". Error-Code ', num2str(err))
     831             :       end if
     832             : 
     833             :       ! header
     834           0 :       write(formHeader, *) '( 5a8, ', domain_mrm(iDomain)%nGauges, '(2X, a5, i10.10, 2X, a5, i10.10) )'
     835           0 :       write(usubdaily_discharge, formHeader) 'No', 'Hour', 'Day', 'Mon', 'Year', &
     836           0 :               ('Qobs_', gauge%gaugeId(gg), &
     837           0 :                       'Qsim_', gauge%gaugeId(gg), gg = igauge_start, igauge_end)
     838             : 
     839             :       ! form data
     840           0 :       write(formData, *) '( 5I8, ', domain_mrm(iDomain)%nGauges, '(2X,   f15.7 , 2X,  f15.7  ) )'
     841             : 
     842             :       ! write data
     843           0 :       newTime = real(evalPer(iDomain)%julStart, dp) - 0.5_dp
     844             : 
     845           0 :       do tt = 1, (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * nMeasPerDay
     846           0 :         call dec2date(newTime, yy = year, mm = month, dd = day, hh = hour)
     847           0 :         write(usubdaily_discharge, formData) tt, hour, day, month, year, (Qobs(tt, gg), Qsim(tt, gg), gg = igauge_start, igauge_end)
     848           0 :         newTime = newTime + 1.0_dp / real(nMeasPerDay, dp)
     849             :       end do
     850             : 
     851             :       ! close file
     852           0 :       close(usubdaily_discharge)
     853             : 
     854             :       ! ======================================================================
     855             :       ! write netcdf file
     856             :       ! ======================================================================
     857           0 :       fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(ncfile_subdaily_discharge))
     858           0 :       nc_out = NcDataset(trim(fName), "w")
     859           0 :       tlength = (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * nMeasPerDay
     860             :       ! write time
     861           0 :       allocate(taxis(tlength))
     862             : 
     863           0 :       use_minutes = .false.
     864           0 :       time_unit_factor = 1
     865           0 :       if ( mod(factor, 2) == 1 ) then
     866           0 :         use_minutes = .true.
     867           0 :         time_unit_factor = 60
     868             :       end if
     869             : 
     870             :       ! tt is dependent on the unit of the time axis and is set to hours in mRM
     871             :       select case( output_time_reference_mrm)
     872             :         case(0)
     873           0 :           forall(tt = 1 : tlength) taxis(tt) = (tt-1) * factor * time_unit_factor
     874             :         case(1)
     875           0 :           forall(tt = 1 : tlength) taxis(tt) = tt * factor - factor * time_unit_factor / 2
     876             :         case(2)
     877           0 :           forall(tt = 1 : tlength) taxis(tt) = tt * factor * time_unit_factor
     878             :       end select
     879             : 
     880           0 :       call dec2date(real(evalPer(iDomain)%julStart, dp) - 0.5_dp, yy = year, mm = month, dd = day, hh = hour)
     881           0 :       dim = nc_out%setDimension("time", tlength)
     882           0 :       var = nc_out%setVariable("time", "i32", [dim])
     883           0 :       call var%setData(taxis)
     884           0 :       if (use_minutes) then
     885             :         call var%setAttribute( &
     886             :           "units", &
     887             :           'minutes since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' '// &
     888             :           trim(num2str(hour, '(i2.2)'))//':00:00' &
     889           0 :         )
     890           0 :         call var%setAttribute("long_name", "time in minutes")
     891             :       else
     892             :         call var%setAttribute( &
     893             :           "units", &
     894             :           'hours since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' '// &
     895             :           trim(num2str(hour, '(i2.2)'))//':00:00' &
     896           0 :         )
     897           0 :         call var%setAttribute("long_name", "time in hours")
     898             :       end if
     899           0 :       call var%setAttribute("bounds", "time_bnds")
     900           0 :       call var%setAttribute("axis", "T")
     901           0 :       dim_bnd = nc_out%setDimension("bnds", 2)
     902           0 :       var = nc_out%setVariable("time_bnds", "i32", [dim_bnd, dim])
     903           0 :       do tt = 1, tlength
     904           0 :         call var%setData((tt - 1) * factor * time_unit_factor, (/1, tt/))
     905           0 :         call var%setData(tt * factor * time_unit_factor, (/2, tt/))
     906             :       end do
     907           0 :       deallocate(taxis)
     908             :       ! write gauges
     909           0 :       do gg = igauge_start, igauge_end
     910           0 :         var = nc_out%setVariable('Qsim_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
     911           0 :         call var%setFillValue(nodata_dp)
     912           0 :         call var%setData(Qsim(1 : tlength, gg))
     913           0 :         call var%setAttribute("units", "m3 s-1")
     914           0 :         call var%setAttribute("long_name", 'simulated discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
     915           0 :         call var%setAttribute("missing_value", nodata_dp)
     916             :         ! write observed discharge at that gauge
     917           0 :         var = nc_out%setVariable('Qobs_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
     918           0 :         call var%setFillValue(nodata_dp)
     919           0 :         call var%setData(Qobs(1 : tlength, gg))
     920           0 :         call var%setAttribute("units", "m3 s-1")
     921           0 :         call var%setAttribute("long_name", 'observed discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
     922           0 :         call var%setAttribute("missing_value", nodata_dp)
     923             :       end do
     924           0 :       call nc_out%close()
     925             : 
     926             :       ! ======================================================================
     927             :       ! screen output
     928             :       ! ======================================================================
     929           0 :       call message()
     930           0 :       write(dummy, '(I3)') domainID
     931           0 :       call message('  OUTPUT: saved subdaily discharge file for domain ', trim(adjustl(dummy)))
     932           0 :       call message('    to ', trim(fname))
     933           0 :       do gg = igauge_start, igauge_end
     934           0 :         if (count(ge(Qobs(:, gg), 0.0_dp)) > 1)  then
     935             :           call message('    KGE of subdaily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
     936           0 :                   trim(adjustl(num2str(kge(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
     937             :           call message('    NSE of subdaily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
     938           0 :                   trim(adjustl(num2str(nse(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
     939             :         end if
     940             :       end do
     941             : 
     942             :       ! update igauge_start
     943           0 :       igauge_start = igauge_end + 1
     944             :       !
     945             :     end do
     946             :     !
     947           0 :   end subroutine write_subdaily_obs_sim_discharge
     948             : 
     949             :   ! ------------------------------------------------------------------
     950             : 
     951             :   !    NAME
     952             :   !        mrm_write_optifile
     953             : 
     954             :   !    PURPOSE
     955             :   !>       \brief Write briefly final optimization results.
     956             : 
     957             :   !>       \details Write overall best objective function and the best optimized parameter set to a file_opti.
     958             : 
     959             :   !    INTENT(IN)
     960             :   !>       \param[in] "real(dp) :: best_OF"                             best objective function value as returnedby the
     961             :   !>       optimization routine
     962             :   !>       \param[in] "real(dp), dimension(:) :: best_paramSet"         best associated global parameter setCalled only
     963             :   !>       when optimize is .TRUE.
     964             :   !>       \param[in] "character(len = *), dimension(:) :: param_names"
     965             : 
     966             :   !    HISTORY
     967             :   !>       \authors David Schaefer
     968             : 
     969             :   !>       \date July 2013
     970             : 
     971             :   ! Modifications:
     972             :   ! Rohini Kumar   Aug 2013 - change in structure of the code including call statements
     973             :   ! Juliane Mai    Oct 2013 - clear parameter names added
     974             :   !                         - double precision written
     975             :   ! Stephan Thober Oct 2015 - ported to mRM
     976             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     977             : 
     978           0 :   subroutine mrm_write_optifile(best_OF, best_paramSet, param_names)
     979             : 
     980           0 :     use mo_common_mhm_mrm_file, only : file_opti, uopti
     981             :     use mo_common_variables, only : dirConfigOut
     982             :     use mo_string_utils, only : num2str
     983             : 
     984             :     implicit none
     985             : 
     986             :     ! best objective function value as returnedby the optimization routine
     987             :     real(dp), intent(in) :: best_OF
     988             : 
     989             :     ! best associated global parameter setCalled only when optimize is .TRUE.
     990             :     real(dp), dimension(:), intent(in) :: best_paramSet
     991             : 
     992             :     character(len = *), dimension(:), intent(in) :: param_names
     993             : 
     994             :     character(256) :: fName, formHeader, formParams
     995             : 
     996             :     integer(i4) :: ii, err, n_params
     997             : 
     998             : 
     999             :     ! number of parameters
    1000           0 :     n_params = size(best_paramSet)
    1001             : 
    1002             :     ! open file
    1003           0 :     fName = trim(adjustl(dirConfigOut)) // trim(adjustl(file_opti))
    1004           0 :     open(uopti, file = fName, status = 'unknown', action = 'write', iostat = err, recl = (n_params + 1) * 40)
    1005           0 :     if(err .ne. 0) then
    1006           0 :       call error_message('  IOError while openening "', trim(fName), '" Error-Code ', num2str(err))
    1007             :     end if
    1008             : 
    1009             :     ! header
    1010           0 :     write(formHeader, *) '(a40,', n_params, 'a40)'
    1011             :     ! len(param_names(1))=256 but only 39 characters taken here
    1012             :     ! write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii))), ii=1, n_params)
    1013           0 :     write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii)(1 : 39))), ii = 1, n_params)
    1014             : 
    1015             :     ! output
    1016           0 :     write(formParams, *) '( es40.14, ', n_params, '(es40.14) )'
    1017           0 :     write(uopti, formParams) best_OF, (best_paramSet(ii), ii = 1, n_params)
    1018             : 
    1019             :     ! close file
    1020           0 :     close(uopti)
    1021             : 
    1022             :     ! screen output
    1023           0 :     call message()
    1024           0 :     call message(' Optimized parameters written to ', trim(fName))
    1025             : 
    1026           0 :   end subroutine mrm_write_optifile
    1027             : 
    1028             :   ! ------------------------------------------------------------------
    1029             : 
    1030             :   !    NAME
    1031             :   !        mrm_write_optinamelist
    1032             : 
    1033             :   !    PURPOSE
    1034             :   !>       \brief Write final, optimized parameter set in a namelist format.
    1035             : 
    1036             :   !>       \details Write final, optimized parameter set in a namelist format.
    1037             : 
    1038             :   !    INTENT(IN)
    1039             :   !>       \param[in] "real(dp), dimension(:, :) :: parameters"                               (min, max, opti)
    1040             :   !>       \param[in] "logical, dimension(size(parameters, 1)) :: maskpara"                   .true. if parameter was
    1041             :   !>       calibrated
    1042             :   !>       \param[in] "character(len = *), dimension(size(parameters, 1)) :: parameters_name" clear names of parameters
    1043             : 
    1044             :   !    HISTORY
    1045             :   !>       \authors Juliane Mai
    1046             : 
    1047             :   !>       \date Dec 2013
    1048             : 
    1049             :   ! Modifications:
    1050             :   ! Stephan Thober Oct 2015 - adapted to mRM
    1051             :   ! Stephan Thober Nov 2016 - adapt header to routing process
    1052             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1053             : 
    1054           0 :   subroutine mrm_write_optinamelist(parameters, maskpara, parameters_name)
    1055             : 
    1056           0 :     use mo_common_mhm_mrm_file, only : file_opti_nml, uopti_nml
    1057             :     use mo_common_variables, only : dirConfigOut, processMatrix
    1058             :     use mo_string_utils, only : num2str
    1059             : 
    1060             :     implicit none
    1061             : 
    1062             :     ! (min, max, opti)
    1063             :     real(dp), dimension(:, :), intent(in) :: parameters
    1064             : 
    1065             :     ! .true. if parameter was calibrated
    1066             :     logical, dimension(size(parameters, 1)), intent(in) :: maskpara
    1067             : 
    1068             :     ! clear names of parameters
    1069             :     character(len = *), dimension(size(parameters, 1)), intent(in) :: parameters_name
    1070             : 
    1071             :     character(256) :: fName
    1072             : 
    1073             :     character(3) :: flag
    1074             : 
    1075             :     integer(i4) :: err
    1076             : 
    1077             :     integer(i4) :: iPar
    1078             : 
    1079             : 
    1080             :     ! open file
    1081           0 :     fName = trim(adjustl(dirConfigOut)) // trim(adjustl(file_opti_nml))
    1082           0 :     open(uopti_nml, file = fName, status = 'unknown', action = 'write', iostat = err)
    1083           0 :     if(err .ne. 0) then
    1084           0 :       call message ('  IOError while openening "', trim(fName), '" Error-Code ', num2str(err))
    1085             :     end if
    1086             : 
    1087           0 :     write(uopti_nml, *) '!global_parameters'
    1088           0 :     write(uopti_nml, *) '!PARAMETER                       lower_bound  upper_bound          value   FLAG  SCALING'
    1089             : 
    1090           0 :     write(uopti_nml, *) '! ', trim(adjustl('routing'))
    1091             : 
    1092           0 :     if (processMatrix(8, 1) .eq. 1_i4) write(uopti_nml, *) '&routing1'
    1093           0 :     if (ProcessMatrix(8, 1) .eq. 2_i4) write(uopti_nml, *) '&routing2'
    1094           0 :     if (ProcessMatrix(8, 1) .eq. 3_i4) write(uopti_nml, *) '&routing3'
    1095             : 
    1096           0 :     do iPar = 1, size(parameters, 1)
    1097           0 :       if (maskpara(iPar)) then
    1098           0 :         flag = ' 1 '
    1099             :       else
    1100           0 :         flag = ' 0 '
    1101             :       end if
    1102           0 :       write(uopti_nml, *) trim(adjustl(parameters_name(iPar))), ' = ', &
    1103           0 :               parameters(iPar, 1), ' , ', &
    1104           0 :               parameters(iPar, 2), ' , ', &
    1105           0 :               parameters(iPar, 3), ' , ', &
    1106           0 :               flag, ', 1 '
    1107             :     end do
    1108             : 
    1109           0 :     write(uopti_nml, *) '/'
    1110           0 :     write(uopti_nml, *) ' '
    1111             : 
    1112             :     ! close file
    1113           0 :     close(uopti_nml)
    1114             : 
    1115             :     ! screen output
    1116           0 :     call message()
    1117           0 :     call message(' Optimized parameters written in namelist format to ', trim(fName))
    1118             : 
    1119           0 :   end subroutine mrm_write_optinamelist
    1120             : 
    1121             : 
    1122             : end module mo_mrm_write

Generated by: LCOV version 1.16