!! Perform discharge routing
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  2.3 - 2nd December 2024  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 07/Oct/2016 | Original code |
! | 1.1      | 08/Feb/2023 | DischargePointInit and DischargePointExport added |
! | 1.2      | 12/Feb/2023 | routing parameters are no more stored in stream network. No more need of HydroNetwork module |
! | 1.3      | 25/Jan/2024 | module name changed from SurfaceRouting to DischargeRouting |
! | 1.4      | 26/Jan/2024 | module adjusted to cope with new module Diversions |
! | 1.5      | 10/Apr/2024 | observed reservoir downstream discharge read from file |
! | 1.6      | 11/Apr/2024 | observed diverted discharge read from file |
! | 1.7      | 18/Apr/2024 | discharge routing parameters assigned for different subbasins |
! | 1.8      | 08/May/2024 | Qin written to output instead of Qout, to fix the bug of zero discharge in outlet cell |
! | 1.9      | 06/Jun/2024 | log message bug fixed when diversion is outside channel|
! | 2.0      | 14/Jun/2024 | back to Qout written to output file. Qout is defined in last downstream cell|
! | 2.1      | 03/Jul/2024 | Diverted discharge from diversions included in Qlateral of Muskingum |
! | 2.2      | 04/Sep/2024 | pools variable moved to Reservoirs module |
! | 2.3      | 02/Dec/2024 | excess of water retained in snow added to Qlateral of Muskingum |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Route discharge in river network.
!
! River network includes different elements:
!
!  -  hillslope reaches
!  -  channel reaches
!  -  reservoirs
!  -  diversions
!
! The user sets the elements to simulate and their properties 
! with a configuration file like in the following example
!
!```
!export-channel-grid = 0
!
!masks-number = 2 # number of masks to assign channel parameters (at least 1 for the base-mask must exist)
!
![reservoirs]
!  file = ./conf/reservoirs.ini
!  dt = 10 # if dt = 0 reservoirs are not solved
!  dt-out = 3600 #optional, default = 0
 ! 
![diversions]
!  file = ./conf/diversions.ini
!  dt = 0 # 0 = suppresses diversion simulation. otherwise dt is set equals to discharge routing dt
!  dt-out = 3600
!  
![discharge-in]
!  scalar = 0.0  
!  
![discharge-out]
!  scalar = 0.0 
!
![discharge-lat]
!  scalar = 0.0
!  
![base-mask]
!  channel-initiation-method = area #area, ask, curvature
!  channel-initiation-threshold = 4000000 # ask 600000 m^2 area 3500000
!  hillslope-width = 200 #  cross section width (m)
!  hillslope-alpha = 45. # slope of trapezoidal section side bank (degree)
!  hillslope-ks = 2 #Strickler roughness coefficient m^1/3 s^-1
!
![mask1]
!  file = ./data/subbasin1.asc
!  format = esri-ascii
!  epsg = 32633
!  channel-initiation-method = area #area, ask, curvature
!  channel-initiation-threshold = 400000 # ask 600000 m^2 area 3500000
!  hillslope-width = 200 #  cross section width (m)
!  hillslope-alpha = 45. # slope of trapezoidal section side bank (degree)
!  hillslope-ks = 3 #Strickler roughness coefficient m^1/3 s^-1
!
!Table Start
!Title: channel properties
!Id: base-mask
!#threshold indicates the basin area below which values of parameters are applied
!Columns:    [count]        [threshold]        [width]   [alpha]         [ks]             
!Units:        [-]            [m^2]             [m]      [deg]          [m^1/3s^-1]             
!			1            5000000                5         45            20           
!			2           10000000                7         45            25               
!			3           15000000               10         45            30               
!			4           20000000               20         45            35               
!			5           30000000               25         45            40           
!			6          100000000               35         45            45
!			7          500000000               50         45            45
!			8         1000000000               65         45            45
!			9         2000000000               80         45            45
!Table End
!
!Table Start
!Title: channel properties
!Id: mask1
!#threshold indicates the basin area below which values of parameters are applied
!Columns:    [count]        [threshold]        [width]   [alpha]         [ks]             
!Units:        [-]            [m^2]             [m]      [deg]          [m^1/3s^-1]             
!			1            5000000               5          45            15           
!			2           10000000               9          45            20               
!			3           15000000               20         45            25               
!			4           20000000               25         45            30               
!			5           30000000               30         45            35           
!			6          100000000               35         45            40
!			7          500000000               40         45            40
!			8         2000000000               45         45            45
!Table End         
!
!```
!
MODULE DischargeRouting      


! Modules used: 

USE DataTypeSizes, ONLY : &  
  ! Imported Parameters:
  float, &
  short

USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE GridLib, ONLY : &
  !imported definitions:
  grid_real, &
  grid_integer, &
  !imported routines:
  NewGrid, &
  GridDestroy, &
  ExportGrid, &
  ESRI_ASCII

USE GridOperations, ONLY : &
  !Imported routines:
  GridByIni, &
  CRSisEqual,  &
  CellArea,  & 
  !Imported operators:
  ASSIGNMENT (=)


