!     FILE writgen.f
!
!     SUBROUTINES IN THIS FILE CALLED FROM
!     subroutine wr_daily main
!     subroutine wr_month(mo1) main
!     subroutine wr_annual main
module output

  use utilities, only : &
    dp, &
    log_debug, &
    log_info, &
    log_error, &
    identifier_max_length, &
    output_max_variables, &
    path_max_length, &
    output_label_max_length

  implicit none

  integer, parameter :: nih = 7
  integer :: nvrch = 18
  integer :: nvsub = 30
  integer :: nns = 30
  !
  integer, save :: ieapu
  ! area covered by snow [m2]
  real(dp), save :: area_tot_snow
  ! average snow depth [mm]
  real(dp), save :: depth_ave_snow
  ! area covered by glaciers [m2]
  real(dp), save :: area_tot_glacier
  !
  real(dp), save :: depth_ave_glacier ! average glacier depth [mm], averaged over the entire catchment
  ! annual SUM of precipitation in subbasin
  real(dp), save, dimension(:), allocatable :: sbpy
  ! annual SUBBASIN outputs (dif. components)
  real(dp), save, dimension(:, :), allocatable :: sysub
  ! total SUBBASIN outputs = SUM(sysub)
  real(dp), save, dimension(:, :), allocatable :: stsub
  integer, save :: giscounter = 1;
  ! annual REACH outputs (dif. components)
  real(dp), save, dimension(:, :), allocatable :: syrch
  ! total REACH outputs = SUM(syrch)
  real(dp), save, dimension(:, :), allocatable :: strch
  integer :: nsb = 30

  ! Directory to write output to
  character(len=path_max_length) :: output_dir = "output"
  ! Output specs namelist path
  character(len=path_max_length) :: output_specs_path = "input/output.nml"
  ! Write out interval in timesteps: daily D, monthly M, annual Y
  character(len=1) :: output_write_interval = "M"
  ! Float csv output format
  character(len=identifier_max_length) :: output_float_format = "f12.4"
  ! Format of space index in csv output
  character(len=identifier_max_length) :: output_space_index_format = "i6"
  ! Default file format if none given in output file nml
  character(len=identifier_max_length) :: output_default_format = "csv"

  ! Log levels and master logfile
  character(len=path_max_length) :: master_logfile = "swim.log"
  character(len=identifier_max_length) :: &
    log_stderr_level = "warning", &
    log_stdout_level = "info", &
    log_file_level = "info"

  ! Output file attributes (mostly as integer indeces of corresponding module variables)
  type output_file
    ! Suffix after <space>_<time>_<name>.<format>
    character(len=identifier_max_length) :: name
    ! Requested variables as columns
    character(len=identifier_max_length) :: variables(output_max_variables)
    ! Indeces of dimensions
    integer :: format, space, time
    ! Number and indeces of variables in output_requested_vars
    integer :: nvars, variable_ix(output_max_variables)
    ! ID (and record for binary output)
    integer :: file_id, current_record = 1
  end type
  ! Output files array
  type(output_file), dimension(:), allocatable :: output_files

  ! Variables requested and registered in the code
  character(len=identifier_max_length), dimension(output_max_variables) :: &
    output_requested_vars = "", &
    output_registered_vars = "", &
    output_hydrotope_vars = "", &
    output_subbasin_vars = ""

  ! Dimensions
  character(len=identifier_max_length), parameter :: &
    output_space_dim(5) = (/"hydrotope      ", "subbasin       ", "catchment      ", &
                            "hydrotope_label", "subbasin_label "/), &
    output_time_dim(3) = (/"daily  ", "monthly", "annual "/), &
    output_format_dim(2) = (/"csv", "bin"/)

  ! Counters
  integer :: &
    output_id_hydrotope_counter = 0, &
    output_id_subbasin_counter = 0, &
    output_store_day_counter = 1, &
    output_ndays_per_month = 0, &
    output_store_month_counter = 1, &
    output_ndays_per_year = 0, &
    output_ndays = 0, &
    output_nmonths = 0, &
    output_nyears = 0, &
    output_nvars = 0, &
    output_ncatchments = 0, &
    output_nsubbasins = 0, &
    output_nhydrotopes = 0, &
    output_nfiles = 0, &
    output_specs_fid = 0


  ! Storage arrays for daily raw variable values (vars, space, time)
  real, allocatable, dimension(:, :, :) :: &
    output_storage_hydrotope, &
    output_storage_subbasin

  ! Space and time aggregated storage arrays (vars, space, time)
  real, allocatable, dimension(:, :, :) :: &
    output_agg_hydrotope_monthly, &
    output_agg_hydrotope_annual, &
    output_agg_subbasin_daily, &
    output_agg_subbasin_monthly, &
    output_agg_subbasin_annual, &
    output_agg_catchment_daily, &
    output_agg_catchment_monthly, &
    output_agg_catchment_annual

  integer, allocatable :: &
    ! Requested output convenience matrix (vars, space, time)
    output_is_requested(:, :, :), &
    ! Variable hydrotope/subbasin index in output_requested_vars, 0 if not
    output_hydrotope_var_ix(:), &
    output_subbasin_var_ix(:), &
    ! Indeces of hydrotope raw storage for hydrotope_daily output
    output_hydrotope_requested_ix(:), &
    ! Indeces variables in storage arrays (nfiles, nvars)
    output_storage_index(:, :), &
    ! Hydrotope/subbasin indeces
    output_hydrotope_id(:), &
    output_subbasin_id(:), &
    output_catchment_id(:), &
    ! Number of hydrotopes in each subbasin to unsparse 2D hydrotope arrays
    output_subbasin_nhydrotopes(:), &
    ! Start index of hydrotope to unsparse 2D hydrotope arrays
    output_subbasin_hydrotope_stix(:), &
    ! Hydrotope-subbasin-catchment maps
    output_hydrotope_subbasin_ix(:), &
    output_subbasin_catchment_ix(:), &
    ! Labelled values indeces in storage
    output_hydrotope_label_ix(:), &
    output_subbasin_label_ix(:)

  ! Unit areas for weighted averages
  real, allocatable, dimension(:) :: &
    output_hydrotope_subbasin_share, &
    output_subbasin_catchment_share, &
    output_catchment_basin_share

  ! Time indeces for csv output
  character(len=10), allocatable :: output_day_ix(:)
  character(len=7), allocatable :: output_month_ix(:)

  ! Hydrotope/subbasin labels
  character(len=output_label_max_length), allocatable, dimension(:) :: &
    output_hydrotope_label, &
    output_subbasin_label

  ! Temporal aggregation per variable: average if true, sum if false
  logical, allocatable :: output_var_state(:)

  namelist / OUTPUT_PARAMETERS / &
    nns, &
    nsb, &
    nvrch, &
    nvsub, &
    output_dir, &
    output_write_interval, &
    output_float_format, &
    output_space_index_format, &
    output_default_format, &
    master_logfile, &
    log_stderr_level, &
    log_stdout_level, &
    log_file_level

