module time

  use utilities, only : dp, log_info, log_debug, log_progress

  implicit none

  ! total number of days of simulation
  integer, save :: nDaysSim
  ! number of years of simulations
  integer, save :: nbyr = 10
  ! current years (YYYY) real year
  integer, save :: iyr = 2000
  !
  integer, save :: ieap
  ! current year starting with 1
  integer, save :: iy = 0
  ! current month
  integer, save :: mo = 0
  ! current day
  integer, save :: ida = 0
  ! number of days in the current year
  integer, save :: nd = 0
  ! number of days
  integer, save :: ndsum = 0
  ! ??
  integer, save :: nt = 0
  integer, save, dimension(13) :: nc = (/ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366/)
  !
  integer, save :: iday
  !
  integer, save :: inday
  ! counter of the current day of simulation (not of current year)
  integer, save :: daycounter = 0

  namelist / TIME_PARAMETERS / &
    iyr, &
    nbyr


contains

  subroutine time_initialise
    use utilities, only : is_leap_year
    use input, only : get_config_fid
    integer i

    read(get_config_fid(), TIME_PARAMETERS)

    ! *** count total number of days in simulation
    nDaysSim = 0
    do i = iyr, iyr + nbyr - 1
      if (is_leap_year(i) ) then
        nDaysSim = nDaysSim + 366
      else
        nDaysSim = nDaysSim + 365
      end if
    end do
  end subroutine time_initialise

  subroutine time_day_length(daylen, ida, mb, pit, ylc, yls)
    integer, intent(in) :: ida
    integer, intent(in) :: mb
    real(dp), intent(in) :: pit
    real(dp), dimension(:), intent(in) :: ylc
    real(dp), dimension(:), intent(in) :: yls
    real(dp), dimension(:), intent(out) :: daylen
    integer k
    real(dp) ch, h, sd, xi
    !**** CALC day length (from clgen)
    do k = 1, mb
      xi = ida
      sd = .4102 * sin((xi - 80.25) / pit)
      ch = - yls(k) * tan(sd) / ylc(k)
      if (ch .lt. 1.) then
        if (ch .le. - 1.) then
          h = 3.1416
        else
          h = acos(ch)
        end if
      else
        h = 0.
      end if
      daylen(k) = 7.72 * h
    end do
  end subroutine time_day_length

  subroutine time_process_years
    integer year, month

    call log_info("time_process_years", 'SWIM starts simulation')

    do year = 1, nbyr
      call time_initialise_year(year)
      call log_info('time_process_years', 'Processing year', int=iyr)
      do month = 1, 12
        call time_process_month(month)
      end do
      call time_finish_year()
    end do
  end subroutine time_process_years

  ! ----------------------------------------------------------------------
  subroutine time_initialise_year(year)
  ! ----------------------------------------------------------------------
    use input, only : read_real_column
    use utilities, only : is_leap_year
    use crop, only : &
      arylda, &
      avylda, &
      mcrdb, &
      iyrrot, &
      nrotyrs
    use subbasin, only : &
      bRunoffdat, &
      discharge_input_file_id, &
      obs_discharge, &
      nqobs
    integer year, jj, k, skip

    !     iy - current year
    !     mo - current month
    !     nt - =1 if 365 days, =0 if 366 days
    !     nd - number of days in the current year

    iy = year

    mo = 1
    nt = 1
    if (is_leap_year(iyr)) nt = 0
    ida = 1
    nd = 366 - nt
    ndsum = ndsum + nd

    !#### INITIALIZE CROPS: CALL INITCROP

    do jj = 1, mcrdb
      avylda(iy, jj) = 0.
      arylda(iy, jj) = 0.
    end do

    if (bRunoffdat) then
      ! read all columns by index (+1 for date column)
      ! skipping back up before reading except for 1st
      do k = 1, nqobs
        skip = -nd
        if (k == 1) skip = 0
        call read_real_column(discharge_input_file_id, array=obs_discharge(:nd, k), index=k+1, skip=skip)
      end do
    end if

    ! count year of rotation (1-3)
    if (iyrrot < nrotyrs) then
      iyrrot = iyrrot + 1
    else
      iyrrot = 1
    end if
    if (iy == 1) iyrrot = 0

  end subroutine time_initialise_year


  ! ----------------------------------------------------------------------
  subroutine time_finish_year()
  ! ----------------------------------------------------------------------

    use output, only : output_nashsutcliffe_efficiency, output_year
    use crop, only : istyr
    use subbasin, only : &
      runs, &
      bRunoffDat, &
      obs_discharge

    if (bRunoffDat) then
      call output_nashsutcliffe_efficiency(obs_discharge(:nd, 1), runs(:nd), 1, istyr, iy)
    end if

    call output_year

    call log_debug("time_finish_year", 'Completed year =', int=iyr)

    return
  end subroutine time_finish_year

  ! ----------------------------------------------------------------------
  subroutine time_process_month(month)
  ! ----------------------------------------------------------------------

    use output, only : output_month
    use vegetation, only : ialpha
    use nutrient, only : &
      degNgrw, &
      degNsub, &
      degNsur, &
      degPsur, &
      retNgrw, &
      retNsub, &
      retNsur, &
      retPsur
    use snow, only : bSnowModule
    use soil, only : psp, rtn
    use output, only : nsb, nvsub
    use crop, only : icc, mfe
    use evapotranspiration, only : pit, radiation_switch
    use catchment, only : bSubcatch

    integer, intent(in) :: month

    mo = month

    !ls** compute number of days of the month
    if (mo .eq. 2) then
      inday = nc(mo + 1) - nc(mo) - nt
    else
      inday = nc(mo + 1) - nc(mo)
    endif

    do iday = 1, inday
      call time_process_day(month, iday)
    end do

    call output_month(iyr, month)

    if (month .eq. 12) iyr = iyr + 1

    call log_debug("time_process_month", "Completed month =", int=month)

    return
  end subroutine time_process_month

  subroutine time_process_day(mo1, iday)
    !**** CALLED IN MAIN
    !**** THIS SUBROUTINE COMPUTES ONE DAY
    !****
    !**** CALLED: IN MAIN
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !     ida - current day
    !     ieap - index for GRASS output in subbasin (yield)
    !     ieapu - index for GRASS output in subbasin (annual sums)
    !     xxswind - soil water index for basin
    !     xwysb - water yield for basin
    !     xnflow() - N flows for a chosen hydrotop
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    use catchment, only : bSubcatch, da, da9, wy
    use crop, only : icc, mfe
    use erosion, only : &
      chc, &
      chxk, &
      conn, &
      cpp, &
      er, &
      xcklsp, &
      xnorg, &
      xnorgp, &
      xporg, &
      xpsedp
    use evapotranspiration, only : &
      canev, &
      ecal, &
      evapotranspiration_radiation, &
      humi, &
      omega, &
      pit, &
      ra, &
      radiation_switch, &
      snoev, &
      tmn, &
      ylc, &
      yls
    use groundwater, only : &
      additionalGwUptake, &
      disLastDay, &
      revapst, &
      xet
    use hydrotope, only : smm