USE IniLib, ONLY: &
  !Imported routines:
  IniOpen, &
  IniClose, &
  SectionIsPresent, &
  KeyIsPresent, &
  IniReadString, &
  IniReadInt,&
  IniReadReal, &
  !Imported type definitions:
  IniList

USE StringManipulation, ONLY : &
  !Imported routines:
  ToString, &
  StringCompact, &
  StringToUpper

USE Morphology, ONLY : &
  !Imported rituines:
  DownstreamCell, &
  CheckOutlet

USE Reservoirs, ONLY : &
  !imported definition:
  Reservoir, &
  !imported routines:
  InitReservoirs, &
  GetReservoirById, &
  CountReservoirs, &
  OutReservoirsInit, &
  OutReservoirs, &
  !imported variables:
  pools

USE Diversions, ONLY : &
    !imported definition:
    Diversion, &
    !Imported routines:
    InitDiversions, &
    OutDiversionsInit, &
    GetDiversionById, &
    OutDiversions, &
    !Imported variables:
    dtDiversion, &
    dtOutDiversion, &
    diversionChannels

USE Chronos, ONLY : &
  !imported definition:
  DateTime, &
  !Imported operators:
  OPERATOR (==), &
  OPERATOR (+), &
  ASSIGNMENT( = )

USE ObservationalNetworks, ONLY : &
  !Imported routines:
  ReadData, &
  ReadMetadata, &
  WriteData, &
  WriteMetadata, &
  CopyNetwork, &
  DestroyNetwork, &
  AssignDataFromGrid, &
  !Imported defined variable:
  ObservationalNetwork

USE RoutingModels, ONLY : &
  !Imported routines:
  LevelPool, &
  DivertFlow, &
  MuskingumCungeTodini, &
  MuskingumCunge

USE TableLib, ONLY: &
  !imported definitions:
  Table, &
  TableCollection, &
  !Imported routines:
  TableNew, &
  TableGetValue

USE Irrigation, ONLY: &
    !imported variables:
    Qirrigation

USE Snow, ONLY: &
    !imported variables:
    excessWaterSnowFlux

USE Utilities, ONLY : &
  !Imported routines:
  GetUnit

USE MorphologicalProperties, ONLY : &
  !imported variables:
  flowAccumulation, &
  flowDirection, &
  dem, &
  streamNetwork, &
  streamNetworkCreated

USE RiverDrainage, ONLY : &
  !imported routines:
  ChannelInitiation 



IMPLICIT NONE 


!Public routines
PUBLIC :: DischargeRoute
PUBLIC :: InitDischargeRouting
PUBLIC :: DischargePointInit

!public declarations:
TYPE(grid_real) :: Qin !! inlet discharge at time t (current)
TYPE(grid_real) :: Qout !! outlet discharge at time t (current)
TYPE(grid_real) :: Pin !! inlet discharge at time t-1 (previous)
TYPE(grid_real) :: Pout !! outlet discharge at time t-1 (previous)
TYPE(grid_real) :: Plat !! lateral discharge at time t-1 (previous)
TYPE(grid_real) :: waterDepth !! water depth from thalweg (m)
TYPE(grid_real) :: topWidth !! channel top width (m)

INTEGER (KIND = short) :: dtDischargeRouting


!Local routines
PRIVATE :: DischargePointExport

!Local declarations:
INTEGER (KIND = short), PRIVATE, &
    PARAMETER :: MISSING_DEF_INT  = -9999
REAL (KIND = float)   , PRIVATE, &
    PARAMETER :: MISSING_DEF_REAL = -9999.9

INTEGER (KIND = short), PRIVATE, &
    PARAMETER :: SLOPES           = 0
INTEGER (KIND = short), PRIVATE, &
    PARAMETER :: MCT              = 1 !Muskingum-Cunge-Todini
INTEGER (KIND = short), PRIVATE, &
    PARAMETER :: MUSKINGUM        = 2 !Muskingum
INTEGER (KIND = short), PRIVATE, &
    PARAMETER :: INSTANTANEOUS    = 3 !Infinite celerity (Lake)
INTEGER (KIND = short), PRIVATE, &
    PARAMETER :: MC               = 4 !Muskingum-Cunge

TYPE (ObservationalNetwork), PRIVATE :: sites
TYPE (ObservationalNetwork), PRIVATE :: sitesDischarge
INTEGER (KIND = short), PRIVATE :: fileUnitPointDischarge
TYPE (DateTime), PRIVATE :: timePointExport
TYPE (DateTime), PRIVATE :: timePoolsExport
TYPE (DateTime), PRIVATE :: timeDiversionsExport
INTEGER (KIND = short) :: dtOutReservoirs
TYPE (grid_integer), PRIVATE :: channel
TYPE (grid_integer), PRIVATE :: dams !!location of reservoirs in raster map
TYPE (grid_integer), PRIVATE :: divChannels !!location of diversions in raster map
TYPE (grid_real), PRIVATE :: manning !! Manning roughness (s / m^1/3)
TYPE (grid_real), PRIVATE :: bankSlope !!river bank slope  (degree)
TYPE (grid_real), PRIVATE :: sectionWidth !!river section width  (m)
TYPE (TableCollection), PRIVATE :: channelParameters