contains

  subroutine output_allocate(mb, mch)
    integer, intent(in) :: mb, mch
    allocate(sbpy(mb))
    allocate(strch(18, mch))
    allocate(stsub(30, mb))
    allocate(syrch(18, mch))
    allocate(sysub(30, mb))
    sbpy = 0.
    strch = 0.
    stsub = 0.
    syrch = 0.
    sysub = 0.
  end subroutine output_allocate

  subroutine dealloc_output
   deallocate(sbpy)
   deallocate(strch)
   deallocate(stsub)
   deallocate(syrch)
   deallocate(sysub)
  end subroutine dealloc_output


  integer function output_open_file(file_name) result(file_id)
    ! Abstract output file opener (so far only for formatted output)
    use utilities, only : open_file
    character(len=*), intent(in) :: file_name
    file_id = open_file(trim(output_dir)//"/"//trim(file_name), "w")
  end function output_open_file

  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  !     PARAMETERS & VARIABLES
  !
  !     susb(30, j) = monthly SUBBASIN outputs (dif. components) sysub
  !     sysub(30, j) = annual SUBBASIN outputs (dif. components)
  !     stsub(30, j) = total SUBBASIN outputs = SUM(sysub)

  !     sub(30) = daily BASIN outputs (dif. components, weighted sums)
  !     smm(30) = monthly BASIN outputs (dif. components)
  !     smy(30) = annual BASIN outputs (dif. components)
  !     sm(30) = average annual BASIN outputs (dif. components)

  !     srch(18, ir) = monthly REACH outputs (dif. components)
  !     syrch(18, ir)= annual REACH outputs (dif. components)
  !     strch(18, ir)= total REACH outputs = SUM(syrch)

  !     sbp(j) = monthly SUM of precipitation in subbasin
  !     smq(j) = monthly SUM of surface runoff for subbasin
  !     smsq(j) = monthly SUM of subsurface runoff for subbasin
  !     sym(j) = monthly SUM of sediment yield for subbasins
  !     syq(j) = SUM(smq) annual SUM surface runoff for subbasins
  !     sysq(j) = SUM(smsq) annual SUM sub-surface runoff for subbasins
  !     syy(j) = SUM(sym) annual SUM sed. yield for subbasins
  !     sq(j) = SUM(syq) total SUM (whole period) surf. runoff in subb.
  !     ssq(j) = SUM(sysq) total SUM (whole period) sub-surf. runoff in subb.
  !     sy(j) = SUM(syy) total SUM (whole period) sed. yield in subb.

  ! DESCRIPTION of BASIN OUTPUTS sub() (Analogue - for susb, smm, smy):
  ! sub(1) = SUM(precip) precipitation, mm
  ! sub(2) = SUM(precip), if tx()<0. snowfall, mm
  ! sub(3) = SUM(sml) snowmelt(mm)
  ! sub(4) = SUM(tmx*flu) average "weighted" max temp.
  ! sub(5) = SUM(tmn*flu) average "weighted" min temp.
  ! sub(6) = SUM(ra*flu) radiation
  ! sub(7) = SUM(sumcn*flu) curve number
  ! sub(8) = SUM(xqi*flu) surface runoff, mm,
  ! sub(9) = SUM(xssf*flu) subsurface runoff, mm,
  ! sub(10) = SUM(qtl*flu) channel transm. losses, mm
  ! sub(11) = SUM(xsep*flu) percolation, mm,
  ! sub(12) = SUM(xeo*flu) pot evapotr., mm
  ! sub(13) = SUM(xet*flu) evapotranspiration, mm
  ! sub(14) = SUM(snoev*flu) snow evaporation, mm
  ! sub(15) = SUM(gwq*flu) gr. water q, mm
  ! sub(16) = SUM(revap*flu) revap from g-w to soil prof., mm
  ! sub(17) = SUM(gwseep*flu) g-w percolation
  ! sub(18) = SUM(gwchrg*flu) g-w recharge
  ! sub(19) = SUM(xswimd*flu) soil water
  ! sub(20) = SUM(wysb*flu) wysb=qi+ssf+gwq(j)-qtl, water yield, mm
  ! sub(21) = SUM(yd/(100*da*flu) sed yield
  ! sub(22) = SUM(yon*flu) org. N loss with sed
  ! sub(23) = SUM(yph*flu) org P loss with sed
  ! sub(24) = SUM(ysp*flu) sol P
  ! sub(25) = SUM(xyno3*flu) no3 in sur.
  ! sub(26) = SUM(xssfn*flu) no3 in sub sur
  ! sub(27) = SUM(xpercn*flu) no3 leached
  ! sub(28) = SUM(xuno3*flu) N uptake by plants
  !         Analogue - for susb, smm, smy:
  !         SUB(1)= DAILY PRECIPITATION
  !         SMM(1)= TOTAL MONTHLY PRECIPITATION
  !         SMY(1)= TOTAL YEARLY PRECIPITATION
  ! DESCRIPTION OF monthly REACH OUTPUTS
  ! srch(1, ir) = SUM(varoute(2, )/86400.) surface flow, inflow
  ! srch(2, ir) = SUM(xxqd/86400.) surface flow, outflow
  ! srch(3, ir) = SUM(varoute(3, )) sediment in
  ! srch(4, ir) = SUM(yd) sediment out
  ! srch(5, ir) = SUM(sedcon) sediment conc.
  ! srch(6, ir) = SUM(varoute(4, )) organic N in
  ! srch(7, ir) = SUM(yon*dart(ihout)*100) organic N out
  ! srch(8, ir) = SUM(varoute(5, 2)) organic P in
  ! srch(9, ir) = SUM(yph*dart(ihout)*100) organic P out
  ! srch(10, ir) = SUM(evp/86400.) evaporation (not active!)
  ! srch(11, ir) = SUM(tlc/86400. transmission losses (not active!)
  ! srch(12, ir) = SUM(rl/86400. seepage (not active!)
  ! srch(13, ir) = SUM(diver/86400.) diversion (not active!)
  ! srch(14, ir) = SUM(rflow/86400.) return flow (not active!)
  ! srch(15, ir) = SUM(varoute(6, 2)) nitrate N in
  ! srch(16, ir) = SUM(xxnit) nitrate N out
  ! srch(17, ir) = SUM(varoute(7, 2)*100) soluble P in
  ! srch(18, ir) = SUM(xysp*dart(ihout)*100) soluble P out
  !          Analogue - for syrch = SUM(srch) - annual sums
  !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine output_nashsutcliffe_efficiency(qo, qs, icd, istyr, iy)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES CRITERIA OF FIT:
    !              difference in calulated water balance, relative difference**2,
    !              Nash&Sutcliffe Efficiency & LOG-Nash&Sutcliffe Efficiency
    !**** CALLED IN: MAIN
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> PARAMETERS & VARIABLES in TITLE
    !      qo(1:inn) = observed water discharge, m3/sec.
    !      qs(1:inn) = simulated water discharge, m3/sec.
    !      inn = number of days
    !      icd = code: 1 - one year, 2 - whole period

    !      >>>>> COMMON PARAMETERS
    !      istyr = starting year
    !      iy = current year
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      akk(20) = accumulated values, internal
    !      crdif = difference in calulated water balance
    !      crdifp = difference in calulated water balance, %
    !      dif = difference between qs and qo
    !      eff = Nash & Sutcliffe Efficiency for variables qo and qs
    !      efflog = Nash & Sutcliffe Efficiency for log(qo) and log(qs)
    !      f00 = local par
    !      f11 = local par
    !      ik = local par
    !      im = local par
    !      qmid = local par
    !      qmid2 = local par
    !      reldif2 = relative difference **2
    !      xdif = xqo - xqs
    !      xqo = log(qo)
    !      xqs = log(qs)
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** Include common parameters & descriptions

    integer, intent(in) :: istyr, icd, iy
    real(dp), dimension(:), intent(in) :: qo, qs

    integer ik, im
    real(dp) efflog, f00, f11, qmid, qmid2, reldif2, xdif, xqo, xqs, akk(20), crdif, crdifp, dif, eff

    eff = 0.
    efflog = 0.

    do im = 1, 20
      akk(im) = 0.
    end do

    do ik = 1, size(qo)
      dif = qs(ik) - qo (ik)
      qmid = 0.5 * (abs(qo(ik)) + qs(ik))
      qmid2 = qmid * qmid
      akk(1) = akk(1) + 1.
      akk(2) = akk(2) + qo(ik)
      akk(3) = akk(3) + qo(ik) * qo(ik)
      akk(4) = akk(4) + dif
      akk(5) = akk(5) + dif * dif
      if (qmid2 .gt. 0.001) then
        akk(6) = akk(6) + dif * dif / qmid2
      else
      endif

      if (qo(ik) .gt. 0. .and. qs(ik) .gt. 0.) then
        xqo = log(qo(ik))
        xqs = log(qs(ik))
        xdif = xqo - xqs
        akk(7) = akk(7) + 1.
        akk(8) = akk(8) + xqo
        akk(9) = akk(9) + xqo * xqo
        akk(10) = akk(10) + xdif * xdif
      endif
    end do

    crdif = akk(4)

    if (akk(2) .ne. 0) then
      crdifp = akk(4) * 100 / akk(2)
    else
      crdifp = 0.
    endif

    reldif2 = akk(6)

    if (akk(1) .gt. 0.) then
      f00 = akk(3) - akk(2) * akk(2) / akk(1)
      if (f00 .gt. 0.0001) eff = (f00 - akk(5)) / f00
    endif

    if (akk(7) .gt. 0.) then
      f11 = akk(9) - akk(8) * akk(8) / akk(7)
      if (f11 .gt. 0.0001) efflog = (f11 - akk(10)) / f11
    endif

    if (icd .eq. 1) then
      call log_debug("output_nashsutcliffe_efficiency", &
        'CRITERIA of fit in (year): Abs. difference, Difference %, Rel.dif**2, L-Efficiency, Efficiency:', &
        i1=istyr+iy-1, reals=(/crdif, crdifp, reldif2, efflog, eff/))
    else
      call log_debug("output_nashsutcliffe_efficiency", &
        'CRITERIA of fit TOTAL: Abs. difference, Difference %, Rel.dif**2, L-Efficiency, Efficiency:', &
        reals=(/crdif, crdifp, reldif2, efflog, eff/))
    endif

    return
  end subroutine output_nashsutcliffe_efficiency

  subroutine output_initialise(hydrotope_id, hydrotope_area, hydrotope_subbasin, &
    subbasin_id, subbasin_nhydrotopes, subbasin_catchment, catchment_id, &
    hydrotope_input_fid, subbasin_input_fid)
    ! Process user input before variables are registered
    use utilities, only: string_index
    integer, dimension(:), intent(in) :: &
      ! Hydrotope IDs and subbasin mapping (must be same length)
      hydrotope_id, hydrotope_subbasin, &
      ! Subbasin IDs and catchment mapping (must be same length)
      subbasin_id, subbasin_nhydrotopes, subbasin_catchment, &
      !Catchment IDs
      catchment_id
    ! Hydrotope area for weighted averages
    real, dimension(:), intent(in) :: hydrotope_area
    ! File ids to read output labels
    integer, intent(in) :: hydrotope_input_fid
    integer, intent(in) :: subbasin_input_fid
    integer s

    ! Hydrotope arrays
    output_nhydrotopes = size(hydrotope_id)
    allocate(output_hydrotope_id(output_nhydrotopes))
    allocate(output_hydrotope_subbasin_ix(output_nhydrotopes))
    allocate(output_hydrotope_subbasin_share(output_nhydrotopes))
    output_hydrotope_id = hydrotope_id

    ! Subbasin arrays
    output_nsubbasins = size(subbasin_id)
    allocate(output_subbasin_id(output_nsubbasins))
    allocate(output_subbasin_catchment_ix(output_nsubbasins))
    ! initialised with 2 because if catchment ids not read with subcatch = False
    ! 2 is the index of the default catchment ID, 1 being the total catchment
    output_subbasin_catchment_ix = 2
    allocate(output_subbasin_catchment_share(output_nsubbasins))
    output_subbasin_id = subbasin_id
    allocate(output_subbasin_nhydrotopes(output_nsubbasins))
    output_subbasin_nhydrotopes = subbasin_nhydrotopes
    allocate(output_subbasin_hydrotope_stix(output_nsubbasins))
    output_subbasin_hydrotope_stix = 0
    ! Make cumsum out of sparse hydrotope numbers for quick start index
    do s = 1, output_nsubbasins-1
      output_subbasin_hydrotope_stix(s+1) = &
        output_subbasin_hydrotope_stix(s) + subbasin_nhydrotopes(s)
    end do

    ! Catchment arrays
    output_ncatchments = size(catchment_id) + 1  ! one extra for catchment total
    allocate(output_catchment_id(output_ncatchments))
    allocate(output_catchment_basin_share(output_ncatchments))
    output_catchment_id(1) = 0
    output_catchment_id(2:output_ncatchments) = catchment_id

    allocate(output_storage_index(output_nfiles, output_nvars))
    output_storage_index = 0

    call output_initialise_unit_shares(hydrotope_subbasin, subbasin_catchment, hydrotope_area)
    call output_initialise_time
    ! Fill convenience matrix
    call output_initialise_is_requested
    call output_initialise_storage(hydrotope_input_fid, subbasin_input_fid)

    call output_allocate(output_nsubbasins, output_nsubbasins)

  end subroutine output_initialise

  subroutine output_initialise_unit_shares( &
    hydrotope_subbasin, subbasin_catchment, hydrotope_area)
    ! Get array indeces of subbasin and catchment maps and create areal weights
    integer, dimension(:), intent(in) :: hydrotope_subbasin, subbasin_catchment
    real, dimension(:), intent(in) :: hydrotope_area
    integer h, s, c

    output_catchment_basin_share = 0
    output_subbasin_catchment_share = 0
    ! indeces and area sums
    do s = 1, output_nsubbasins
      do h = 1, output_nhydrotopes
        if (hydrotope_subbasin(h) == output_subbasin_id(s)) then
          output_hydrotope_subbasin_ix(h) = s
          output_subbasin_catchment_share(s) = &
            output_subbasin_catchment_share(s) + hydrotope_area(h)
        end if
      end do
    end do
    do c = 2, output_ncatchments
      do s = 1, output_nsubbasins
        if (subbasin_catchment(s) == output_catchment_id(c)) then
          output_subbasin_catchment_ix(s) = c
          output_catchment_basin_share(c) = output_catchment_basin_share(c) &
            + output_subbasin_catchment_share(s)
        end if
      end do
    end do
  ! Make shares
  output_hydrotope_subbasin_share = hydrotope_area / &
    output_subbasin_catchment_share(output_hydrotope_subbasin_ix)
  output_subbasin_catchment_share = output_subbasin_catchment_share &
    / output_catchment_basin_share(output_subbasin_catchment_ix)
  output_catchment_basin_share(1) = sum(hydrotope_area)  ! entire catchment
  output_catchment_basin_share = output_catchment_basin_share &
    / output_catchment_basin_share(1)
  end subroutine output_initialise_unit_shares

  subroutine output_initialise_user_input
    ! Read and interpret output namelist parameters
    use utilities, only: &
      extend_unique_string, &
      string_index, &
      open_file, &
      log_create, &
      master_log
    use input, only: get_config_fid
    integer ios, i, v
    character(len=identifier_max_length) :: &
      name, format, variables(output_max_variables), space, time
    namelist / file / &
      name, format, variables, space, time

    ! General options from central config nml
    read(get_config_fid(), nml=OUTPUT_PARAMETERS)

    ! Create master log file
    master_log = log_create(trim(output_dir)//"/"//trim(master_logfile), &
      log_stderr_level, log_stdout_level, log_file_level)

    ! count requested output files
    output_specs_fid = open_file(output_specs_path)
    output_nfiles = -1
    ios = 0
    do while (ios /= -1)
      output_nfiles = output_nfiles + 1
      read(output_specs_fid, nml=file, iostat=ios)
    end do
    rewind(output_specs_fid)
    if (output_nfiles == 0) &
      call log_error("output_initialise_user_input", &
        "No output file definition found in "//output_specs_path)
    ! Allocate arrays of nfiles length
    allocate(output_files(output_nfiles))
    do i = 1, output_nfiles
      name = ''
      format = output_default_format
      variables = ''
      space = ''
      time = ''
      read(output_specs_fid, nml=file)
      ! Check input and store
      if (trim(name) == '') &
        call log_error("output_initialise_user_input", "No name given for output file.")
      output_files(i)%name = name
      output_files(i)%nvars = string_index('', variables) - 1
      if (output_files(i)%nvars == 0) &
        call log_error("output_initialise_user_input", &
          "No variables define for output file "//trim(name))
      output_files(i)%variables = variables
      output_files(i)%space = string_index(space, output_space_dim)
      if (output_files(i)%space == 0) call log_error("output_initialise_user_input", &
        trim(space)//" is not a valid output space domain in file "//trim(name))
      output_files(i)%time = string_index(time, output_time_dim)
      if (output_files(i)%time == 0) call log_error("output_initialise_user_input", &
        trim(time)//" is not a valid output time domain in file "//trim(name))
      output_files(i)%format = string_index(format, output_format_dim)
      if (output_files(i)%format == 0) call log_error("output_initialise_user_input", &
        trim(format)//" is not a valid output format in file "//trim(name))
      ! Store unique variables
      call extend_unique_string(variables, output_requested_vars)
    end do
    ! Get variable indeces for each file%variables in output_requested_vars
    do i = 1, output_nfiles
      do v = 1, output_files(i)%nvars
        output_files(i)%variable_ix(v) = string_index(output_files(i)%variables(v), &
          output_requested_vars)
      end do
    end do
    ! Count all requested variables
    output_nvars = string_index("", output_requested_vars) - 1
    allocate(output_var_state(output_nvars))
    output_var_state = .True.
  end subroutine output_initialise_user_input

  subroutine output_initialise_is_requested
    type(output_file) :: file
    integer i, v, sp

    ! Allocate without labelled space units as they are just copied
    allocate(output_is_requested(output_nvars, 3, size(output_time_dim)))
    output_is_requested = 0
    ! Marke 1 where needed in (var, space, time) matrix
    do i = 1, output_nfiles
      file = output_files(i)
      do v = 1, file%nvars
        sp = file%space
        ! Make labelled count the same as unlabelled
        if (sp > 3) sp = sp - 3
        output_is_requested(file%variable_ix(v), sp, file%time) = 1
      end do
    end do
  end subroutine output_initialise_is_requested

  integer function output_register_hydrotope_var(name, state)
    ! Register hydrotope variable name and optionally "state" (default .True.).
    ! Returns output ID if variable requested for output, otherwise 0.
    ! If state = .False. temporal aggregation is sum not average.
    character(*), intent(in) :: name
    logical, intent(in), optional :: state
    logical s

    s = .True.
    if (present(state)) s = state
    call output_register_var(name, s, output_id_hydrotope_counter, &
      output_hydrotope_vars, output_register_hydrotope_var)
  end function output_register_hydrotope_var

  integer function output_register_subbasin_var(name, state)
    ! Register subbasin variable name and optionally "state" (default .True.).
    ! Returns output ID if variable requested for output, otherwise 0.
    ! If state = .False. temporal aggregation is sum not average.
    character(*), intent(in) :: name
    logical, intent(in), optional :: state
    logical s

    s = .True.
    if (present(state)) s = state
    call output_register_var(name, s, output_id_subbasin_counter, &
      output_subbasin_vars, output_register_subbasin_var)
  end function output_register_subbasin_var

  subroutine output_register_var(name, state, counter, vars, id)
    ! Abstract subroutine of register_hydrotope/subbasin_var functions.
    use utilities, only: string_index
    character(*), intent(in) :: name
    logical, intent(in) :: state
    integer, intent(inout) :: counter
    character(*), intent(out) :: vars(:)
    integer, intent(out) :: id
    integer vi

    ! make sure name hasnt been registered already
    if (string_index(name, vars) > 0) &
      call log_error("output_register_var", trim(name)//" has already been registered.")
    ! is variable requested?
    vi = string_index(name, output_requested_vars)
    output_registered_vars(string_index("", output_registered_vars)) = name
    ! increase counter if requested
    if (vi > 0) then
      counter = counter + 1
      vars(counter) = name
      id = counter
      output_var_state(vi) = state
    else
      id = 0
    end if
  end subroutine output_register_var

  subroutine output_initialise_storage(hydrotope_input_fid, subbasin_input_fid)
    ! Initialisation storage after all variables have been registered registration.
    integer, intent(in) :: hydrotope_input_fid, subbasin_input_fid
    integer :: nd, nm, ny
    integer :: nspacetime(3, 3)

    call output_initialise_labelled(hydrotope_input_fid, subbasin_input_fid)
    call output_check_requested_vars
    call output_open_files
    call output_initialise_storage_index

    nd = output_ndays
    nm = output_nmonths
    ny = output_nyears
    ! Raw daily stores
    allocate(output_storage_hydrotope(output_id_hydrotope_counter, output_nhydrotopes, nd))
    allocate(output_storage_subbasin(output_id_subbasin_counter, output_nsubbasins, nd))
    ! Total variables requested for (space, time), maxval because [csv, bin] = 1
    nspacetime = sum(output_is_requested, dim=1)
    allocate(output_agg_hydrotope_monthly(nspacetime(1, 2), output_nhydrotopes, nm))
    allocate(output_agg_hydrotope_annual(nspacetime(1, 3), output_nhydrotopes, ny))
    allocate(output_agg_subbasin_daily(nspacetime(2, 1), output_nsubbasins, nd))
    allocate(output_agg_subbasin_monthly(nspacetime(2, 2), output_nsubbasins, nm))
    allocate(output_agg_subbasin_annual(nspacetime(2, 3), output_nsubbasins, ny))
    allocate(output_agg_catchment_daily(nspacetime(3, 1), output_ncatchments, nd))
    allocate(output_agg_catchment_monthly(nspacetime(3, 2), output_ncatchments, nm))
    allocate(output_agg_catchment_annual(nspacetime(3, 3), output_ncatchments, ny))
    ! set all to 0
    output_storage_hydrotope = 0
    output_storage_subbasin = 0
    output_agg_hydrotope_monthly = 0
    output_agg_hydrotope_annual = 0
    output_agg_subbasin_daily = 0
    output_agg_subbasin_monthly = 0
    output_agg_subbasin_annual = 0
    output_agg_catchment_daily = 0
    output_agg_catchment_monthly = 0
    output_agg_catchment_annual = 0
  end subroutine output_initialise_storage

  subroutine output_initialise_time
    ! decide time dimension/index
    output_ndays = 1
    output_nmonths = 1
    output_nyears = 1
    if (output_write_interval == "M") then
      output_ndays = 31
    else if (output_write_interval == "Y") then
      output_ndays = 366
      output_nmonths = 12
    else if (output_write_interval /= "D") then
      call log_error("output_initialise_time", &
        "output_write_interval must be one of D, M, Y not "//trim(output_write_interval))
    end if
    allocate(output_day_ix(output_ndays))
    allocate(output_month_ix(output_nmonths))
  end subroutine output_initialise_time

  subroutine output_initialise_storage_index
    use utilities, only: string_index, indeces_in_array
    integer fi, sp, i, cumsum(output_nvars)
    type(output_file) :: fl

    do fi = 1, output_nfiles
      fl = output_files(fi)
      sp = fl%space
      if (sp > 3) sp = sp - 3
      cumsum = 0
      do i = 1, output_nvars
        cumsum(i) = sum(output_is_requested(:i, sp, fl%time))
      end do
      output_storage_index(fi, :fl%nvars) = cumsum(fl%variable_ix(:fl%nvars))
    end do
    ! hydrotope/subbasin storage indeces
    allocate(output_hydrotope_var_ix(output_nvars))
    output_hydrotope_var_ix = 0
    allocate(output_subbasin_var_ix(output_nvars))
    output_subbasin_var_ix = 0
    allocate(output_hydrotope_requested_ix(sum(output_is_requested(:, 1, 1))))
    output_hydrotope_requested_ix = 0

    ! Hydrotope/subbasin indeces
    output_hydrotope_requested_ix = indeces_in_array( &
      pack(output_requested_vars(:output_nvars), output_is_requested(:, 1, 1) == 1), &
      output_hydrotope_vars)
    output_hydrotope_var_ix = indeces_in_array( &
      output_requested_vars(:output_nvars), output_hydrotope_vars)
    output_subbasin_var_ix = indeces_in_array(&
      output_requested_vars(:output_nvars), output_subbasin_vars)

  end subroutine output_initialise_storage_index

  subroutine output_initialise_labelled(hydrotope_input_fid, subbasin_input_fid)
    use utilities, only: string_index
    !use input, only: read_string_column
    use input, only: read_string_column, has_column
    integer, intent(in) :: hydrotope_input_fid, subbasin_input_fid
    integer i, ucount
    logical ir

    ir = any(output_files(:)%space == 4)
    if (ir .and. has_column("output_label", hydrotope_input_fid)) then
      allocate(output_hydrotope_label(output_nhydrotopes))
      output_hydrotope_label = ''
      call read_string_column(hydrotope_input_fid, 'output_label', output_hydrotope_label)
      ucount = 0
      do i = 1, output_nhydrotopes
        if (trim(output_hydrotope_label(i)) /= '') ucount = ucount + 1
      end do
      allocate(output_hydrotope_label_ix(ucount))
      ucount = 0
      do i = 1, output_nhydrotopes
        if (trim(output_hydrotope_label(i)) /= '') then
          ucount = ucount + 1
          output_hydrotope_label_ix(ucount) = i
        end if
      end do
    else if (ir) then
      call log_error("output_initialise_labelled", &
        "hydrotope_label output requested but output_label column not found.")
    end if

    ir = any(output_files(:)%space == 5)
    if (ir .and. has_column("output_label", subbasin_input_fid)) then
      allocate(output_subbasin_label(output_nsubbasins))
      output_subbasin_label = ''
      call read_string_column(subbasin_input_fid, 'output_label', output_subbasin_label)
      ucount = 0
      do i = 1, output_nsubbasins
        if (trim(output_subbasin_label(i)) /= '') ucount = ucount + 1
      end do
      allocate(output_subbasin_label_ix(ucount))
      ucount = 0
      do i = 1, output_nsubbasins
        if (trim(output_subbasin_label(i)) /= '') then
          ucount = ucount + 1
          output_subbasin_label_ix(ucount) = i
        end if
      end do
    else if (ir) then
      call log_error("output_initialise_labelled", &
        "subbasin_label output requested but output_label column not found.")
    end if
  end subroutine output_initialise_labelled

  subroutine output_open_files
    ! Open all needed output files and write csv headers
    use utilities, only : open_file
    integer fi, fid, rl
    type(output_file) :: fl
    character(len=path_max_length) :: fname
    integer, parameter :: nbytes_per_real = 4

    do fi = 1, output_nfiles
      fl = output_files(fi)
      fname = trim(output_space_dim(fl%space))//"_"//&
        trim(output_time_dim(fl%time))//"_"//trim(fl%name)// &
        "."//trim(output_format_dim(fl%format))
      if (fl%format == 1) then
        fid = output_open_file(fname)
      else if (fl%format == 2) then
        rl = output_nhydrotopes
        if (fl%space == 2) rl = output_nsubbasins
        if (fl%space == 3) rl = output_ncatchments
        if (fl%space == 4) rl = size(output_hydrotope_label_ix)
        if (fl%space == 5) rl = size(output_subbasin_label_ix)
        rl = rl * nbytes_per_real * fl%nvars
        ! TODO implement binary open file in open_file
        open(fid, file=fname, form="unformatted", access="direct", recl=rl)
      end if
      output_files(fi)%file_id = fid
      ! Write header if needed
      if (fl%format == 1) &
        call output_write_csv_header(fid, output_space_dim(fl%space), &
          output_requested_vars(fl%variable_ix(:fl%nvars)))
    end do
  end subroutine output_open_files

  subroutine output_write_csv_header(fid, space, variables)
    ! Write csv header <time> <space> <variables>
    integer, intent(in) :: fid
    character(len=identifier_max_length), intent(in) :: space, variables(:)
    character(len=identifier_max_length) :: fmt
    integer i, ffl

    ! Format lengths as int
    read(output_float_format(2:index(output_float_format, ".") - 1), *) ffl
    ! time
    write(fid, "(a)", advance="no") "time, "
    ! space
    write(fid, "(a)", advance="no") trim(space)//","
    ! variables
    do i = 1, size(variables) - 1
      ! Make sure variable name fits into the size of real format
      write(fmt, "('(a',i4,')')") max(ffl, len(trim(variables(i)))+1) + 1
      write(fid, trim(fmt), advance="no") trim(variables(i))//","
    end do
    ! Last column
    write(fmt, "('(a',i4,')')") max(ffl, len(trim(variables(size(variables))))+1)
    write(fid, fmt) trim(variables(size(variables)))
  end subroutine output_write_csv_header

  subroutine output_check_requested_vars
    ! Make sure all requested variables have been registered.
    use utilities, only: string_index
    character(len=identifier_max_length) vs
    integer i, ih, is

    do i = 1, output_nvars
      vs = output_requested_vars(i)
      ih = string_index(vs, output_hydrotope_vars)
      is = string_index(vs, output_subbasin_vars)
      ! Is requested variable registered?
      if (ih + is == 0) call log_error("output_check_requested_vars", &
        "Requested output variable invalid: "//vs)
      ! Is subbasin variable requested at hydrotope level?
      if (ih == 0 .and. maxval(output_is_requested(i, 1, :)) > 0) then
        call log_error("output_check_requested_vars", &
          "Requested variable at hydrotope level is not a hydrotope variable: "//vs)
      end if
    end do
  end subroutine output_check_requested_vars

  subroutine output_store_hydrotope_values(id, array)
    ! Store hydrotope values at the end of timestep.
    integer, intent(in) :: id
    real, dimension(:, :), intent(in) :: array
    integer s, ns, nh

    if (id > 0) then
      ! Unsparse hydrotope data
      do s = 1, output_nsubbasins
        ns = output_subbasin_hydrotope_stix(s)
        nh = output_subbasin_nhydrotopes(s)
        output_storage_hydrotope(id, ns:ns+nh-1, output_store_day_counter) = array(s, :nh)
      end do
    end if
  end subroutine output_store_hydrotope_values

  subroutine output_store_subbasin_values(id, array)
    ! Store subbasin values at the end of timestep.
    integer, intent(in) :: id
    real, dimension(:), intent(in) :: array

    if (id > 0) output_storage_subbasin(id, :, output_store_day_counter) = array
  end subroutine output_store_subbasin_values

  subroutine output_store_hydrotope_value(id, subbasinix, hydrotopeix, value)
    ! Store hydrotope values at the end of timestep.
    integer, intent(in) :: id
    ! Subbasin and hydrotope (within subbasin) indeces
    integer, intent(in) :: subbasinix, hydrotopeix
    real, intent(in) :: value
    integer ns

    if (id > 0) then
      ns = output_subbasin_hydrotope_stix(subbasinix)
      output_storage_hydrotope(id, ns+hydrotopeix, output_store_day_counter) = value
    end if
  end subroutine output_store_hydrotope_value

  subroutine output_store_subbasin_value(id, subbasinix, value)
    ! Store a single subbasin value at the end of timestep.
    integer, intent(in) :: id, subbasinix
    real, intent(in) :: value

    if (id > 0) output_storage_subbasin(id, subbasinix, output_store_day_counter) = value
  end subroutine output_store_subbasin_value

  function output_hydrotope_to_subbasin(hydrotope_values) result(subbasin_values)
    ! Average hydrotope values to subbasins
    real, intent(in) :: hydrotope_values(output_nhydrotopes)
    real :: subbasin_values(output_nsubbasins)
    integer h, hsix
    real share

    subbasin_values = 0
    do h = 1, output_nhydrotopes
      hsix = output_hydrotope_subbasin_ix(h)
      share = hydrotope_values(h) * output_hydrotope_subbasin_share(h)
      subbasin_values(hsix) = subbasin_values(hsix) + share
    end do
  end function output_hydrotope_to_subbasin

  function output_hydrotope_to_catchment(hydrotope_values) result(catchment_values)
    ! Average hydrotope values to catchment
    real, intent(in) :: hydrotope_values(output_nhydrotopes)
    real :: catchment_values(output_ncatchments)
    integer h, hsix, scix
    real share

    catchment_values = 0
    do h = 1, output_nhydrotopes
      hsix = output_hydrotope_subbasin_ix(h)
      scix = output_subbasin_catchment_ix(hsix)
      ! Share of hydrotope in catchment
      share = hydrotope_values(h) * output_hydrotope_subbasin_share(h) &
        * output_subbasin_catchment_share(hsix)
      catchment_values(scix) = catchment_values(scix) + share
      ! Entire basin
      catchment_values(1) = &
        catchment_values(1) + share * output_catchment_basin_share(scix)
    end do
  end function output_hydrotope_to_catchment

  function output_subbasin_to_catchment(subbasin_values) result(catchment_values)
    ! Average subbasin values to catchment
    real, intent(in) :: subbasin_values(output_nsubbasins)
    real :: catchment_values(output_ncatchments)
    integer s, scix
    real share


    catchment_values = 0
    do s = 1, output_nsubbasins
      scix = output_subbasin_catchment_ix(s)
      share = subbasin_values(s) * output_subbasin_catchment_share(s)
      catchment_values(scix) = catchment_values(scix) + share
      catchment_values(1) = &
        catchment_values(1) + share * output_catchment_basin_share(scix)
    end do
  end function output_subbasin_to_catchment

  subroutine output_day(year, month, day)
    ! Aggregate space/time and write to csv/bin, if interval is daily.
    integer, intent(in) :: year, month, day
    integer i, dc, mc
    integer :: sti(3, 3), stb(3, 3)
    real :: hydval(output_nhydrotopes), subval(output_nsubbasins), &
      catval(output_ncatchments)

    dc = output_store_day_counter
    mc = output_store_month_counter
    ! Format day label for csv output
    write(output_day_ix(dc), "(i0.4, '-', i0.2, '-', i0.2)") year, month, day
    ! aggregate day to month
    do i = 1, output_nvars
      ! 0/1 (space, time) for current variable
      stb = output_is_requested(i, :, :)
      ! indeces in aggregation arrays of current variable
      sti = sum(output_is_requested(:i, :, :), dim=1)
      ! current hydrotope/subbasin variables
      if (output_hydrotope_var_ix(i) > 0) then
        hydval = output_storage_hydrotope(output_hydrotope_var_ix(i), :, dc)
        subval = output_hydrotope_to_subbasin(hydval)
        catval = output_hydrotope_to_catchment(hydval)
      end if
      if (output_subbasin_var_ix(i) > 0) then
        hydval = 0
        subval = output_storage_subbasin(output_subbasin_var_ix(i), :, dc)
        catval = output_subbasin_to_catchment(subval)
      end if
      ! Add daily values to aggregation arrays as needed
      if (stb(1, 2) == 1) output_agg_hydrotope_monthly(sti(1, 2), :, mc) = &
        output_agg_hydrotope_monthly(sti(1, 2), :, mc) + hydval
      if (stb(1, 3) == 1) output_agg_hydrotope_annual(sti(1, 3), :, 1) = &
        output_agg_hydrotope_annual(sti(1, 3), :, 1) + hydval
      if (stb(2, 1) == 1) output_agg_subbasin_daily(sti(2, 1), :, dc) = &
        output_agg_subbasin_daily(sti(2, 1), :, dc) + subval
      if (stb(2, 2) == 1) output_agg_subbasin_monthly(sti(2, 2), :, mc) = &
        output_agg_subbasin_monthly(sti(2, 2), :, mc) + subval
      if (stb(2, 3) == 1) output_agg_subbasin_annual(sti(2, 3), :, 1) = &
        output_agg_subbasin_annual(sti(2, 3), :, 1) + subval
      if (stb(3, 1) == 1) output_agg_catchment_daily(sti(3, 1), :, dc) = &
        output_agg_catchment_daily(sti(3, 1), :, dc) + catval
      if (stb(3, 2) == 1) output_agg_catchment_monthly(sti(3, 2), :, mc) = &
        output_agg_catchment_monthly(sti(3, 2), :, mc) + catval
      if (stb(3, 3) == 1) output_agg_catchment_annual(sti(3, 3), :, 1) = &
        output_agg_catchment_annual(sti(3, 3), :, 1) + catval
    end do
    ! Increase/reset day counter
    if (output_write_interval == "D") then
      call output_write_daily
      output_store_day_counter = 0
    end if
    output_store_day_counter = output_store_day_counter + 1
    output_ndays_per_month = output_ndays_per_month + 1
  end subroutine output_day

  subroutine output_month(year, month)
    ! Take monthly averages and write to csv/bin, if interval is monthly.
    integer, intent(in) :: year, month
    integer :: sti(3, 3), stb(3, 3)
    integer i, nd

    ! Format month label for csv output
    write(output_month_ix(output_store_month_counter), "(i0.4, '-', i0.2)") year, month
    ! Take average of monthly values as needed
    nd = output_ndays_per_month
    do i = 1, output_nvars
      ! 0/1 (space, time) for current variable
      stb = output_is_requested(i, :, :)
      ! indeces in aggregation arrays of current variable
      sti = sum(output_is_requested(:i, :, :), dim=1)
      if (output_var_state(i)) then
        if (stb(1, 2) == 1) output_agg_hydrotope_monthly(sti(1, 2), :, :) = &
          output_agg_hydrotope_monthly(sti(1, 2), :, :) / nd
        if (stb(2, 2) == 1) output_agg_subbasin_monthly(sti(2, 2), :, :) = &
          output_agg_subbasin_monthly(sti(2, 2), :, :) / nd
        if (stb(3, 2) == 1) output_agg_catchment_monthly(sti(3, 2), :, :) = &
          output_agg_catchment_monthly(sti(3, 2), :, :) / nd
      end if
    end do
    ! Increase/reset counters and write out if necessary
    if (output_write_interval == "D") then
      call output_write_monthly
    else if (output_write_interval == "M") then
      call output_write_daily
      call output_write_monthly
      output_store_day_counter = 1
    else if (output_write_interval == "Y") then
      output_store_month_counter = output_store_month_counter + 1
    end if
    output_ndays_per_year = output_ndays_per_year + output_ndays_per_month
    output_ndays_per_month = 0

  end subroutine output_month

  subroutine output_year
    ! Take annual averages and write to csv/bin, if interval is annual.
    integer :: sti(3, 3), stb(3, 3)
    integer i, nd

    ! Take average of annual values as needed
    nd = output_ndays_per_year
    do i = 1, output_nvars
      ! 0/1 (space, time) for current variable
      stb = output_is_requested(i, :, :)
      ! indeces in aggregation arrays of current variable
      sti = sum(output_is_requested(:i, :, :), dim=1)
      if (output_var_state(i)) then
        if (stb(1, 3) == 1) output_agg_hydrotope_annual(sti(1, 3), :, :) = &
          output_agg_hydrotope_annual(sti(1, 3), :, :) / nd
        if (stb(2, 3) == 1) output_agg_subbasin_annual(sti(2, 3), :, :) = &
          output_agg_subbasin_annual(sti(2, 3), :, :) / nd
        if (stb(3, 3) == 1) output_agg_catchment_annual(sti(3, 3), :, :) = &
          output_agg_catchment_annual(sti(3, 3), :, :) / nd
      end if
    end do
    ! Write out before reset counters
    if (output_write_interval == "Y") then
      call output_write_daily
      call output_write_monthly
      output_store_day_counter = 1
      output_store_month_counter = 1
    end if
    call output_write_annual
    output_ndays_per_year = 0
  end subroutine output_year

  subroutine output_write_daily
    ! Write out storage/agg arrays as requested and reset to 0
    integer nd

    ! For varying year and month lengths
    nd = output_store_day_counter
    ! Undo day increase of previous output_day call
    if (output_write_interval /= "D") nd = nd - 1

    call output_write_time(1, output_day_ix(:nd), &
      output_storage_hydrotope(output_hydrotope_requested_ix, :, :nd), &
      output_agg_subbasin_daily(:, :, :nd), &
      output_agg_catchment_daily(:, :, :nd))
    ! Set all storage arrays to 0
    output_storage_hydrotope = 0
    output_storage_subbasin = 0
    output_agg_subbasin_daily = 0
    output_agg_catchment_daily = 0
  end subroutine output_write_daily

  subroutine output_write_monthly
    ! Write out storage/agg monthly arrays CSV output
    call output_write_time(2, output_month_ix, output_agg_hydrotope_monthly, &
      output_agg_subbasin_monthly, output_agg_catchment_monthly)
    ! Set all storage arrays to 0
    output_agg_hydrotope_monthly = 0
    output_agg_subbasin_monthly = 0
    output_agg_catchment_monthly = 0
  end subroutine output_write_monthly

  subroutine output_write_annual
    ! Write out storage/agg annual arrays CSV output
    call output_write_time(3, (/output_month_ix(1)(:4)/), output_agg_hydrotope_annual, &
      output_agg_subbasin_annual, output_agg_catchment_annual)
    ! Set all storage arrays to 0
    output_agg_hydrotope_annual = 0
    output_agg_subbasin_annual = 0
    output_agg_catchment_annual = 0
  end subroutine output_write_annual

  subroutine output_write_time(timeix, time_label, hydrotope, subbasin, catchment)
    ! Abstraction of write out storage/agg arrays
    integer, intent(in) :: timeix
    character(len=*), dimension(:), intent(in) :: time_label
    real, dimension(:, :, :), intent(in) :: hydrotope, subbasin, catchment
    type(output_file) :: fl
    integer i, t, nv, ntm, osix(output_nvars)

    ntm = size(time_label)
    do i = 1, output_nfiles
      fl = output_files(i)
      ! Only continue with files of the timeix frequency
      if (fl%time /= timeix) cycle
      ! Abbreviations
      nv = fl%nvars
      osix = output_storage_index(i, :)
      ! CSV output
      if (fl%format == 1) then
        select case (fl%space)
          case (1)
            call output_write_space_time_csv(fl%file_id, hydrotope(osix(:nv), :, :), &
              output_hydrotope_id, time_label)
          case (2)
            call output_write_space_time_csv(fl%file_id, subbasin(osix(:nv), :, :), &
              output_subbasin_id, time_label)
          case (3)
            call output_write_space_time_csv(fl%file_id, catchment(osix(:nv), :, :), &
              output_catchment_id, time_label)
          case (4)
            call output_write_space_time_csv(fl%file_id, hydrotope(osix(:nv), &
              output_hydrotope_label_ix, :), time_label=time_label, &
              space_label=output_hydrotope_label(output_hydrotope_label_ix))
          case (5)
            call output_write_space_time_csv(fl%file_id, subbasin(osix(:nv), &
              output_subbasin_label_ix, :), time_label=time_label, &
              space_label=output_subbasin_label(output_subbasin_label_ix))
        end select
      ! Binary output
      else if (fl%format == 2) then
        select case (fl%space)
          case (1)
            do t = 1, ntm
              call output_array_to_bin(fl%file_id, hydrotope(osix(:nv), :, t), fl%current_record)
              output_files(i)%current_record = output_files(i)%current_record + 1
            end do
          case (2)
            do t = 1, ntm
              call output_array_to_bin(fl%file_id, subbasin(osix(:nv), :, t), fl%current_record)
              output_files(i)%current_record = output_files(i)%current_record + 1
            end do
          case (3)
            do t = 1, ntm
              call output_array_to_bin(fl%file_id, catchment(osix(:nv), :, t), fl%current_record)
              output_files(i)%current_record = output_files(i)%current_record + 1
            end do
          case (4)
            do t = 1, ntm
              call output_array_to_bin(fl%file_id, hydrotope(osix(:nv), &
                output_hydrotope_label_ix, t), fl%current_record)
              output_files(i)%current_record = output_files(i)%current_record + 1
            end do
          case (5)
            do t = 1, ntm
              call output_array_to_bin(fl%file_id, subbasin(osix(:nv), &
                output_subbasin_label_ix, t), fl%current_record)
              output_files(i)%current_record = output_files(i)%current_record + 1
            end do
        end select
      end if
    end do
  end subroutine output_write_time

  subroutine output_close
    integer fi
    ! Close files
    do fi = 1, output_nfiles
      close(output_files(fi)%file_id)
    end do
    ! TODO are deallocate arrays really necessary?
  end subroutine output_close

  subroutine output_write_space_time_csv(file_id, array, space_ix, time_label, &
      space_label)
    ! Write a 3D (vars, space, time) array with time, space index to tidy csv
    integer, intent(in) :: file_id
    integer, intent(in), optional :: space_ix(:)
    character(len=*), intent(in) :: time_label(:)
    character(len=*), intent(in), optional :: space_label(:)
    real, intent(in) :: array(:, :, :)
    integer it

    if (present(space_ix)) then
      do it = 1, size(time_label)
        call output_array_to_csv(file_id, array(:, :, it), space_ix, time_label(it:it))
      end do
    else if (present(space_label)) then
      do it = 1, size(time_label)
        call output_array_to_csv(file_id, array(:, :, it), &
          time_label=time_label(it:it), space_label=space_label)
      end do
    else
      call log_error("output_write_space_time_csv", &
        "space_ix or space_label have to be given.")
    end if
  end subroutine output_write_space_time_csv

  subroutine output_array_to_csv(file_id, array, space_ix, time_label, space_label)
    ! Write 2D real array to csv in file_id with a space and time index
    ! Both space_ix and time_label must be either the same length as array or of length 1
    ! to repeat over all rows.
    integer, intent(in) :: file_id
    real, intent(in) :: array(:, :)
    integer, intent(in), optional :: space_ix(:)
    character(len=*), intent(in) :: time_label(:)
    character(len=*), intent(in), optional :: space_label(:)
    character(len=identifier_max_length) :: tsfmt
    integer :: sh(2), i, ii, lens, lent

    ! Check input
    if (rank(array) /= 2) call log_error("output_array_to_csv", &
      "output_array_to_csv array is not 2D.")
    sh = shape(array)
    lens = 1
    if (present(space_ix)) then
      lens = size(space_ix)
      if (lens /= sh(2) .and. lens /= 1) call log_error("output_array_to_csv", &
        "space_ix not the right length.")
    else if (present(space_label)) then
      lens = size(space_label)
      if (lens /= sh(2) .and. lens /= 1) &
        call log_error("output_array_to_csv", "space_label not the right length.")
    else
      call log_error("output_array_to_csv", "space_ix or space_label have to be given.")
    end if
    lent = size(time_label)
    if (lent /= sh(2) .and. lent /= 1) &
      call log_error("output_array_to_csv", "time_label not the right length or 1.")
    ! Time, space format
    tsfmt = "(a, ',', "//trim(output_space_index_format)//", ',')"
    do i = 1, sh(2)
      ! Time, space indeces
      if (present(space_ix)) then
        write(file_id, tsfmt, advance="no") time_label(min(i, lent)), &
          space_ix(min(i, lens))
      else
        write(file_id, "(a, ', ', a, ', ')", advance="no") time_label(min(i, lent)), &
          trim(space_label(min(i, lens)))
      end if
      do ii = 1, sh(1)-1
        write(file_id, "("//trim(output_float_format)//", ',')", advance="no") &
          array(ii, i)
      end do
      ! Last column without ,
      write(file_id, "("//trim(output_float_format)//")") array(sh(1), i)
    end do
  end subroutine output_array_to_csv

  subroutine output_array_to_bin(file_id, array, record)
    ! Write array to unformatted (binary) output in file_id
    integer, intent(in) :: file_id, record
    real, intent(in) :: array(:, :)

    write(file_id, rec=record) array
  end subroutine output_array_to_bin

  subroutine output_print_variables
    ! Print registered subbasin/hydrotope output variables
    use utilities, only : string_index
    write(*, "(A)") output_registered_vars(:string_index("", output_registered_vars)-1)
    stop
  end subroutine output_print_variables

end module output