#ifdef with_netcdf
    use input, only : input_nc_read_climate
#endif
    use management, only : &
      bWAM_Module, &
      management_external_supply, &
      management_is_transfer_subbasin, &
      management_total_demand, &
      wam_d, &
      wam_y
    use nutrient, only : &
      degNgrw, &
      degNsub, &
      degNsur, &
      degPsur, &
      retNgrw, &
      retNsub, &
      retNsur, &
      retPsur, &
      xnflow, &
      yon, &
      yph
    use output, only : &
      area_tot_glacier, &
      area_tot_snow, &
      depth_ave_glacier, &
      depth_ave_snow, &
      ieapu, &
      nsb, &
      nvsub, &
      output_day
    use reservoir, only : &
      bRsvModule, &
      bRsvHydrograph, &
      reservoir_is_operational, &
      reservoir_process, &
      reservoir_subbasin, &
      rsvSubbasin, &
      xwysb, &
      xxswind
    use river, only : &
      accf, &
      river_route, &
      river_route_add, &
      river_transfer, &
      varoute
    use snow, only : &
      bSnowModule, &
      sml, &
      snowVal, &
      tmx, &
      ieapg
    use soil, only : pr, psp, rtn, xqd
    use subbasin, only : &
      aff, &
      bRunoffdat, &
      dart, &
      flu, &
      icodes, &
      ihouts, &
      inum1s, &
      inum2s, &
      mb, &
      mhyd, &
      neap, &
      nqobs, &
      obs_discharge, &
      obssb, &
      precip, &
      qtl, &
      runs, &
      runsubbasin, &
      sbar, &
      sbp, &
      sda, &
      smq, &
      smsq, &
      sub, &
      subbasin_initialise_subbasin, &
      subbasin_read_climate, &
      subcatch_id, &
      subouthyd, &
      subp, &
      sumcn, &
      susb, &
      tx, &
      wysb, &
      xeo, &
      xpercn, &
      xqi, &
      xsep, &
      xssf, &
      xssfn
    ! avoid continuation limit warning
    use subbasin, only : &
      xswind, &
      xyno3, &
      xysp, &
      yd
    use vegetation, only : daylen

    integer, intent(in) :: mo1, iday

    integer ii, iik, k
    integer inum1, inum2
    integer idum
    integer icode, icodep, ihout
    logical :: bRoute = .true. !#### RESERVOIR MODULE ####


    !ls** compute day number in the year
    if (mo1 .gt. 2) then
      ida = iday + nc(mo1) - nt
    else
      ida = iday + nc(mo1)
    endif

    daycounter = daycounter + 1

    area_tot_snow = 0.
    depth_ave_snow = 0.
    !###########################
    !#### SNOW MODULE ####
    !###########################
    if (bSnowModule) then
      ieapg = 1
      area_tot_glacier = 0.
      depth_ave_glacier = 0.
    end if
    !###########################

    xxswind = 0.
    xwysb = 0.
    ieap = 1
    ieapu = 1
    do ii = 1, 20
      xnflow(ii) = 0.
    end do
    do k = 1, nsb
      sub(k) = 0.
    end do