!=======
CONTAINS
!=======
! Define procedures contained in this module. 
  
!==============================================================================
!| Description:
!   Initialize surface routing
SUBROUTINE InitDischargeRouting  &
!
(fileini, time, path_out, flowdirection, domain, & 
  storage, netRainfall, dtrk )

IMPLICIT NONE

!Arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: fileini !!name of configuration file
TYPE(DateTime), INTENT(IN)      :: time
CHARACTER (LEN = *), INTENT(IN) :: path_out  !!path of output folder
TYPE(grid_integer), INTENT(IN)  :: flowdirection
TYPE(grid_integer), INTENT(IN)  :: domain

!Arguments with intent out:
TYPE(grid_real), INTENT(OUT)          :: storage  !water stored in channel cell
TYPE(grid_real), INTENT(OUT)          :: netRainfall ! net rainfall rate [m/s]
INTEGER(KIND = short), INTENT(OUT)    :: dtrk

!local declarations:
TYPE (IniList) :: ini
INTEGER (KIND = short) :: i, j, k 
!CHARACTER (LEN=5)     :: ic !!option for initial condition
TYPE(Reservoir), POINTER :: p !!point to current reservoir
TYPE(Diversion), POINTER :: d !!point to current diversion channel
CHARACTER (LEN = 30) :: channelInitiationMethod !!method to assign channel initiation
REAL (KIND = float) :: channelInitiationThreshold !!channel initiation threshold (m2)
REAL (KIND = float) :: hillslopelWidth !!hillslope section width (m)
REAL (KIND = float) :: hillslopeKs !!hillslope roughness (m^1/3 s-1)
REAL (KIND = float) :: hillslopeBankSlope !!hillslope bank slope (degree)
REAL (KIND = float) :: area !! basin area (m2)
REAL (KIND = float) :: channelValue
INTEGER (KIND = short) :: exportChannelGrid !!set channel grid export
REAL (KIND = float) :: scalar
INTEGER (KIND = short) :: nmasks !! number of masks to assign routing parameters
TYPE (grid_integer) :: maskTmp
CHARACTER (LEN = 100) :: stringTmp

!-------------------------------end of declaration-----------------------------

!check stream network

IF ( .NOT. streamNetworkCreated ) THEN
   CALL Catch ('error', 'DischargeRouting',   &
          'stream network is not available, check morphological properties' )
END IF

!load configuration file
CALL IniOpen (fileini, ini)

!number of masks to assign discharge routing parameters
IF ( KeyIsPresent ('masks-number', ini) ) THEN
    nmasks = IniReadInt ('masks-number', ini)
ELSE
    CALL Catch ('error', 'DischargeRouting',   &
			   'masks-number missing in configuration file' )
END IF


!channel initiation
!------------------

!allocate channel grid
CALL NewGrid ( channel, domain, 0)

!set option to export channel grid
exportChannelGrid = IniReadInt ('export-channel-grid', ini)

!first use base mask

channelInitiationMethod = IniReadString ('channel-initiation-method', &
                                         ini, section = 'base-mask')
channelInitiationThreshold = IniReadReal ('channel-initiation-threshold', &
                                         ini, section = 'base-mask')

CALL ChannelInitiation (method = channelInitiationMethod, &
                           threshold = channelInitiationThreshold, &
                           mask = domain, flowAcc = flowAccumulation,&
                           flowDir = flowDirection, dem = dem, &
                           channel = channel ) 


! modify channel initiation grid if additional masks are required
DO k = 1, nmasks - 1
    
    stringTmp = 'mask' // ToString (k)
    
    channelInitiationMethod = IniReadString ('channel-initiation-method', &
                                         ini, section = stringTmp)
    channelInitiationThreshold = IniReadReal ('channel-initiation-threshold', &
                                         ini, section = stringTmp)
    
    CALL GridByIni (ini, maskTmp, section = stringTmp)

    CALL ChannelInitiation (method = channelInitiationMethod, &
                           threshold = channelInitiationThreshold, &
                           mask = maskTmp, flowAcc = flowAccumulation,&
                           flowDir = flowDirection, dem = dem, &
                           channel = channel ) 
    CALL GridDestroy (maskTmp)
    
END DO


!export channel grid
IF ( exportChannelGrid == 1 ) THEN
    CALL ExportGrid (channel, 'channel.asc', ESRI_ASCII)
END IF

!create discharge parameter maps
!-------------------------------

!initialize empty maps
CALL NewGrid (layer = manning, grid = domain, initial = 0. )
CALL NewGrid (layer = bankSlope, grid = domain, initial = 0. )
CALL NewGrid (layer = sectionWidth, grid = domain, initial = 0. )

!read channel parameter table(s)
CALL TableNew ( fileini, channelParameters )

!assign parameters to base mask
hillslopelWidth = IniReadReal ('hillslope-width', ini, section = 'base-mask')
hillslopeKs = IniReadReal ('hillslope-ks', ini, section = 'base-mask')
hillslopeBankSlope = IniReadReal ('hillslope-alpha', ini, section = 'base-mask')
 

