
module input

  use utilities, only : &
    open_file, &
    log_debug, &
    log_info, &
    log_warn, &
    log_error, &
    dp, &
    check_range, &
    check_int_range, &
    path_max_length

  implicit none


  ! Path to config file set via commandline argument
  character(len=path_max_length) :: config_file_path = ""

  integer :: &
    ! Incremental line reading size in characters, must fit max. column name length
    read_line_buffer_size = 128

  integer, parameter :: SEEK_SET = 0, SEEK_CUR = 1, SEEK_END = 2

  ! Input prefix
  character(len=path_max_length), save :: input_dir = "input/"

  ! Input parameters
  logical :: correct_precipitation = .false.

  ! Counter for offset in discharge.csv column, as leap year logic etc is
  ! in outer loops
  integer :: discharge_offset = 1

  namelist / INPUT_PARAMETERS / &
    input_dir, &
    correct_precipitation

#ifdef with_netcdf

  ! INPUT FILES AND VARIABLE NAMES
  ! variables read from ncinfo file
  character(len=path_max_length) :: NC_GRID = "subbasin_climate_grid.csv"
  character(len=path_max_length) :: NC_FNAMES(6) = ''  ! netCDF file names for 6 variables
  character(len=path_max_length) :: NC_VNAMES(6) = ''  ! netCDF variable names for 6 variables
  real :: NC_OFFSET(6) = 0  ! added to value at runtime for unit conversion
  real :: NC_SCALE_FACTOR(6) = 1  ! factor applied at runtime for unit conversion

  character(len = 100) :: NC_LAT_VNAME = "lat"
  character(len = 100) :: NC_LON_VNAME = "lon"
  character(len = 100) :: NC_TIME_VNAME = "time"
  real :: NC_LON_RANGE(2) = (/-180., 180./)
  real :: NC_LAT_RANGE(2) = (/-90., 90./)
  real :: NC_WEIGHT_RANGE(2) = (/0., 1.e6/)

  ! reference year of integer time variable (uninterpreted if < 0)
  integer :: NC_REF_YEAR = - 1
  ! number of days to offset irregardless of NC_REF_YEAR
  integer :: NC_OFFSET_DAYS = 0
  logical :: NC_DEBUG = .true.

  namelist / NC_CLIMATE_PARAMETERS / &
    NC_FNAMES, &
    NC_VNAMES, &
    NC_LAT_VNAME, &
    NC_LON_VNAME, &
    NC_TIME_VNAME, &
    NC_REF_YEAR, &
    NC_OFFSET_DAYS, &
    NC_GRID, &
    NC_OFFSET, &
    NC_SCALE_FACTOR, &
    NC_LON_RANGE, &
    NC_LAT_RANGE, &
    NC_WEIGHT_RANGE, &
    NC_DEBUG
  ! Attributes
  character(len=14) :: nc_attr_names(4) = (/'scale_factor ', 'add_offset   ', &
                                            '_FillValue   ', 'missing_value'/) ! all need to be the same length
  real(dp) :: nc_attr_vals(6, 4) = 0

  ! IDs and counters
  integer :: nc_ids(6), nc_var_ids(6)
  integer :: nc_nlons, nc_nlats, nc_nrecs
  integer :: nc_xmin = 10000000
  integer :: nc_xmax = 0
  integer :: nc_ymin = 10000000
  integer :: nc_ymax = 0
  integer :: nc_nx, nc_ny
  ! current record number, ie. day counter from begining of file (reset in input_nc_find_time)
  integer :: nc_nday = 1
  ! number of rows in ncinfo file, set later
  integer :: nc_nrows = 0
  ! max cells per subbasin
  integer :: nc_mxc = 0
  ! if humidity is not provided hum = -999
  logical :: nc_blank(6) = .false.

  ! allocated with nc_mxc in readinfo
  ! weights read from input file
    real(dp), dimension(:, :), allocatable :: nc_weight_sb
  ! lon/lat read from input file
    real(dp), dimension(:, :), allocatable :: nc_lon_sb, nc_lat_sb
  ! indecies of lon / lat in nc files
  integer, dimension(:, :), allocatable :: nc_x_sb, nc_y_sb

  ! allocatable data arrays
  ! sum of weights per subbasin
    real(dp), dimension(:), allocatable :: nc_wsum_sb
  ! number of cells per subbasin
  integer, dimension(:), allocatable :: nc_ncells_sb
  ! to store lat / nc_lons
    real(dp), dimension(:), allocatable :: nc_lons, nc_lats
  ! time index
  integer, dimension(:), allocatable :: nc_time_ix
  ! 2D array for data of one day
    real(dp), dimension(:, :), allocatable :: nc_var_in

#endif