!#### CALL (NC_) READCLI - to read climate data daily
#ifdef with_netcdf
    call input_nc_read_climate(flu, humi, mb, ra, subp, tmn, tmx, tx)
#else
    call subbasin_read_climate(humi, mb, ra, subp, tmn, tmx, tx)
#endif
    call time_day_length(daylen, ida, mb, pit, ylc, yls)
    !### VA
    !*********** For Radiation Data generated by Hargreaves Samani
    if (radiation_switch > 0) then
      call evapotranspiration_radiation(ida, mb, tmx)
    end if
    !### VA

    !ls ndmo(mo) = ndmo(mo) + 1
    !ls dtot = dtot + 1.
    snoev = 0.
    pr = 0.

    !########################################################### START ROUTING
    !#### CALL subbasin, route, transfer, add
    do idum = 1, mhyd
      icode = icodes(idum)
      ihout = ihouts(idum)
      inum1 = inum1s(idum)
      inum2 = inum2s(idum)
      if (icode .gt. 0) then
        select case (icode)
          case (1) ! SUBBASIN command

            !#################################
            !#### WATER MANAGEMENT MODULE ####
            !#################################
            if (bWAM_Module) then
              if (ihout == 1) then
                wam_y(daycounter) = iyr
                wam_d(daycounter) = ida
              end if

              if (management_is_transfer_subbasin(inum1) ) then

                ! Summarise total inflow only from external sources and transmission losses
                ! Current day inflows from external sources are added to subbasins %inflow
                ! NOTE: Adding to routing variables sda(2, j) takes place in subroutine subbasin
                call management_external_supply(inum1, daycounter, ida, iyr)

                ! calculate total water demand of water user(s)
                ! NOTE: nothing is removed here, just computation of total demand
                call management_total_demand(inum1, daycounter, ida, iyr)

                ! Depending on whether the subbasin is a headwater or not
                ! inputs and outputs are added/removed in subroutines:
                ! 'subbasin' or 'add'

              end if
            end if
            !###########################

            !###########################
            !#### RESERVOIR MODULE ####
            !###########################
            if (bRsvModule) then
              ! if actual subbasin is reservoir skip subbasin call
              ! TODO: include res_active function
              if (rsvSubbasin(inum1) == 0 ) then
                if (subcatch_id(inum1) .ne. 0) then
                  !write(*,*) 'Inside swim.f95::: olai(1, 1) = ', olai(1, 1)
                  call runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd)
                  !write(*,*) 'Inside swim.f95::: olai(1, 1) After = ', olai(1,1)
                end if
              else
                if (reservoir_is_operational(iyr, ida, inum1) ) then
                  call subbasin_initialise_subbasin(canev, sml, xcklsp, xet, xnorg, xnorgp, xporg, xpsedp, xqd, yon, yph)
                  call reservoir_subbasin(inum1, da, aff, dart, flu, neap, precip, qtl, sbp, sda, smq, smsq, snowval, sub, subp, sumcn, susb, tx, varoute, wysb, xeo, xpercn, xqi, xsep, xssf, xssfn, xswind, xyno3, xysp, yd, yon, yph)
                else if (subcatch_id(inum1) .ne. 0) then
                  call runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd)
                end if
              end if
            ! calculate subbasin but only if subbasin is listed in subcatch.def
            ! for convenient subsetting of model eg. by subcatchment
            else if (subcatch_id(inum1) .ne. 0) then
              call runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd)
            end if
            !###########################

          case (2) ! ROUTE command
            bRoute = .true.
            !###########################
            !#### RESERVOIR MODULE ####
            !###########################
            if (bRsvModule) then
              ! If inum1 is a reservoir subbasin then...
              if (bRsvHydrograph(ihout) ) then
                ! reservoir input is stored in varoute(2, inum2) and varoute(8, inum2)
                ! skip routing
                ! modify varoute(2, inum1), varoute(8, inum1) in reservoir functions
                ! DO NOT call route(icode, ihout, inum1, inum2)
                if (reservoir_is_operational(iyr, ida, inum1) ) then
                  call reservoir_process(ihout, inum1, inum2, ecal, humi, ida, iyr, mo, nc, omega, ra, subp, tx, varoute)
                  bRoute = .false.
                end if
              else
                bRoute = .true.
                !call route(icode, ihout, inum1, inum2)
              end if
            end if
            !###########################

            if (bRoute) then
              call river_route(ihout, inum1, inum2, chc, chxk, conn, cpp, da9, dart, er, flu, ida, iy, iyr, revapst, runs, sbar, sub, susb, xysp, yd, yon, yph)
            end if

          case (3) ! not implemented
            ! do nothing

          case (4) ! not implemented
            ihout = ihouts (idum - 1)
            icodep = icodes (idum - 1)
            call river_transfer() ! not implemented

          case (5) ! ADD command
            call river_route_add(bRunoffdat, ihout, inum1, inum2, disLastDay, additionalGwUptake, bWAM_Module, daycounter, ida, iyr, mb, nqobs, obssb, obs_discharge, runs, subouthyd, inum1s, bRsvModule, reservoir_is_operational(iyr, ida, inum2), rsvSubbasin)

          case default ! do nothing
        end select
      end if ! (icode .gt. 0)
    end do ! do idum = 1, mhyd
    !########################################################### END ROUTING

    !**** Calc wy - water yield
    wy = sub(8) + sub(9) - sub(10) + sub(15) + 1.e-20

    !**** Correction - sediment yield
    sub(21) = sub(21) / da9

    !**** Calculation of monthly sums
    do iik = 1, nvsub
      smm(iik) = smm(iik) + sub(iik)
    end do

    !**** Calc monthly water discharge
    if (bRunoffdat) accf(7) = accf(7) + obs_discharge(ida, 1)
    accf(8) = accf(8) + runs(ida)

    call output_day(iyr, mo1, iday)

    call log_progress('time_process_day', ida, nd)

    call log_debug("time_process_day", "Completed day =", int=iday)
    return
  end subroutine time_process_day

!   subroutine time_process_day(year, month, day)
!     use hydrotope, only: hydrotope_process, hydrotope_initialise_vars_timestep, n_hyrotopes
!     use subbasin, only: subbasin_process, subbasin_initialise_timestep, n_subbasins

!     integer, intent(in) :: year, month, day
!     integer hydrotope, subbasin

!     do hydrotope = 1, n_hyrotopes
!       call hydrotope_process(hydrotope)
!     end do

!     do subbasin = 1, n_subbasins
!       call subbasin_process(subbasin)
!     end do

!     call catchment_process(year, month, day)
!     call output_day(year, day, month)

!   end subroutine time_process_day

end module time