DO i = 1, domain % idim
    DO j = 1, domain % jdim
        IF ( domain % mat (i,j) /= domain % nodata ) THEN
            area = flowAccumulation % mat (i,j) * CellArea (domain,i,j)
            IF ( channel % mat (i,j) == 0 ) THEN !hillslope cell
                manning % mat (i,j) = 1. / hillslopeKs
                bankSlope % mat (i,j) = hillslopeBankSlope
                sectionWidth % mat (i,j) = hillslopelWidth
            ELSE !channel cell
                
                !set manning roughness coefficient
                CALL TableGetValue ( valueIn =  area, &
                                     tables = channelParameters, &
                                     id = 'base-mask',  keyIn = 'threshold', &
                                     keyOut ='ks', match = 'linear', &
                                     valueOut = channelValue, &
                                     bound = 'extendlinear' )
                
                manning % mat (i,j) = 1. / channelValue
                
                !set river bank slope
                CALL TableGetValue ( valueIn =  area, &
                                     tables = channelParameters, &
                                     id = 'base-mask',  keyIn = 'threshold', &
                                     keyOut ='alpha', match = 'linear', &
                                     valueOut = channelValue, &
                                     bound = 'extendlinear' )
                
                bankSlope % mat (i,j) = channelValue
                
                !set section width
                CALL TableGetValue ( valueIn =  area, &
                                     tables = channelParameters, &
                                     id = 'base-mask',  keyIn = 'threshold', &
                                     keyOut ='width', match = 'linear', &
                                     valueOut = channelValue, &
                                     bound = 'extendlinear' )
                
                sectionWidth % mat (i,j) = channelValue
                
            END IF
        END IF
    END DO
END DO


!assign parameters to sub masks

DO k = 1, nmasks - 1
    
    stringTmp = 'mask' // ToString (k)

    hillslopelWidth = IniReadReal ('hillslope-width', ini, &
                                   section = stringTmp)
    hillslopeKs = IniReadReal ('hillslope-ks', ini, &
                               section = stringTmp)
    hillslopeBankSlope = IniReadReal ('hillslope-alpha', ini, &
                                       section = stringTmp)
    
    CALL GridByIni (ini, maskTmp, section = stringTmp)
 
    DO i = 1, domain % idim
        DO j = 1, domain % jdim
            IF ( maskTmp % mat (i,j) /= maskTmp % nodata ) THEN
                area = flowAccumulation % mat (i,j) * CellArea (domain,i,j)
                IF ( channel % mat (i,j) == 0 ) THEN !hillslope cell
                    manning % mat (i,j) = 1. / hillslopeKs
                    bankSlope % mat (i,j) = hillslopeBankSlope
                    sectionWidth % mat (i,j) = hillslopelWidth
                ELSE !channel cell
                
                    !set manning roughness coefficient
                    CALL TableGetValue ( valueIn =  area, &
                                         tables = channelParameters, &
                                         id = stringTmp,  keyIn = 'threshold', &
                                         keyOut ='ks', match = 'linear', &
                                         valueOut = channelValue, &
                                         bound = 'extendlinear' )
                
                    manning % mat (i,j) = 1. / channelValue
                
                    !set river bank slope
                    CALL TableGetValue ( valueIn =  area, &
                                         tables = channelParameters, &
                                         id = stringTmp,  keyIn = 'threshold', &
                                         keyOut ='alpha', match = 'linear', &
                                         valueOut = channelValue, &
                                         bound = 'extendlinear' )
                
                    bankSlope % mat (i,j) = channelValue
                
                    !set section width
                    CALL TableGetValue ( valueIn =  area, &
                                         tables = channelParameters, &
                                         id = stringTmp,  keyIn = 'threshold', &
                                         keyOut ='width', match = 'linear', &
                                         valueOut = channelValue, &
                                         bound = 'extendlinear' )
                
                    sectionWidth % mat (i,j) = channelValue
                
                END IF
            END IF
        END DO
    END DO
    
    
    CALL GridDestroy (maskTmp)
    
END DO

!initial condition
!-----------------

!input discharge
IF (SectionIsPresent('discharge-in', ini )) THEN
     IF (KeyIsPresent ('scalar', ini, 'discharge-in') ) THEN
        scalar = IniReadReal ('scalar', ini, 'discharge-in')
        CALL NewGrid (Pin, domain, scalar)
    ELSE
        CALL GridByIni (ini, Pin, section = 'discharge-in')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( Pin, domain, 0. )
END IF

!output discharge
IF (SectionIsPresent('discharge-out', ini )) THEN
     IF (KeyIsPresent ('scalar', ini, 'discharge-out') ) THEN
        scalar = IniReadReal ('scalar', ini, 'discharge-out')
        CALL NewGrid (Pout, domain, scalar)
    ELSE
        CALL GridByIni (ini, Pout, section = 'discharge-out')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( Pout, domain, 0. )
END IF