contains

  subroutine input_initialise

    read(get_config_fid(), INPUT_PARAMETERS)

  end subroutine input_initialise

  integer function input_open_file(file_name) result(file_id)
    ! Abstract input file opener (ascii only)
    character(len=*), intent(in) :: file_name

    file_id = open_file(trim(input_dir)//"/"//trim(file_name), 'r')
  end function input_open_file

#ifdef with_netcdf

  subroutine input_nc_read_climate(flu, humi, mb, ra, subp, tmn, tmx, tx)
    real(dp), dimension(:), intent(in) :: flu
    real(dp), dimension(:), intent(out) :: humi
    integer, intent(in) :: mb
    real(dp), dimension(:), intent(out) :: ra
    real(dp), dimension(:), intent(out) :: subp
    real(dp), dimension(:), intent(out) :: tmn
    real(dp), dimension(:), intent(out) :: tmx
    real(dp), dimension(:), intent(out) :: tx
    ! counters need for day length
    integer j
    ! reals need for day length
    real(dp) :: mm(6) = 0

    call input_nc_weighted_mean(1, tx, mb)
    call input_nc_weighted_mean(2, tmn, mb)
    call input_nc_weighted_mean(3, tmx, mb)
    call input_nc_weighted_mean(4, subp, mb)
    if (nc_blank(5)) then
      ra = 1.
    else
      call input_nc_weighted_mean(5, ra, mb)
    end if
    if (nc_blank(6)) then
      humi = -999.
    else
      call input_nc_weighted_mean(6, humi, mb)
    end if

    if (NC_DEBUG) then
      mm = 0
      do j = 1, mb
        mm(1) = mm(1) + tx(j) * flu(j)
        mm(2) = mm(2) + tmn(j) * flu(j)
        mm(3) = mm(3) + tmx(j) * flu(j)
        mm(4) = mm(4) + subp(j) * flu(j)
        mm(5) = mm(5) + ra(j) * flu(j)
        mm(6) = mm(6) + humi(j) * flu(j)
      enddo
      call log_debug("input_nc_read_climate", &
        'ncdf climate (tmean, tmin, tmax, prec, rad, hum) avarage on day:', &
        int=nc_nday, reals=mm)
    endif

    ! update current record number
    nc_nday = nc_nday + 1

  end subroutine input_nc_read_climate

  subroutine input_nc_weighted_mean(ivar, v_sb, mb)
    use netcdf, only : nf90_get_var
    integer, intent(in) :: mb
    ! variable index
    integer, intent (in) :: ivar
    real(dp), intent (out) :: v_sb(:) ! the subbasin variable to be filled
    integer :: start(3), count(3) ! to supply to nc_get_var
    ! looping ints
    integer sb, i

    ! Read 1 record of NLONS*NLATS*NLVLS values, starting at the beginning
    ! of the record (the (1, 1, rec) element in the netCDF file).
    count = (/ nc_nx, nc_ny, 1 /)
    start = (/ nc_xmin, nc_ymin, nc_nday /)

    ! read one day/record
    call input_nc_check_error( nf90_get_var(nc_ids(ivar), nc_var_ids(ivar), &
                                            nc_var_in, start, count) )
    ! apply offset and scale_factor
    nc_var_in = nc_var_in * nc_attr_vals(ivar, 1) * NC_SCALE_FACTOR(ivar) &
      + nc_attr_vals(ivar, 2) + NC_OFFSET(ivar)

    ! interpolate and fill v_sb
    v_sb = 0
    do sb = 1, mb
      do i = 1, nc_ncells_sb(sb)
        v_sb(sb) = v_sb(sb) + nc_var_in(nc_x_sb(sb, i), nc_y_sb(sb, i)) * nc_weight_sb(sb, i)
      end do
    enddo

  end subroutine input_nc_weighted_mean

  subroutine input_nc_initialise(iyr, mb)
    use netcdf, only : nf90_open, nf90_nowrite, nf90_inq_varid

    integer, intent(in) :: iyr
    integer, intent(in) :: mb
    integer :: i
    logical :: fexists
    character(len=path_max_length) :: fpaths(6)

    call log_info('input_nc_initialise', 'Opening netCDF files...')

    ! allocate arrays with subbasin count
    allocate(nc_wsum_sb(mb))
    nc_wsum_sb = 0
    allocate(nc_ncells_sb(mb))
    nc_ncells_sb = 0

    ! read namelist parameters
    read(get_config_fid(), nml=NC_CLIMATE_PARAMETERS)
    ! open nc files, get var ids and compare the dimensions of first file (0), check if others are the same (1)
    do i = 1, 6
      ! check if file exists in input_dir, otherwise interpret as absolute path
      INQUIRE(file = trim(input_dir) // trim(NC_FNAMES(i)), EXIST = fexists)
      if (fexists) then
        fpaths(i) = trim(input_dir) // trim(NC_FNAMES(i))
      else
        fpaths(i) = trim(NC_FNAMES(i))
      endif
      if (NC_DEBUG) call log_debug("input_nc_initialise", trim(fpaths(i))//' ('//trim(NC_VNAMES(i))//')')
      if (i >= 5 .and. trim(NC_FNAMES(i)) == '') then
        nc_blank(i) = .TRUE.
        call log_debug("input_nc_initialise", "No nc file provided for variable "//NC_FNAMES(i), int=i)
        cycle
      end if
      ! open file, get fileID, check dims and attributes
      call input_nc_check_error( nf90_open(trim(fpaths(i)), nf90_nowrite, nc_ids(i)) )
      call input_nc_check_error( nf90_inq_varid(nc_ids(i), trim(NC_VNAMES(i)), nc_var_ids(i)) )
      call input_nc_check_dims(nc_ids(i), i /= 1, trim(NC_VNAMES(i)))
      call input_nc_check_attr(i)
    enddo

    ! check ncgrid.dat, count rows (nc_nrows) and max cell per subbasins (nc_mxc)
    call input_nc_check_grid(mb)

    ! allocated arrays with nc_mxc
    allocate(nc_weight_sb(mb, nc_mxc))
    nc_weight_sb = 0
    allocate(nc_lon_sb(mb, nc_mxc))
    nc_lon_sb = 0
    allocate(nc_lat_sb(mb, nc_mxc))
    nc_lat_sb = 0
    allocate(nc_x_sb(mb, nc_mxc))
    nc_x_sb = 0
    allocate(nc_y_sb(mb, nc_mxc))
    nc_y_sb = 0

    ! read ncinfo.dat
    call input_nc_read_grid(mb)

    ! convert lon/lat to x/y indecies
    call input_nc_convert_coordinates(mb)
    ! find time offset
    if (NC_REF_YEAR >= 0) then
        call input_nc_find_time(iyr) ! resets nc_nday
    else if (NC_OFFSET_DAYS /= 0) then
        call input_nc_offset_time() ! adds to nc_nday
    endif

    ! allocate nc_var_in to read in
    allocate(nc_var_in(nc_nx, nc_ny))
    nc_var_in = 0

    ! check for missing values
    if (NC_DEBUG) call input_nc_check_missing(mb)

  end subroutine input_nc_initialise

  subroutine input_nc_check_dims(ncid, icheck, vname)
    use netcdf, only : &
      nf90_inq_dimid, nf90_inquire_dimension, nf90_inq_varid, nf90_get_var

    ! file and variable id to read from
    integer, intent (in) :: ncid
    ! check if same as already set values
    logical, intent (in) :: icheck
    ! name of variable being checked
    character(len=*), intent(in) :: vname
    ! dimension ids to read length
    integer :: lon_dimid, lat_dimid, rec_dimid
    ! variable ids to read time, lon/lat values
    integer :: lon_varid, lat_varid, time_varid
    integer :: nlo, nla, nre

    ! Get var ids of lon/lat and time
    call input_nc_check_error( nf90_inq_dimid(ncid, trim(NC_LON_VNAME), lon_dimid) )
    call input_nc_check_error( nf90_inquire_dimension(ncid, lon_dimid, len = nlo) )
    call input_nc_check_error( nf90_inq_dimid(ncid, trim(NC_LAT_VNAME), lat_dimid) )
    call input_nc_check_error( nf90_inquire_dimension(ncid, lat_dimid, len = nla) )
    call input_nc_check_error( nf90_inq_dimid(ncid, trim(NC_TIME_VNAME), rec_dimid) )
    call input_nc_check_error( nf90_inquire_dimension(ncid, rec_dimid, len = nre) )

    ! report dimensions
    if (NC_DEBUG) call log_debug("input_nc_check_dims", &
      'Number of nc_lons, nc_lats, timesteps:', ints=(/nlo, nla, nre/))

    ! check or assign
    if (icheck) then
      if (nlo /= nc_nlons) then
        call log_error("input_nc_check_dims", vname//' has different number of nc_lons', &
          ints=(/nlo, nc_nlons/))
      endif
      if (nla /= nc_nlats) then
        call log_error("input_nc_check_dims", vname//' has different number of nc_lats', &
          ints=(/nla, nc_nlats/))
      endif
      if (nre /= nc_nrecs) then
        call log_error("input_nc_check_dims", vname//' has different number of time steps', &
          ints=(/nre, nc_nrecs/))
      endif
    else
      ! assign
      nc_nlons = nlo
      nc_nlats = nla
      nc_nrecs = nre
      ! allocate nc_lons, nc_lats, time
      allocate(nc_lons(nc_nlons))
      allocate(nc_lats(nc_nlats))
      allocate(nc_time_ix(nc_nrecs))
      ! read time, nc_lons and nc_lats
      call input_nc_check_error( nf90_inq_varid(ncid, trim(NC_LON_VNAME), lon_varid) )
      call input_nc_check_error( nf90_get_var(ncid, lon_varid, nc_lons) )
      call input_nc_check_error( nf90_inq_varid(ncid, trim(NC_LAT_VNAME), lat_varid) )
      call input_nc_check_error( nf90_get_var(ncid, lat_varid, nc_lats) )
      call input_nc_check_error( nf90_inq_varid(ncid, trim(NC_TIME_VNAME), time_varid) )
      call input_nc_check_error( nf90_get_var(ncid, time_varid, nc_time_ix) )
    endif

  end subroutine input_nc_check_dims

  subroutine input_nc_check_attr(ivar)
    use netcdf, only : nf90_get_att, NF90_NOERR
    ! variable counter
    integer, intent (in) :: ivar
    integer :: err, i
    real(dp) :: attval

    if (NC_DEBUG) then
      call log_debug("input_nc_check_attr", &
        "SWIM applies offset to "//trim(NC_VNAMES(ivar))//":", real=NC_OFFSET(ivar))
      call log_debug("input_nc_check_attr", &
        "SWIM applies scale factor to "//trim(NC_VNAMES(ivar))//":", &
        real=NC_SCALE_FACTOR(ivar))
    end if
    ! check for variable attributes
    do i = 1, 4
      err = nf90_get_att(nc_ids(ivar), nc_var_ids(ivar), trim(nc_attr_names(i)) , attval)
      if (err == NF90_NOERR) then
        nc_attr_vals(ivar, i) = attval
      ! Attribute default values if not set in addition to initialisation (=0)
      elseif (i == 1) then ! scale_factor
        nc_attr_vals(ivar, i) = 1
      elseif (i == 4) then ! missing value
        nc_attr_vals(ivar, i) = nc_attr_vals(ivar, 3)
      endif
      if (NC_DEBUG .and. err == NF90_NOERR) then
        call log_debug("input_nc_check_attr", trim(NC_VNAMES(ivar))//' has a '// &
          nc_attr_names(i)//' attribute of:', real=nc_attr_vals(ivar, i))
      elseif (NC_DEBUG) then
        call log_debug("input_nc_check_attr", trim(NC_VNAMES(ivar))//' has a '// &
          nc_attr_names(i)//' default value of:', real=nc_attr_vals(ivar, i))
      endif
    enddo

  endsubroutine input_nc_check_attr

  subroutine input_nc_check_grid(mb)
    integer, intent(in) :: mb
    integer :: fui, sb
    integer :: lastsb = 0
    integer :: imxc = 0

    ! check ncinfo file for length and maxcells
    fui = input_open_file(NC_GRID)
    read(fui, *) ! header line
    do
      read(fui, *, END=10) sb
      nc_nrows = nc_nrows + 1
      ! count cells per subbasin and update max counter
      if (sb /= lastsb) imxc = 0
      lastsb = sb
      imxc = imxc + 1
      nc_mxc = max(nc_mxc, imxc)
      ! check subbasin count
      if (sb > mb) then
        call log_warn("input_nc_check_grid", "SubbasinID in (line) larger than number of subbasins:", &
          i1=sb, int=nc_nrows)
      endif
    enddo
    10 close(fui)
    ! make sure lastsb==mb
    if (lastsb /= mb) then
      call log_warn("input_nc_check_grid", 'Last subbasinID in '//trim(NC_GRID)// &
        ' does not match the last subbasinID', i1=lastsb, int=mb)
    endif
    ! report file length and maxcells per subbasin
    if (NC_DEBUG) call log_debug("input_nc_check_grid", 'Lines in ncgrid file:', int=nc_nrows)
    if (NC_DEBUG) call log_debug("input_nc_check_grid", 'Maximum cells per subbasin:', int=nc_mxc)
  end subroutine input_nc_check_grid


  subroutine input_nc_read_grid(mb)
    integer, intent(in) :: mb
    integer :: fui, i, ii, sb
    real(dp) :: ilon, ilat, iweight

    ! open weights and coordinates file
    fui = input_open_file(NC_GRID)
    read(fui, *) ! header line
    do i = 1, nc_nrows
      ! read line
      read(fui, *) sb, ilon, ilat, iweight
      ! count cells
      nc_ncells_sb(sb) = nc_ncells_sb(sb) + 1
      ii = nc_ncells_sb(sb)
      ! check input and assign
      call check_range((/ilon/), 'longitude', NC_LON_RANGE, sb)
      nc_lon_sb(sb, ii) = ilon
      call check_range((/ilat/), 'latitude', NC_LAT_RANGE, sb)
      nc_lat_sb(sb, ii) = ilat
      call check_range((/iweight/), 'weight', NC_WEIGHT_RANGE, sb)
      nc_weight_sb(sb, ii) = iweight
      ! weights sum
      nc_wsum_sb(sb) = nc_wsum_sb(sb) + nc_weight_sb(sb, ii)
    end do
    ! close info file again
    close(fui)

    ! adjust/make sure weights add up to 1
    do sb = 1, mb
      do ii = 1, nc_ncells_sb(sb)
        nc_weight_sb(sb, ii) = nc_weight_sb(sb, ii) / nc_wsum_sb(sb)
      enddo
    enddo

  end subroutine input_nc_read_grid

  subroutine input_nc_convert_coordinates(mb)
    integer, intent(in) :: mb
    integer sb, c, i
    ! loop over subbasins, then cells, then lon/lat stopping at the matching one
    do sb = 1, mb
      do c = 1, nc_ncells_sb(sb)
        do i = 1, nc_nlons
          if (nc_lons(i) == nc_lon_sb(sb, c)) then
            nc_x_sb(sb, c) = i
            nc_xmin = min(nc_xmin, i)
            nc_xmax = max(nc_xmax, i)
            exit
          endif
          ! if it hasn't found one and got to the end, stop and warn
          if (i == nc_nlons) then
            call log_error("input_nc_convert_coordinates", 'Longitude (at subbasin, cell) not in ncdf lon!', &
              real=nc_lon_sb(sb, c), i1=sb, i2=c)
          endif
        end do
        do i = 1, nc_nlats
          if (nc_lats(i) == nc_lat_sb(sb, c)) then
            nc_y_sb(sb, c) = i
            nc_ymin = min(nc_ymin, i)
            nc_ymax = max(nc_ymax, i)
            exit
          endif
          ! if it hasn't found one ang got to the end, stop and warn
          if (i == nc_nlats) then
            call log_error("input_nc_convert_coordinates", 'Latitude (at subbasin, cell) not in ncdf lat!', &
              real=nc_lon_sb(sb, c), i1=sb, i2=c)
          endif
        end do
      enddo
    enddo
    ! adjust x/y indicies for minimal window
    nc_nx = nc_xmax - (nc_xmin - 1)
    nc_ny = nc_ymax - (nc_ymin - 1)
    do sb = 1, mb
      do c = 1, nc_ncells_sb(sb)
        nc_x_sb(sb, c) = nc_x_sb(sb, c) - (nc_xmin - 1)
        nc_y_sb(sb, c) = nc_y_sb(sb, c) - (nc_ymin - 1)
      enddo
    enddo
    if (NC_DEBUG) call log_debug("input_nc_convert_coordinates", &
      'Minimum window size (n nc_lons, n nc_lats):', i1=nc_nx, i2=nc_ny)
  end subroutine input_nc_convert_coordinates

  subroutine input_nc_find_time(iyr)
    integer, intent(in) :: iyr
    integer ndays_offset, y, d
    logical :: time_not_found = .true.

    ! should be replaced by common switch from cod - file
    integer :: ileap = 1

    ndays_offset = (iyr - NC_REF_YEAR) * 365 - NC_OFFSET_DAYS
    ! add leap years
    if (ileap == 1) then
      do y = NC_REF_YEAR, iyr - 1
        if ((MOD(y, 4) == 0 .AND. MOD(y, 100) > 0) .OR. MOD(y, 400) == 0) then
          ndays_offset = ndays_offset + 1
        endif
      enddo
    endif

    if (NC_DEBUG) call log_debug("input_nc_find_time", &
      'Searching for netCDF time index', int=ndays_offset)
    do d = 1, nc_nrecs
      if (ndays_offset == nc_time_ix(d)) then
        nc_nday = d
        time_not_found = .false.
        if (NC_DEBUG) call log_debug("input_nc_find_time", 'Found time index at', int=d)
        exit ! break from do loop
      endif
    enddo
    ! stop if time index not found
    if (time_not_found) then
      call log_error("input_nc_find_time", &
        'Cant find time index in climate netCDF file (start, ref. year):', i1=iyr, &
        i2=NC_REF_YEAR, int=ndays_offset)
    endif
  end subroutine input_nc_find_time

  subroutine input_nc_offset_time()
    nc_nday = nc_nday + NC_OFFSET_DAYS
    if (nc_nday <= 0 .OR. nc_nday > nc_nrecs) then
      call log_error("input_nc_offset_time", &
        'Time NC_OFFSET_DAYS leads to out of range index ', int=nc_nday)
    endif
    if (NC_DEBUG) call log_error("input_nc_offset_time", &
      'Offsetting time index (by) to:', i1=NC_OFFSET_DAYS, int=nc_nday)
  end subroutine input_nc_offset_time

  subroutine input_nc_check_missing(mb)
    use netcdf, only : nf90_get_var
    integer, intent(in) :: mb
    integer :: sb, i, ii, nmissing
    real(dp) :: v

    call log_debug("input_nc_check_missing", 'Checking for missing values...')
    nmissing = 0
    do i = 1, 6
      if (nc_blank(i)) cycle
      ! read one day/record
      call input_nc_check_error( nf90_get_var(nc_ids(i), nc_var_ids(i), nc_var_in, (/ nc_xmin, nc_ymin, nc_nday /), (/ nc_nx, nc_ny, 1 /)) )
      ! loop over nc_var_in again and check for missing or fill values
      do sb = 1, mb
        do ii = 1, nc_ncells_sb(sb)
          v = nc_var_in(nc_x_sb(sb, ii), nc_y_sb(sb, ii))
          if (v == nc_attr_vals(i, 3) .or. v == nc_attr_vals(i, 4)) then
            call log_warn("input_nc_check_missing", trim(NC_FNAMES(i))//' ('// &
              trim(NC_VNAMES(i))//') missing in (subbasin):', i1=sb, &
              reals=(/nc_lon_sb(sb, ii), nc_lat_sb(sb, ii)/))
            nmissing = nmissing + 1
          endif
        enddo
      enddo
    enddo
    if (nmissing > 0) then
      call log_error("input_nc_check_missing", 'Encountered N missing values:', &
        int=nmissing)
    endif
  endsubroutine input_nc_check_missing

  subroutine input_nc_close
    use netcdf, only : nf90_close
    integer :: i
    ! deallocate arrays and close open files
    do i = 1, 6
      if (.not. nc_blank(i)) call input_nc_check_error( nf90_close(nc_ids(i)) )
      !call input_nc_check_error( nf90_close(nc_ids(i)) )
    enddo

    deallocate(nc_lons)
    deallocate(nc_lats)
    deallocate(nc_var_in)

  endsubroutine input_nc_close

  subroutine input_nc_check_error(status)
    use netcdf, only : nf90_noerr, nf90_strerror

    integer, intent ( in) :: status

    if (status /= nf90_noerr) then
      call log_error('input_nc_check_error', trim(nf90_strerror(status)))
    end if
  end subroutine input_nc_check_error

#endif

  integer function get_config_fid()
    integer unit

    unit = -1
    inquire(file=config_file_path, number=unit)
    if (unit >= 0) then
      close(unit)
    end if
    get_config_fid = open_file(config_file_path)
  end function get_config_fid

  function input_count_rows(funit, header) result(nlines)
    ! Counts the number of rows in input file with given file unit
    ! The file must already be open!!!
    ! header (bool, default=.true.) : If false, dont count first line.
    ! file unit
    integer, intent(in) :: funit
    ! true if input file contains a header
    logical, intent(in), optional :: header
    ! counter
    integer :: nlines
    ! file status
    integer :: status
    ! dummy for reading
    character(len=read_line_buffer_size) :: a
    character(len=path_max_length) fname

    status = 0
    nlines = -1
    if (present(header)) then
      if (.not. header) nlines = nlines + 1
    end if
    ! rewind file in order to be sure to read the file from the first line
    rewind(funit)

    ! count number of rows in file: fname
    do
      read(funit, *, IOSTAT=status) a
      if (status > 0) then  ! read error (unlikely since we are only skipping lines)
        inquire(unit=funit, name=fname)
        call log_error("input_count_rows", 'Unable to count lines of file '//trim(fname))
      else if (status == -1 .or. len_trim(a) == 0) then  ! end of file or empty line
        EXIT
      end if
      nlines = nlines + 1
    end do

    rewind(funit)
  end function input_count_rows

  elemental function precipitation_correction(precip, tmean) result(corrected_p)
    ! Description needed
    real(dp), intent(in) :: precip, tmean
    real(dp) :: corrected_p
  
    real(dp) :: k

    if (precip .gt. 0.1) then
      if (tmean .lt. 0.) then
        k = exp(0.0459 + 0.2367 * 2 + 0.018 * tmean - 0.0154 * tmean)
      else
        k = exp(-0.0010 * log(precip) - 0.0122 * 2 * log(precip) + 0.0343 * 2 + 0.0077 + 0.01)
      endif
      corrected_p = k * precip
    else
      corrected_p = 0.
    endif

  end function precipitation_correction

  subroutine read_logical_column(file_id, column, array, default, index, skip)
    ! Read column of logical values to fill a 1d array
    !
    ! Optional arguments:
    ! -------------------
    ! column : Name of column to read
    ! default : Value to fill array if column not found, raises error otherwise
    ! index : Read column by position/index, since the header is not read
    !   the column is read from the current line onwards
    ! skip : Skip rows after header (column parsed) or current position (index)
    integer, intent(in) :: file_id
    character(len=*), intent(in), optional :: column
    logical, dimension(:), intent(out) :: array
    logical, intent(in), optional :: default
    integer, intent(in), optional :: index, skip
    integer i, column_ix
    character(len=read_line_buffer_size) str_value
    integer temp, ios

    if (present(column)) then
      ! Find column index
      column_ix = header_column_index(column, file_id)
    else if (present(index)) then
      column_ix = index
    else
      call log_error("read_logical_column", "column or index argument required.")
    end if
    ! Skip lines after the header
    if (present(skip)) call move_lines(file_id, skip)
    ! Set default or error if not found
    if (column_ix == 0 .and. present(default)) then
      array = default
      return
    else if (column_ix == 0) then
      call input_error_column_not_found(column, file_id)
    end if
    ! Fill array
    do i = 1, size(array)
      str_value = read_csv_item(column_ix, file_id)
      read(str_value, *, iostat=ios) temp
      if (ios /= 0 .or. .not. (temp == 1 .or. temp == 0)) &
        call input_type_conversion_error(file_id, str_value, "0/1 logical", i+1, &
          column_ix)
      array(i) = temp == 1
    end do
  end subroutine read_logical_column

  subroutine read_real_column(file_id, column, array, default, index, skip, range, closed)
    ! Read column of real values to fill a 1d array
    !
    ! Optional arguments:
    ! -------------------
    ! column : Name of column to read
    ! default : Value to fill array if column not found, raises error otherwise
    ! index : Read column by position/index, since the header is not read
    !   the column is read from the current line onwards
    ! skip : Skip rows after header (column parsed) or current position (index)
    ! range(2) : Check (/lower, upper/) range, column name is required
    ! closed : Open or closed range, "l"/"u"/"n" lower/upper/neiter, default both
    integer, intent(in) :: file_id
    character(len=*), intent(in), optional :: column
    real(dp), dimension(:), intent(out) :: array
    real(dp), intent(in), optional :: default, range(2)
    integer, intent(in), optional :: index, skip
    character(len=1), intent(in), optional :: closed
    integer i, column_ix, ios
    character(len=read_line_buffer_size) str_value
    character(len=1) clsd

    if (present(column)) then
      ! Find column index
      column_ix = header_column_index(column, file_id)
    else if (present(index)) then
      column_ix = index
    else
      call log_error("read_real_column", "column or index argument required.")
    end if
    ! Skip lines after the header
    if (present(skip)) call move_lines(file_id, skip)
    ! Set default or error if not found
    if (column_ix == 0 .and. present(default)) then
      call log_debug("read_real_column", "using default for column "//trim(column), &
        real=default)
      array = default
      return
    else if (column_ix == 0) then
      call input_error_column_not_found(column, file_id)
    end if
    ! Fill array
    do i = 1, size(array)
      str_value = read_csv_item(column_ix, file_id)
      read(str_value, *, iostat=ios) array(i)
      if (ios /= 0) &
        call input_type_conversion_error(file_id, str_value, "real", i+1, column_ix)
    end do
    ! Check ranges
    if (present(range) .and. present(column)) then
      clsd = "b"
      if (present(closed)) clsd = closed
      call check_range(array, trim(column), range, closed=clsd)
    end if
  end subroutine read_real_column

  subroutine read_integer_column(file_id, column, array, default, index, skip, range)
    ! Read column of integer values to fill a 1d array
    !
    ! Optional arguments: see read_real_column
    integer, intent(in) :: file_id
    character(len=*), intent(in), optional :: column
    integer, dimension(:), intent(out) :: array
    integer, intent(in), optional :: default, range(2)
    integer, intent(in), optional :: index, skip
    integer i, column_ix, ios
    character(len=read_line_buffer_size) str_value

    if (present(column)) then
      ! Find column index
      column_ix = header_column_index(column, file_id)
    else if (present(index)) then
      column_ix = index
    else
      call log_error("read_integer_column", "column or index argument required.")
    end if
    ! Skip lines after the header
    if (present(skip)) call move_lines(file_id, skip)
    ! Set default or error if not found
    if (column_ix == 0 .and. present(default)) then
      call log_debug("read_integer_column", "using default for column "//trim(column), &
        int=default)
      array = default
      return
    else if (column_ix == 0) then
      call input_error_column_not_found(column, file_id)
    end if
    ! Fill array
    do i = 1, size(array)
      str_value = read_csv_item(column_ix, file_id)
      ! Accept fewer items, as newlines at the end of a csv file will break the reader
      if (len(trim(adjustl(str_value))) .eq. 0) then
        return
      end if
      read(str_value, *, iostat=ios) array(i)
      if (ios /= 0) &
        call input_type_conversion_error(file_id, str_value, "integer", i+1, column_ix)
    end do
    ! Check ranges
    if (present(range) .and. present(column)) then
      call check_int_range(array, trim(column), range)
    end if
  end subroutine read_integer_column

  subroutine read_string_column(file_id, column, array, default, index, skip)
    ! Read column of strings to fill a 1d array
    !
    ! Optional arguments: see read_real_column
    integer, intent(in) :: file_id
    character(len=*), intent(in), optional :: column
    character(len=*), dimension(:), intent(out) :: array
    character(len=*), intent(in), optional :: default
    integer, intent(in), optional :: index, skip
    integer i, column_ix

    if (present(column)) then
      ! Find column index
      column_ix = header_column_index(column, file_id)
    else if (present(index)) then
      column_ix = index
    else
      call log_error("read_string_column", "column or index argument required.")
    end if
    ! Skip lines after the header
    if (present(skip)) call move_lines(file_id, skip)
    ! Set default or error if not found
    if (column_ix == 0 .and. present(default)) then
      call log_debug("read_integer_column", "using default for column "// &
        trim(column)//": "//trim(default))
      array = default
      return
    else if (column_ix == 0) then
      call input_error_column_not_found(column, file_id)
    end if
    ! Fill array
    do i = 1, size(array)
      array(i) = read_csv_item(column_ix, file_id)
    end do
  end subroutine read_string_column

  logical function has_column(column, file_id)
    character(len=*), intent(in) :: column
    integer, intent(in) :: file_id
    has_column = header_column_index(trim(adjustl(column)), file_id) >0
  end function has_column

  integer function header_column_index(column, file_id, required)
    ! Return index of column in first line of file_id, return 0 if not found
    character(len=*), intent(in) :: column
    integer, intent(in) :: file_id
    logical, intent(in), optional :: required
    integer iostat, ie
    character(len=read_line_buffer_size) buffer
    character(len=read_line_buffer_size*2) scanpart, current_col
    logical column_not_found

    rewind(file_id)
    scanpart = ''
    column_not_found = .True.
    header_column_index = 0
    iostat = 0
    do while (column_not_found .and. iostat == 0)
      read(file_id, "(a)", advance="no", iostat=iostat) buffer
      scanpart = trim(adjustl(scanpart))//buffer
      ie = 1
      do while (column_not_found .and. ie > 0)
        ie = scan(scanpart, ",")
        if (ie > 0) then
          current_col = trim(adjustl(scanpart(:ie-1)))
          scanpart = scanpart(ie+1:)
          header_column_index = header_column_index + 1
          column_not_found = trim(adjustl(current_col)) /= trim(adjustl(column))
        else if (iostat < 0) then  ! last column
          current_col = scanpart
          header_column_index = header_column_index + 1
          column_not_found = trim(adjustl(current_col)) /= trim(adjustl(column))
        end if
      end do
    end do
    if (column_not_found) then
      header_column_index = 0
      if (present(required)) then
        if (required) call input_error_column_not_found(column, file_id)
      end if
    end if
    ! Make sure we move to next line if not already
    if (iostat == 0) call move_lines(file_id, 1)
  end function header_column_index

  subroutine move_lines(file_id, n)
    ! Move n lines forwards or backwards (if negative) in formatted sequential file
    integer, intent(in) :: file_id, n
    integer i
    if (n > 0) then  ! move forwards
      do i = 1, n
        read(file_id, '()')
      end do
    else if (n < 0) then  ! move backwards
      do i = 1, -1*n
        backspace file_id
      end do
    end if
  end subroutine move_lines

  function read_csv_item(index, file_id) result(column)
    ! Return the index'th item from the next line of file_id and move to next line
    ! Raise error if number of items not found
    integer, intent(in) :: index, file_id
    character(len=read_line_buffer_size) :: column
    character(len=read_line_buffer_size) buffer
    character(len=read_line_buffer_size*2) scanpart
    logical column_not_found
    integer ic, iostat, ie

    scanpart = ''
    column_not_found = .True.
    ic = 0
    iostat = 0
    do while (column_not_found .and. iostat == 0)
      read(file_id, "(a)", advance="no", iostat=iostat) buffer
      scanpart = trim(adjustl(scanpart))//buffer
      ie = 1
      do while (column_not_found .and. ie > 0)
        ie = scan(scanpart, ",")
        if (ie > 0) then  ! comma found
          ic = ic + 1
          column_not_found = ic /= index
          if (.not. column_not_found) column = trim(adjustl(scanpart(:ie-1)))
          scanpart = scanpart(ie+1:)
        else if (iostat < 0 .and. ic+1 == index) then  ! last item
          column = trim(adjustl(scanpart))
          column_not_found = .False.
        end if
      end do
    end do
    if (column_not_found) then
      ! TODO error if not found, e.g. because end of file or line
    else
      ! Move to next line if not already
      if (iostat == 0) read(file_id, "(a)")
    end if
  end function read_csv_item

  subroutine input_error_column_not_found(column, file_id)
    character(len=*), intent(in) :: column
    integer, intent(in) ::  file_id
    character(len=path_max_length) :: fname

    inquire(file_id, name=fname)
    call log_error("input_error_column_not_found", &
      "Required column "//trim(column)//" not found in "//trim(fname))
  end subroutine input_error_column_not_found

  subroutine input_type_conversion_error(file_id, str, type, line_index, column_index)
    character(len=*), intent(in) :: str, type
    integer, intent(in) ::  file_id, line_index, column_index
    character(len=path_max_length) :: fname

    inquire(file_id, name=fname)
    call log_error("input_type_conversion_error", "Failed to convert "//trim(str)// &
      " to "//trim(type)//" of "//trim(fname)//" in (line, column)", i1=line_index, i2=column_index)
  end subroutine input_type_conversion_error

  subroutine print_splash(version)
    use utilities, only: colourise
    character(len=*), intent(in) :: version
    integer :: nlsp, nrsp
    character(len=12) :: l1
    character(len=14) :: l2, l4
    character(len=16) :: l3
    l1 = colourise("~~~", 'yellow')
    l2 = colourise("~~~~~", 'green')
    l3 = colourise("~~~~~~~", 'cyan')
    l4 = colourise("~~~~~", 'blue')

    nlsp = floor((49 - len(version)) / 2.)
    nrsp = ceiling((49 - len(version)) / 2.)
    write(*, "(A)") "", &
      "  + . . . . . . . . . . . . . . . . . . . . . . . . . . +", &
      "  .                                                     .", &
      "  .      .      //////// \\    \\      // ||  ||    ||  .", &
      "  .     / \     \\        \\    \\    //  ||  ||\\//||  .", &
      "  .    /"//l1//"\     \\\\\\\   \\  //\\  //   ||  || \/ ||  .", &
      "  .   /"//l2//"\         //    \\//  \\//    ||  ||    ||  .", &
      "  .  (" //l3// ")  ///////      \\    \\     ||  ||    ||  .", &
      "  .   \"//l4//"/   --------------------------------------  .", &
      "  .    -----    Soil   and   Water  Integrated   Model  .", &
      "  .                                                     .", &
      "  ."//repeat(" ", nlsp)//"* "//version//" *"//repeat(" ", nrsp)//".", &
      "  + . . . . . . . . . . . . . . . . . . . . . . . . . . +", ""
  end subroutine print_splash

  subroutine print_help(errormsg)
    use utilities, only: log_error
    character(len=*), optional :: errormsg

    write(*, "(A)") &
      "usage: swim [options] [parameter-nml]", &
      "", &
      "  parameter-nml               path to the parameter namelist", &
      "", &
      "options:", &
      "  -h, --help                  show this help message", &
      "  -v, --version               print version", &
      "  -d, --defaults [module]     print default parameters", &
      "  -o, --output-variables      print output variables (requires parameter-nml)"
    ! End program either with error or silently
    if (present(errormsg)) then
      write(*, *)  ! empty line btw help and error
      call log_error("print_help", errormsg)
    end if
    stop
  end subroutine print_help

  subroutine parse_commandline_arguments(version, print_output_variables, print_defaults)
    use utilities, only : master_log, log_error_level

    character(len=*), intent(in) :: version
    logical, intent(out) :: print_output_variables
    character(len=*), intent(out) :: print_defaults
    character(len=path_max_length) :: arg
    integer i, iarg

    ! Number of commandline arguments
    iarg = iargc()
    ! without arguments only show help
    if (iarg == 0) call print_help("parameter-nml is missing")
    print_output_variables = .False.
    print_defaults = " "
    do i = 1, iarg
        call getarg(i, arg)
        ! Options
        if (arg(:1) == "-") then
          ! remove first - if long form
          if (arg(1:2) == "--") arg = arg(2:)
          select case (trim(arg(2:)))
            case ("h", "help")
              call print_splash(version)
              call print_help
            case ("d", "defaults")
              ! Interpret next argument as module
              call getarg(i + 1, print_defaults)
              if (trim(print_defaults) == "") print_defaults = "all"
            case ("o", "output-variables")
              print_output_variables = .True.
              ! only print errors to stderr and never log to stdout
              master_log%stderr_threshold = log_error_level
              master_log%stdout_threshold = huge(1)
            case ("v", "version")
              write(master_log%stdout, "(A)") version
              stop
            case default
              call print_help("Unknown option: "//trim(arg))
          end select
        else  ! Positional argument
          config_file_path = arg
        end if
    end do
  end subroutine parse_commandline_arguments

end module input
