module catchment

  use utilities, only : dp, path_max_length, log_debug, log_info, log_warn
  use input, only : get_config_fid, input_count_rows

  implicit none

  integer, save :: n_subcatch = 1
  ! Catchment ids read from catchment.csv
  integer, dimension(:), allocatable :: catchment_id
  ! area of subcatchment in [m2]
  real(dp), save, dimension(:), allocatable :: subcatch_area
  ! array storing annual subcatch output (nbyr, nSubcatch, 30)
  real(dp), save, dimension(:, :, :), allocatable :: subcatch_an
  ! basin area, km2
  real(dp), save :: da
  ! 1000. * da, basin area (1000 * km ** 2)
  real(dp), save :: af
  ! 100. * da = basin area in ha, from readbas
  real(dp), save :: da9
  !
  real(dp), save :: wy
  ! Subcatchment parameters on/off
  logical, save :: bSubcatch

  integer :: catchment_input_file_id

  character(len=path_max_length) :: catchment_input_file = "catchment.csv"

  namelist / CATCHMENT_PARAMETERS / &
    catchment_input_file, &
    bSubcatch

contains

  subroutine catchment_initialise(sbar, flu, subcatch_id, subcatch_idx)
    use input, only : input_open_file, get_config_fid, read_integer_column
    use utilities, only : int_index
    real(dp), intent(in) :: sbar(:)
    real(dp), intent(out) :: flu(:)
    integer, intent(in) :: subcatch_id(:)
    integer, intent(out) :: subcatch_idx(:)
    integer i, ii, minc, maxc

    read(get_config_fid(), CATCHMENT_PARAMETERS)

    if (bSubcatch) then
      catchment_input_file_id = input_open_file(catchment_input_file)
      n_subcatch = input_count_rows(catchment_input_file_id)
      allocate(catchment_id(n_subcatch))
      call read_integer_column(catchment_input_file_id, "catchment_id", catchment_id)
      call log_debug("catchment_initialise", "Catchment ids read from catchment file", \
        ints=catchment_id)
    else
      ! get catchment ids from subbasin catchment_id column, unique values
      minc = minval(subcatch_id)
      maxc = maxval(subcatch_id)
      n_subcatch = 0
      do i = minc, maxc
        if (any(subcatch_id == i)) n_subcatch = n_subcatch + 1
      end do
      allocate(catchment_id(n_subcatch))
      ii = 0
      do i = 1, n_subcatch
        if (any(subcatch_id == i)) then
          ii = ii + 1
          catchment_id(ii) = i
        end if
      end do
      call log_debug("catchment_initialise", \
        "Catchment ids inferred from subbasin catchment_id", ints=catchment_id)
    end if
    call log_info("catchment_initialise", "Number of catchments:", int=n_subcatch)

    ! assign indeces for subbasin-catchment mapping
    ! only translater ids if not 0
    subcatch_idx = 0
    do i = 1, size(subcatch_id)
      if (subcatch_id(i) > 0) &
        subcatch_idx(i) = int_index(subcatch_id(i), catchment_id)
    end do
    call catchment_allocate

    ! calculate total catchment (drainage) area "da" in km2
    da = sum(sbar(:)) / 10 ** 6
    da9 = 100. * da
    ! fraction of subbasin area of total catchment area
    flu(:) = sbar(:) / (da * 10 ** 6)
  end subroutine catchment_initialise

  subroutine catchment_allocate

  end subroutine catchment_allocate

  subroutine dealloc_catchment

  end subroutine dealloc_catchment

  subroutine catchment_initialise_parameters(mb, nbyr, sbar, subcatch_idx)
    !-------------------------------------------------------------------------------
    ! Author : stefan.liersch@pik-potsdam.de
    ! Date : 2010-02-24
    ! modified: 2010-02-25
    !
    ! PURPOSE : Reading file subcatch.def
    !           count number of subcatchments (user defined aggregation of subbasins)
    !
    ! CALLED : from subroutine main program
    !
    ! ToDo : - Writing output like pcp etc. at subcatchment level to specific output files
    !
    !-------------------------------------------------------------------------------
    integer, intent(in) :: mb
    integer, intent(in) :: nbyr
    real(dp), dimension(:), intent(in) :: sbar
    integer, dimension(:), intent(in) :: subcatch_idx
    integer :: i, n, nosub

    ! skip this routine if no subcatch parameters
    if (.not. bSubcatch) then
      return
    end if

    nosub = 0
    do n = 1, mb
      if (subcatch_idx(n) == 0) nosub = nosub + 1
    end do
    if (nosub > 0) then
      call log_warn("catchment_initialise_parameters", &
        "You are simulating with a subset of n subbasins:", int=mb - nosub)
    end if

    ! allocate arrays
    call catchment_initialise_subcatchm(n_subcatch, nbyr)
    do i = 1, mb
      if (subcatch_idx(i) /= 0) then
        subcatch_area(subcatch_idx(i)) = subcatch_area(subcatch_idx(i)) + sbar(i) !subarea(i)
      end if
    end do

    ! read individual subcatchment bsn parameters
    call catchment_read_subcatch_params

    ! assign parameters at the subcatchment level
    call catchment_assign_subcatch(mb, subcatch_idx)

  end subroutine catchment_initialise_parameters

  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  subroutine catchment_read_subcatch_params
    use input, only : read_integer_column, read_real_column
    use groundwater, only : &
      gw_abf, &
      gw_bff, &
      gw_delay, &
      gw_rchrgc, &
      gw_revapc, &
      gw_revapmn
    use soil, only: bsn_cncor, bsn_sccor
    use river, only: bsn_roc2, bsn_roc4
    use snow, only: &
      bsn_gmrate, &
      bsn_smrate, &
      bsn_tmelt, &
      bsn_tsnfall
    use evapotranspiration, only: bsn_thc, bsn_ecal
    use utilities, only: check_range

    ! read file subcatch.prm
    ! enable assignment of individual subbasin parameters

    call read_real_column(catchment_input_file_id, "ecal", bsn_ecal, 0.0_dp)
    call read_real_column(catchment_input_file_id, "thc", bsn_thc, 0.0_dp)
    call read_real_column(catchment_input_file_id, "roc2", bsn_roc2, 0.0_dp)
    call read_real_column(catchment_input_file_id, "roc4", bsn_roc4, 0.0_dp)
    call read_real_column(catchment_input_file_id, "cncor", bsn_cncor, 0.0_dp, range=(/0.25, 1.25/))
    call read_real_column(catchment_input_file_id, "sccor", bsn_sccor, 0.0_dp)
    call read_real_column(catchment_input_file_id, "tsnfall", bsn_tsnfall, 0.0_dp)
    call read_real_column(catchment_input_file_id, "tmelt", bsn_tmelt, 0.0_dp)
    call read_real_column(catchment_input_file_id, "smrate", bsn_smrate, 0.0_dp)
    call read_real_column(catchment_input_file_id, "gmrate", bsn_gmrate, 0.0_dp)
    call read_real_column(catchment_input_file_id, "bff", gw_bff, 0.0_dp)
    call read_real_column(catchment_input_file_id, "abf", gw_abf, 0.0_dp)
    call read_real_column(catchment_input_file_id, "delay", gw_delay, 0.0_dp)
    call read_real_column(catchment_input_file_id, "revapc", gw_revapc, 0.0_dp)
    call read_real_column(catchment_input_file_id, "rchrgc", gw_rchrgc, 0.0_dp)
    call read_real_column(catchment_input_file_id, "revapmn", gw_revapmn, 0.0_dp)

    call log_debug("catchment_read_subcatch_params", "Subcatchment parameters:")
    call log_debug("catchment_read_subcatch_params", "id", ints=catchment_id)
    call log_debug("catchment_read_subcatch_params", "ecal", reals=bsn_ecal)
    call log_debug("catchment_read_subcatch_params", "thc", reals=bsn_thc)
    call log_debug("catchment_read_subcatch_params", "roc2", reals=bsn_roc2)
    call log_debug("catchment_read_subcatch_params", "roc4", reals=bsn_roc4)
    call log_debug("catchment_read_subcatch_params", "cncor", reals=bsn_cncor)
    call log_debug("catchment_read_subcatch_params", "sccor", reals=bsn_sccor)
    call log_debug("catchment_read_subcatch_params", "tsnfall", reals=bsn_tsnfall)
    call log_debug("catchment_read_subcatch_params", "tmelt", reals=bsn_tmelt)
    call log_debug("catchment_read_subcatch_params", "smrate", reals=bsn_smrate)
    call log_debug("catchment_read_subcatch_params", "gmrate", reals=bsn_gmrate)
    call log_debug("catchment_read_subcatch_params", "bff", reals=gw_bff)
    call log_debug("catchment_read_subcatch_params", "abf", reals=gw_abf)
    call log_debug("catchment_read_subcatch_params", "delay", reals=gw_delay)
    call log_debug("catchment_read_subcatch_params", "revapc", reals=gw_revapc)
    call log_debug("catchment_read_subcatch_params", "rchrgc", reals=gw_rchrgc)
    call log_debug("catchment_read_subcatch_params", "revapmn", reals=gw_revapmn)

  end subroutine catchment_read_subcatch_params

  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  subroutine catchment_assign_subcatch(mb, subcatch_idx)
    !-------------------------------------------------------------------------------
    ! Author : stefan.liersch@pik-potsdam.de
    ! Date : 2010-02-25
    ! modified: 2010-02-26
    !
    ! PURPOSE : Assigning individual subbasin parameters read from subcatch.prm
    !
    ! CALLED : from subroutine read_subcatch_def
    !-------------------------------------------------------------------------------
    use groundwater, only : &
      gw_abf, &
      gw_bff, &
      gw_delay, &
      gw_rchrgc, &
      gw_revapc, &
      gw_revapmn, &
      abf, &
      bff, &
      delay, &
      gwht, &
      gwq, &
      syld, &
      rchrgc, &
      revapc, &
      revapmn
    use soil, only: bsn_cncor, bsn_sccor, cncor, sccor
    use river, only: bsn_roc2, bsn_roc4, roc2, roc4
    use snow, only: &
      bsn_gmrate, &
      bsn_smrate, &
      bsn_tmelt, &
      bsn_tsnfall, &
      smrate, &
      gmrate, &
      tmelt, &
      tsnfall
    use evapotranspiration, only: bsn_thc, bsn_ecal, thc, ecal

    integer, intent(in) :: mb
    integer, dimension(:), intent(in) :: subcatch_idx
    integer :: i, si

    do i = 1, mb
      si = subcatch_idx(i)
      if (si /= 0) then
        ! groundwater parameters
        bff(i) = gw_bff(subcatch_idx(i))
        gwht(i) = .5 ! gw_gwht(subcatch_idx(i))
        gwq(i) = .5 ! gw_gwq(subcatch_idx(i))
        abf(i) = gw_abf(subcatch_idx(i))
        syld(i) = .003 ! gw_syld(subcatch_idx(i))
        delay(i) = gw_delay(subcatch_idx(i))
        revapc(i) = gw_revapc(subcatch_idx(i))
        rchrgc(i) = gw_rchrgc(subcatch_idx(i))
        revapmn(i) = gw_revapmn(subcatch_idx(i))

        ! bsn parameters
        ecal(i) = bsn_ecal(subcatch_idx(i))
        thc(i) = bsn_thc(subcatch_idx(i))
        sccor(i) = bsn_sccor(subcatch_idx(i))
        roc2(i) = bsn_roc2(subcatch_idx(i))
        roc4(i) = bsn_roc4(subcatch_idx(i))
        cncor(i) = bsn_cncor(subcatch_idx(i))
        tsnfall(i) = bsn_tsnfall(subcatch_idx(i))
        tmelt(i) = bsn_tmelt(subcatch_idx(i))
        smrate(i) = bsn_smrate(subcatch_idx(i))
        gmrate(i) = bsn_gmrate(subcatch_idx(i))
      end if
    end do

  end subroutine catchment_assign_subcatch

  subroutine catchment_initialise_subcatchm(n, nbyr)
    use groundwater, only: groundwater_allocate_subcatch
    use evapotranspiration, only: evapotranspiration_allocate_sc
    use river, only: river_allocate_subcatch
    use soil, only: soil_allocate_subcatch
    use snow, only: snow_allocate_subcatch

    integer, intent(in) :: nbyr
    integer, intent(in) :: n

    allocate(subcatch_an(nbyr, n + 1, 30))
    subcatch_an = 0.
    allocate(subcatch_area(n))
    subcatch_area = 0.

    call groundwater_allocate_subcatch(n)
    call evapotranspiration_allocate_sc(n)
    call snow_allocate_subcatch(n)
    call river_allocate_subcatch(n)
    call soil_allocate_subcatch(n)

  end subroutine catchment_initialise_subcatchm

end module catchment