!lateral discharge
IF (SectionIsPresent('discharge-lat', ini )) THEN
     IF (KeyIsPresent ('scalar', ini, 'discharge-lat') ) THEN
        scalar = IniReadReal ('scalar', ini, 'discharge-lat')
        CALL NewGrid (Plat, domain, scalar)
    ELSE
        CALL GridByIni (ini, Plat, section = 'discharge-lat')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( Plat, domain, 0. )
END IF

!allocate variables
CALL NewGrid (layer = Qin, grid = domain, initial = 0. )
CALL NewGrid (layer = Qout, grid = domain, initial = 0. )
CALL NewGrid (layer = storage, grid = domain, initial = 0. )
 
!net rainfall grid
CALL NewGrid (layer = netRainfall, grid = domain, initial = 0.)

!water depth grid
CALL NewGrid (layer = waterDepth, grid = domain, initial = 0.)

!topwidth 
CALL NewGrid (layer = topWidth, grid = domain, initial = 0.)

!diversions
IF (SectionIsPresent('diversions', ini)) THEN
    
    dtDiversion = IniReadInt ('dt', ini, section = 'diversions')
    IF ( dtDiversion /= 0 ) THEN
        !set dtDiversion = dtDischargeRouting
        dtDiversion = dtDischargeRouting
    END IF
    
    IF (dtDiversion > 0) THEN
        CALL InitDiversions ( IniReadString('file', ini, &
                              section = 'diversions') )
        
        !create grid with diversion location
        CALL NewGrid (divChannels, domain, domain % nodata ) 
        
        d => diversionChannels
        DO 
          i = d % r
          j = d % c
          divChannels % mat (i,j) = d % id
          
          IF ( channel % mat (i,j) == 0 ) THEN
             CALL Catch ('warning', 'DischargeRouting',   &
			   'diversion is not on channel: ', argument = ToString(d%id) )
          END IF
          
            
          IF ( .NOT. ASSOCIATED (d % next) ) THEN !found last element of list
              EXIT
          END IF
          
          !set next diversion
          d => d % next
        END DO
        
        !initialize files for output
        IF ( KeyIsPresent ('dt-out', ini, 'diversions' ) ) THEN
            
            dtOutDiversion = IniReadInt ('dt-out', ini, 'diversions')
            
            IF ( dtOutDiversion > 0) THEN
               timeDiversionsExport = time
               CALL OutDiversionsInit ( path_out )
            END IF
            
        END IF
    
    END IF
    
ELSE !section is mandatory: stop the program if missing
   CALL Catch ('error', 'DischargeRouting',   &
			   'error loading diversions: ' ,  &
			    argument = 'section not defined in ini file' )
END IF

!reservoirs
IF (SectionIsPresent('reservoirs', ini)) THEN
    
    !create grid with reservoirs location
    CALL NewGrid (dams, domain, domain % nodata ) 
    
    
    dtrk = IniReadInt ('dt', ini, section = 'reservoirs')
    IF (dtrk > 0) THEN
        CALL InitReservoirs (IniReadString('file', ini, &
                             section = 'reservoirs'), &
                             time, domain, pools)
        
        p => pools
        DO 
          i = p % r
          j = p % c
          dams % mat (i,j) = p % id
          
          IF ( channel % mat (i,j) == 0 ) THEN
             CALL Catch ('warning', 'DischargeRouting',   &
			   'reservoir is not on channel: ', argument = ToString(p%id) )
          END IF
          
            
          IF ( .NOT. ASSOCIATED (p % next) ) THEN !found last element of list
              EXIT
          END IF
  
           p => p % next
        END DO
    
            
        !initialize files for output
        IF ( KeyIsPresent ('dt-out', ini, 'reservoirs' ) ) THEN
            
            dtOutReservoirs = IniReadInt ('dt-out', ini, 'reservoirs')
            
            IF ( dtOutReservoirs > 0) THEN
               timePoolsExport = time
               CALL OutReservoirsInit (pools, path_out)
            END IF
            
        END IF
    
    END IF
    

ELSE !section is optional
    
   dtrk = 0
   
 END IF
 
!finished configuration, deallocate memory
CALL IniClose (ini)


END SUBROUTINE InitDischargeRouting
  

  
  
!==============================================================================
!| Description:
!   route  discharge in surface network
SUBROUTINE DischargeRoute  &
!
(dt, time, flowdir, runoff, dtrk,  riverToGroundwater, &
    groundwaterToRiver, storage )

IMPLICIT NONE

!Arguments with intent(in):
INTEGER(KIND = short), INTENT(IN)    :: dt !!time step [s]
TYPE(DateTime), INTENT(IN)           :: time
TYPE(grid_integer), INTENT(IN)       :: flowdir !!flow direction
TYPE(grid_real), INTENT(IN)          :: runoff !! net rainfall [m/s]
INTEGER(KIND = short), INTENT(IN)    :: dtrk !!time step for reservoirs 
TYPE(grid_real), INTENT(IN)          :: riverToGroundwater  !!discharge 
                                            !!from river to groundwater (m3/s)
TYPE(grid_real), INTENT(IN)          :: groundwaterToRiver !!discharge 
                                             !!from groundwater to river (m3/s)
!Arguments with intent inout
TYPE(grid_real), INTENT(INOUT) :: storage !!volume stored on channel cell 

!local declarations:
INTEGER(KIND = short)    :: k, iin, jin, is, js
REAL (KIND = float)      :: ddx
REAL (KIND = float)      :: Qlateral !!lateral discharge (m3/s)
REAL (KIND = float)      :: QlateralChannel !!lateral discharge in a diversion channel (m3/s)
REAL (KIND = float)      :: Qdownstream !! discharge downstream a diversion channel (m3/s)
REAL (KIND = float)      :: Qdiverted !! discharge diverted from a diversion channel (m3/s)
REAL (KIND = float)      :: runoff_ij !!runoff of single cell 
REAL (KIND = float)      :: wDepth
REAL (KIND = float)      :: tWidth
REAL (KIND = float)      :: area
TYPE (Reservoir), POINTER :: p !!pointer to current reservoir
TYPE (Diversion), POINTER :: d, dnest !!pointer to current diversion
!-------------------------------end of declaration----------------------------- 

!reset Qin 
Qin  = 0.

!loop troughout hydro network
DO k = 1, streamNetwork % nreach
  iin = streamNetwork % branch (k) % i0
  jin = streamNetwork % branch (k) % j0
  
  !follow the branch
  DO WHILE (  .NOT.( jin == streamNetwork % branch (k) % j1 .AND. &
                    iin == streamNetwork % branch (k) % i1  )   )
    !set runoff_ij
    runoff_ij = runoff % mat(iin,jin)
    
     
    !find downstream cell
    CALL DownstreamCell (iin, jin, flowdir % mat(iin,jin), &
                        is,js, ddx, flowdir)
    
    ! looking for reservoir
    IF ( dtrk > 0) THEN
        IF ( dams % mat (iin,jin) /= dams % nodata ) THEN
            
            !set current reservoir
            p => GetReservoirById (list = pools, id = dams % mat (iin,jin) )
            
            IF (time == p % tReadNewStage) THEN
                !read new stage value from file
                CALL ReadData (network = p % network, &
                            fileunit = p % unit, time = time)
                
                !update target level
                p % stageTarget = p % network % obs (1) % value
                
                 p % tReadNewStage = p % network % time
                
                !update time of next reading of target stage
                p % tReadNewStage = p % tReadNewStage + &
                                    p % network % timeIncrement
            
            END IF
            
            
            IF (p % dischargeDownstream .AND. &
                time == p % tReadNewDischargeDownstream ) THEN
                !read new downstream discharge value from file
                CALL ReadData (network = p % networkDischargeDownstream, &
                            fileunit = p % unitDischargeDownstream, time = time)
                
                !update time of next reading of downstream discharge
                p % tReadNewDischargeDownstream = &
                         p % tReadNewDischargeDownstream + &
                         p % networkDischargeDownstream % timeIncrement
            END IF
                
            IF (p % dischargeDiverted .AND. &
                time == p % tReadNewDischargeDiverted ) THEN
                !read new diverted discharge value from file
                CALL ReadData (network = p % networkDischargeDiverted, &
                            fileunit = p % unitDischargeDiverted, time = time)
                
                !update time of next reading of diverted discharge
                p % tReadNewDischargeDiverted = &
                         p % tReadNewDischargeDiverted + &
                         p % networkDischargeDiverted % timeIncrement
            END IF    
           
            !reservoir routing
            Qin % mat(iin,jin) = Qin % mat(iin,jin) + &
                                    runoff_ij * &
                                    CellArea(runoff,iin,jin) 
            CALL LevelPool (time, dt, dtrk, Pin % mat(iin,jin), &
                            Qin % mat(iin,jin), p)
        
            !check and simulate diversion from reservoir
            IF ( p % bypassIsPresent ) THEN
                !divert flow from river
                d => p % bypass
        
                CALL DivertFlow (time, d, p % Qout, p % Qout )
                
                !overwrite diversion inlet discharge when observation is available
                IF ( p % dischargeDiverted ) THEN
                   IF ( p % networkDischargeDiverted % obs (1) % value /= &
                        p % networkDischargeDiverted % nodata ) THEN
                      d % QinChannel = p % networkDischargeDiverted % obs (1) % value 
                   END IF
                END IF
        
               !route discharge into channel
               Qlateral = 0.
               CALL MuskingumCungeTodini ( dt, &
                                d % channelLenght, &
                                d % QinChannel, &
                                d % PinChannel, &
                                d % PoutChannel, &
                                Qlateral, Qlateral, &
                                d % channelWidth, &
                                d % channelBankSlope, &
                                d % channelSlope, &
                                d % channelManning, &
                                d % QoutChannel, &
                                wDepth, tWidth )
           
               !copy data of current step for next step
               d % PinChannel  = d % QinChannel
               d % PoutChannel = d % QoutChannel
    
            END IF
           
            !set Qout
            SELECT CASE (  p % typ )
            CASE ( 'on' )
                 IF ( p % dischargeDownstream ) THEN
                   IF ( p % networkDischargeDownstream % obs (1) % value /= &
                        p % networkDischargeDownstream % nodata ) THEN
                      !overwrite reservoir downstream discharge
                      p % Qout = p % networkDischargeDownstream % obs (1) % value 
                   END IF
                END IF
                Qout % mat(iin,jin) = p % Qout
            CASE ( 'off' )  !off-stream reservoir  DA RIVEDERE!!
                !compute off-stream pool outflow
                CALL TableGetValue ( valueIn =  p % stage, tab = p % geometry, keyIn = 'h', &
                                    keyOut ='Qout', match = 'linear', &
                                    valueOut = p % Qout_off, bound = 'extendlinear' )
                p % Pout_off = p % Qout_off
                Qout % mat(iin,jin) = p % Qout
                !assign outflow to pool out cell
                !IF (p % rout /= MISSING_DEF_INT .AND. p % cout /= MISSING_DEF_INT ) THEN
                !Qin % mat(p % rout,p % cout) = Qin % mat(p % rout,p % cout) + &
                !    p % Qout_off
                !END IF
             
            END SELECT
            Pin % mat(iin,jin) = Qin % mat(iin,jin)
            Pout % mat(iin,jin) = Qout % mat(iin,jin) 
            jin = js
            iin = is
            CYCLE
        END IF
    END IF
        
    !Muskingum-Cunge-Todini
    
    area = CellArea(runoff,iin,jin)
        
    !set Qlateral
    Qlateral = runoff_ij * area
        
    !remove irrigation
    IF ( ALLOCATED (Qirrigation % mat) ) THEN
        Qlateral = Qlateral - Qirrigation % mat (iin, jin)
    END IF
        
    !include river-groundwater exchange
    IF ( ALLOCATED ( riverToGroundwater % mat ) ) THEN
        IF ( riverToGroundwater % mat (iin,jin) /= &
                riverToGroundwater % nodata ) THEN
           
            Qlateral = Qlateral + groundwaterToRiver % mat (iin,jin) - &
                                    riverToGroundwater % mat (iin,jin) 
            
        END IF
    END IF
    
    !add excess of water retained in snow
    IF ( ALLOCATED ( excessWaterSnowFlux % mat ) ) THEN
        Qlateral = Qlateral + excessWaterSnowFlux % mat (iin, jin) * area
    END IF
        
        
    !add outflow from by-pass channel or off-stream basin
    p => pools
    DO WHILE (ASSOCIATED(p)) !loop trough all reservoirs
        IF ( p % typ == 'off' ) THEN
            IF ( iin == p % rout .AND. jin == p % cout ) THEN
                    Qlateral = Qlateral + p % Pout_off
            END IF
        END IF
            
        IF (p % bypassIsPresent) THEN
            IF ( iin == p % bypass % rout .AND. jin == p % bypass % cout ) THEN
                    Qlateral = Qlateral + p % bypass % PoutChannel
            END IF
        END IF
            
        p => p % next
    END DO
    
    !remove diverted discharge from diversion channel
    IF (dtDiversion > 0) THEN
        IF ( divChannels % mat (iin,jin) /= divChannels % nodata ) THEN
            !divert flow from river
            d => GetDiversionById (list = diversionChannels, &
                              id = divChannels % mat (iin,jin) )
        
            CALL DivertFlow (time, d, Qin % mat(iin,jin), Qdownstream )
            
            Qdiverted = Qin % mat (iin, jin) - Qdownstream
           
            !update Qlateral
            Qlateral = Qlateral -  Qdiverted
            
            !route discharge into channel
            QlateralChannel = 0.
            CALL MuskingumCungeTodini ( dtDiversion, &
                                d % channelLenght, &
                                d % QinChannel, &
                                d % PinChannel, &
                                d % PoutChannel, &
                                QlateralChannel, QlateralChannel, &
                                d % channelWidth, &
                                d % channelBankSlope, &
                                d % channelSlope, &
                                d % channelManning, &
                                d % QoutChannel, &
                                wDepth, tWidth )
    
           !copy data of current step for next step
            d % PinChannel  = d % QinChannel
            d % PoutChannel = d % QoutChannel
        END IF
    END IF
    
    
        
    !add outflow from diversion channel
    d => diversionChannels
        
    DO WHILE (ASSOCIATED(d)) !loop trough all diversion channels
            
        IF ( iin == d % rout .AND. jin == d % cout ) THEN
            ! As the flood routing is solved from upstream to downstream
            ! PoutChannel contains the current or the previous time step discharge according to the 
            ! Horton order of the outlet section respect to the horton order of onlet section 
            Qlateral = Qlateral + d % PoutChannel
        END IF
            
        d => d % next
    END DO
        
    CALL MuskingumCungeTodini ( dt, ddx, Qin % mat(iin,jin), &
                                Pin % mat(iin,jin), &
                                Pout % mat(iin,jin), &
                                Qlateral, Plat % mat(iin,jin), &
                                sectionWidth % mat (iin,jin), &
                                bankSlope % mat (iin,jin), &
                                streamNetwork % branch (k) % slope, &
                                manning % mat (iin,jin), &
                                Qout % mat(iin,jin), &
                                waterDepth % mat(iin,jin), &
                                topWidth % mat(iin,jin) )
    
    storage % mat (iin,jin) = storage % mat (iin,jin) + &
                ( Qin % mat(iin,jin) + Qlateral - Qout % mat(iin,jin) ) * &
                dt / CellArea(Qin,iin,jin)
           
    IF ( storage % mat (iin,jin) < 0. ) THEN
        !storage % mat (iin,jin) = 0.
        !Qin % mat(iin,jin) = 0.
        !Qout % mat(iin,jin) = 0.
    END IF
              
    IF (isnan (Qout % mat(iin,jin)) ) THEN
        Qout % mat(iin,jin) = Pout % mat(iin,jin)
    END IF
      
    !computed Qout becomes Qin for the downstream section. 
    !Take account of junctions, sum to existing discharge
    Qin % mat(is,js) = Qin % mat(is,js) + Qout % mat(iin,jin)
       
    !store previous values for next time step
    Pin % mat(iin,jin) = Qin % mat(iin,jin)
    Pout % mat(iin,jin) = Qout % mat(iin,jin)
    Plat % mat(iin,jin) = Qlateral
    
    !store cell indexes and select downstream cell for next loop
    jin = js
    iin = is
    
    !check next cell is outlet cell
    CALL DownstreamCell (iin, jin, flowdir % mat(iin,jin), &
                        is,js, ddx, flowdir)
    
    IF ( CheckOutlet (iin,jin,is,js,flowdir) ) THEN
        !next cell will not be computed, set approximate value of Qout
        Qout % mat (iin, jin) =  Qin % mat (iin, jin)  + &
            runoff % mat(iin,jin) *  CellArea(runoff,iin,jin) 
    END IF
    
    
    
    
  END DO
END DO

!write results
IF ( time == timePointExport ) THEN
   CALL DischargePointExport ( time )
   timePointExport = timePointExport + sitesDischarge % timeIncrement
END IF

IF ( time == timePoolsExport ) THEN
   CALL OutReservoirs ( pools, time, Qin, Qout )
   timePoolsExport = timePoolsExport + dtOutReservoirs
END IF

IF ( time == timeDiversionsExport ) THEN
   CALL OutDiversions ( diversionChannels, time, Qin, Qout )
   timeDiversionsExport = timeDiversionsExport + dtOutDiversion
END IF


RETURN
END SUBROUTINE DischargeRoute


!==============================================================================
!| Description:
!   Initialize export of point site data
SUBROUTINE DischargePointInit &
!
( pointfile, path_out, time )

IMPLICIT NONE

!Arguments with intent (in):
CHARACTER (LEN = *), INTENT(IN) :: pointfile  !!file containing coordinate of points
CHARACTER (LEN = *), INTENT(IN) :: path_out  !!path of output folder
TYPE (DateTime),     INTENT(IN) :: time  !!start time

!local declarations
INTEGER (KIND = short) :: err_io
INTEGER (KIND = short) :: fileunit

!-------------------------end of declarations----------------------------------

timePointExport = time

!open point file
fileunit = GetUnit ()
OPEN ( unit = fileunit, file = pointfile(1:LEN_TRIM(pointfile)), &
      status='OLD', iostat = err_io )

IF ( err_io /= 0 ) THEN
	!file does not exist
    CALL Catch ('error', 'DischargeRouting', 'out point file does not exist')
END IF 

!Read metadata
CALL ReadMetadata (fileunit, sites)

!check dt
IF (.NOT. ( MOD ( sites % timeIncrement, dtDischargeRouting ) == 0 ) ) THEN
  CALL Catch ('error', 'DischargeRouting', &
              'dt out sites must be multiple of dtDischargeRouting ')
END IF

CLOSE ( fileunit )

!create virtual network and initialize file for output

fileUnitPointDischarge = GetUnit ()
OPEN ( unit = fileUnitPointDischarge, &
    file = TRIM(path_out) // 'point_discharge.fts' )
    
CALL CopyNetwork ( sites, sitesDischarge )

sitesDischarge % description = 'discharge data exported from FEST'

sitesDischarge % unit = 'm3/s'

sitesDischarge % offsetZ = 0.

CALL WriteMetadata ( network = sitesDischarge, &
                fileunit = fileUnitPointDischarge )

CALL WriteData (sitesDischarge, fileUnitPointDischarge, .TRUE.)    

! destroy sites
CALL DestroyNetwork ( sites )

RETURN
END SUBROUTINE DischargePointInit


!==============================================================================
!| Description:
!   Export of point site data
SUBROUTINE DischargePointExport &
!
( time )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT (IN) :: time

!local declarations:
INTEGER (KIND = short) :: i

!-------------------------end of declarations----------------------------------


!set current time
sitesDischarge % time = time 
    
!populate data
!CALL AssignDataFromGrid (Qin, sitesDischarge )
!debug
CALL AssignDataFromGrid (Qout, sitesDischarge )
    
!write data
CALL WriteData (sitesDischarge, fileUnitPointDischarge)

RETURN
END SUBROUTINE DischargePointExport


END MODULE DischargeRouting