!! Define river network scheme 
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.2 - 22nd April 2024  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 10/Feb/2023 | Original code |
! | 1.1      | 18/Apr/2024 | ChannelInitiation modified for following DischargeRouting update |
! | 1.2      | 22/Apr/2024 | Flow direction convention imported from Morphology |
    
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! routines to define river network scheme for discharge routing
! The original algorithms were initially implemented in a 
! standalone program used to prepare input file to FEST.
! Starting from 2023 FEST release, river network
! topology orgamìnization is computed within the FEST
! before starting the simulation from the flow direction
! and flow accumulation maps.
!
MODULE RiverDrainage         


!Modules used:

USE GridLib, ONLY: &
    !Imported type definitions:
    grid_real, &
    grid_integer, &
    !Imported routines:
    NewGrid, &
    GridDestroy

USE GridOperations, ONLY: &
    !Imported routines:
    IsOutOfGrid, &
    GetIJ, &
    GetXY, & 
    CellArea

USE DataTypeSizes, ONLY: &
    !Imported definitions:
    short, &
    float, &
    long

USE LogLib, ONLY : &          
    ! Imported Routines:
    Catch

USE StringManipulation, ONLY: &
    !Imported routines:
    ToString

USE Morphology, ONLY: &
    !Imported routines:
    HortonOrders , &
    CheckOutlet, &
    DownstreamCell, &
    DeriveSlope, & 
    CellIsSpring

USE GeoLib, ONLY: &
    !Imported definitions:
    Coordinate, &
    !Imported routines:
    Distance , &
    !Imported assignment:
    ASSIGNMENT (=)

USE Utilities, ONLY: &
    !imported routines:
    GetUnit

USE Morphology, ONLY : &
    !Imported variables:
    NW, & !North-West
    N, &  !North
    NE, & !North-East
    W, &  !West
    E, &  !East
    SW, & !South-West
    S, &  !South
    SE    !South-East


IMPLICIT NONE 

! Local declarations: 

! Local routines:
PRIVATE :: SplitNetwork
PRIVATE :: ExportReaches
PRIVATE :: ExportShape

TYPE Reach
	INTEGER :: id
	!====================================================
	REAL(KIND = float) :: x0 !! beginning of reach x coordinate
	REAL(KIND = float) :: y0 !! beginning of reach y coordinate
	REAL(KIND = float) :: x1 !! end of reach x coordinate
	REAL(KIND = float) :: y1 !! end of reach y coordinate
	!====================================================
	INTEGER :: i0	!! beginning of reach row
	INTEGER :: j0	!! beginning of reach column
	INTEGER :: i1	!! end of reach row
	INTEGER :: j1	!! end of reach column
	!====================================================
	INTEGER :: ncells !!number of cells in a reach
	INTEGER :: order !! Horton-Sthraler order
	REAL(KIND = float) :: slope	!! average reach slope (m/m)
	REAL(KIND = float) :: length	!!reach length (m)
	REAL(KIND = float) :: area	!! area drained by end section (m2)
END TYPE Reach

TYPE ReachNetwork
	TYPE(Reach),POINTER :: branch (:)
	INTEGER :: nreach
END TYPE ReachNetwork

TYPE ReachesList
        INTEGER :: id
        INTEGER :: i0
        INTEGER :: j0
        INTEGER :: i1
        INTEGER :: j1
        REAL (KIND = FLOAT) :: y0
        REAL (KIND = FLOAT) :: x0
        REAL (KIND = FLOAT) :: y1
        REAL (KIND = FLOAT) :: x1
        INTEGER :: n_cells
        REAL (KIND = FLOAT) :: slope ! (L/L)
        REAL (KIND = FLOAT) :: length !(m)
        REAL (KIND = FLOAT) :: area ! (cells)
        INTEGER :: order	!Horton-Strhaler order	
        TYPE (ReachesList), POINTER :: next  !dynamic list
END TYPE ReachesList     

!Global routines:
PUBLIC :: BuildReachNetwork
PUBLIC :: ChannelInitiation

!private 
TYPE (ReachesList), POINTER, PRIVATE     :: list, current, previous


!=======
CONTAINS
!=======
! Define procedures contained in this module. 


!==============================================================================
! Description:
!   Define cells where channels begin. The routine searches for channel cells
!   that have no other input channel cells.
!SUBROUTINE MarkChannelBeginning 
!
!local declaration:
!INTEGER :: i,j
!INTEGER :: row, col
!LOGICAL :: inputFound
!
!-----------------------end of declaration-------------------------------------
!
!DO i = 1, channelBeginning % idim
!  DO j = 1, channelBeginning % jdim
!    IF (channel % mat (i,j) /= channel % nodata) THEN
!      inputFound = .FALSE.
!      
!      !check EAST cell
!      row = i 
!      col = j + 1
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == W .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check SOUTH-EAST cell
!      row = i + 1
!      col = j + 1
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == NW .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check SOUTH cell
!      row = i + 1
!      col = j
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == N .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check SOUTH-WEST cell
!      row = i + 1
!      col = j - 1
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == NE .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check WEST cell
!      row = i 
!      col = j - 1
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == E .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check NORTH-EAST cell
!      row = i - 1
!      col = j - 1
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == SE .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check NORTH cell
!      row = i - 1
!      col = j
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == S .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      !check NORTH-EAST cell
!      row = i - 1
!      col = j + 1
!      IF ( .NOT. IsOutOfGrid(row,col,channelBeginning) ) THEN
!        IF (flowDirection % mat (row,col) == SW .AND. &
!            channel % mat (row,col) /= channel % nodata ) THEN
!           inputFound = .TRUE.
!           channelBeginning % mat (i,j) = channelBeginning % nodata
!           CYCLE
!        END IF
!      END IF
!      
!      IF ( .NOT. inputFound ) THEN
!        channelBeginning % mat (i,j) = 1
!      END IF
!      
!    ELSE
!      channelBeginning % mat (i,j) = channelBeginning % nodata
!    END IF
!    
!    
!  END DO
!END DO
!
!
!END SUBROUTINE MarkChannelBeginning


!==============================================================================
!| Description:
!   Split channel network where a confluence of two different horton order 
!   channels occurs
SUBROUTINE SplitNetwork  &
!
(orders, flowDirection, split)

IMPLICIT NONE

!Arguments with intent in:
TYPE(grid_integer), INTENT(in) :: orders
TYPE(grid_integer), INTENT(in) :: flowDirection

!Arguments with intent out
TYPE(grid_integer),INTENT(inout):: split


!local declaration:
INTEGER :: i,j
INTEGER :: row, col
INTEGER :: hortonOrder
LOGICAL :: splitFound

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

DO i = 1, split % idim
  DO j = 1, split % jdim
    !IF (channel % mat (i,j) /= channel % nodata) THEN
      hortonOrder = orders % mat (i,j) 
      splitFound = .FALSE.
      
      !check EAST cell
      row = i 
      col = j + 1
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == W .AND. &
            !channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check SOUTH-EAST cell
      row = i + 1
      col = j + 1
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == NW .AND. &
           ! channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder ) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check SOUTH cell
      row = i + 1
      col = j
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == N .AND. &
          !  channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder ) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check SOUTH-WEST cell
      row = i + 1
      col = j - 1
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == NE .AND. &
          !  channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder ) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check WEST cell
      row = i 
      col = j - 1
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == E .AND. &
          !  channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check NORTH-EAST cell
      row = i - 1
      col = j - 1
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == SE .AND. &
         !   channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check NORTH cell
      row = i - 1
      col = j
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == S .AND. &
        !    channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder ) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      !check NORTH-EAST cell
      row = i - 1
      col = j + 1
      IF ( .NOT. IsOutOfGrid(row,col,split) ) THEN
        IF (flowDirection % mat (row,col) == SW .AND. &
        !    channel % mat (row,col) /= channel % nodata .AND. &
            orders % mat (row,col) <  hortonOrder) THEN
           splitFound = .TRUE.
           split % mat (i,j) = 1
           CYCLE
        END IF
      END IF
      
      IF ( .NOT. splitFound ) THEN
        split % mat (i,j) = split % nodata
      END IF
      
    !ELSE
    !  split % mat (i,j) = split % nodata
    !END IF
    
    
  END DO
END DO

END SUBROUTINE SplitNetwork 

!==============================================================================
!| Description:
!   Define channel cells. Two methods are possible:
!
!   1. area: channel begins where drainage basin exceedes a given area
!   2. ask: channel begins where AS<sup>k</sup> exceedes a given value, where A is 
!        the basin area (m<sup>2</sup>), S is the local slope (-), and k is the
!        basin fractal dimension (default = 1.7)
SUBROUTINE ChannelInitiation &
!
(method, threshold, mask, flowAcc, flowDir, dem, channel )


IMPLICIT NONE

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: method
REAL (KIND = float), INTENT(IN) :: threshold
TYPE (grid_integer), INTENT(IN) :: mask
TYPE (grid_integer), INTENT(IN) :: flowAcc
TYPE (grid_integer), INTENT(IN) :: flowDir
TYPE (grid_real), INTENT(IN) :: dem

!arguments with intent(inout):
TYPE (grid_integer), INTENT(INOUT) :: channel

!local declarations:
INTEGER :: i,j,iu,ju,id,jd
REAL    :: area !! ( m<sup>2</sup> )
REAL    :: scale = 1.7
REAL    :: ask !!local value of  AS<sup>k</sup>
TYPE (grid_real) :: slope

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

!initialize channel to zero
 DO i = 1, dem % idim
    DO j = 1, dem % jdim
        IF (mask % mat (i,j) /= mask % nodata ) THEN
            channel % mat (i,j) = 0
        END IF
    END DO
 END DO

SELECT CASE (method)
  CASE ('area')
    DO i = 1, mask % idim
      DO j = 1, mask % jdim
        IF (mask % mat (i,j) /= mask % nodata ) THEN
          area = flowAcc % mat (i,j) * CellArea (flowAcc, i, j)
          IF (area >= threshold) THEN
             channel % mat (i,j) = 1
          ELSE
            channel % mat (i,j) = 0
          END IF
        END IF
      END DO
    END DO
  
  CASE ('ask')
     !compute slope in radians
    CALL DeriveSlope (dem, slope)
    
    DO i = 1, mask % idim
      DO j = 1, mask % jdim
        IF (mask % mat (i,j) /= mask % nodata ) THEN
           area = flowAcc % mat (i,j) * CellArea (flowAcc, i, j)
           ask = area * TAN (slope % mat (i,j)) ** scale
           IF (ask >= threshold) THEN !channel initiation found
             !set channel on all downstream cell
             iu = i
             ju = j
             DO
               channel % mat (iu,ju) = 1
               CALL DownstreamCell (iu,ju,flowDir % mat(iu,ju),id,jd)
    
               IF ( CheckOutlet (iu,ju,id,jd,flowDir) ) EXIT
               iu = id
               ju = jd                            
             END DO 
             
           END IF
          
        END IF
      END DO
    END DO 
    CALL GridDestroy (slope) 
 !other methods to be implemented
END SELECT

RETURN
END SUBROUTINE ChannelInitiation



!==============================================================================
!| Description:
!   Build reaches that compose the river network
SUBROUTINE BuildReachNetwork &
!
(maxReachLength, slopeCorrection, flowDirection, flowAcc, dem, &
    fileExport, shpExport, streamNetwork )


IMPLICIT NONE

!arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: maxReachLength  !!max length of a reach (m)
REAL (KIND = float), INTENT (IN) :: slopeCorrection !! slope value to correct negative values
TYPE (grid_integer), INTENT (IN) :: flowDirection !! flow direction
TYPE (grid_integer), INTENT (IN) :: flowAcc !! flow accumulation
TYPE (grid_real),    INTENT (IN) :: dem !! digital elevation model
INTEGER (KIND = short), INTENT (IN) :: fileExport, shpExport

!arguments with intent (out):
TYPE (ReachNetwork), INTENT (OUT) :: streamNetwork

!local declarations
LOGICAL                         :: outletSection
LOGICAL                         :: endOfReach
INTEGER                         :: countReaches
INTEGER                         :: maxOrder
INTEGER                         :: i,j,l !loop index
INTEGER                         :: row,col !current cell
INTEGER                         :: iDown, jDown !downstream cell
INTEGER                         :: n_cells
REAL                            :: reachLength

TYPE (grid_integer)             :: orderBeginning
TYPE (grid_integer)             :: orders
TYPE (grid_integer)             :: split
TYPE (Coordinate)               :: point1, point2

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

!------------------------------------------------------------------------------
!(1.0) Initialize orders, orderBeginning and split grids
!------------------------------------------------------------------------------

CALL NewGrid (orders, flowDirection)

CALL NewGrid (orderBeginning, flowDirection)

CALL NewGrid (split, flowDirection)

!------------------------------------------------------------------------------
!(2.0) Build Horton orders grid
!------------------------------------------------------------------------------

CALL HortonOrders (flowDirection, orders, maxOrder)

!------------------------------------------------------------------------------
!(3.0) Build orderBeginning and split network
!------------------------------------------------------------------------------

CALL MarkOrderBeginning (orders, flowDirection, orderBeginning)

!split network
CALL SplitNetwork (orders, flowDirection, split)

!------------------------------------------------------------------------------
!(4.0) Build reaches
!------------------------------------------------------------------------------

countReaches = 0
NULLIFY(list)  ! at the beginning list is empty

!initialize points to compute reach length
point1 % system = flowDirection % grid_mapping
point2 % system = flowDirection % grid_mapping

DO l =1, maxOrder

    CALL Catch ('info', 'RiverDrainage', 'Elaborating reaches of stream order: ', &
                 argument = ToString(l))


    DO j=1,orderBeginning%jdim         ! scan orderBeginning grid
      DO i=1,orderBeginning%idim

        IF(orderBeginning%mat(i,j) == l) THEN
           countReaches    = countReaches + 1 
           reachLength     = 0.
           n_cells         = 0
           row             = i
           col             = j
           endOfReach      = .FALSE.
           outletSection   = .FALSE.

           IF(.NOT.ASSOCIATED(list)) THEN
               ALLOCATE(list)     !set first reach
               current => list
           ELSE
               ALLOCATE(current%next) !allocate next reach
               NULLIFY(current % next % next)
               current => current%next
           ENDIF

           current % i0    = i
           current % j0    = j
           current % order = l
           current % id    = countReaches

    	  
           !follow the reach until an order change or to the outlet, 
           !if it is the maximum order    
           DO WHILE (.NOT. endOfReach) 
           
          
             CALL DownstreamCell (row,col,flowDirection%mat(row,col), &
                                 iDown,jDown)
             n_cells      = n_cells + 1
             
             CALL GetXY (row, col, flowDirection, point1 % easting, point1 % northing)
             CALL GetXY (iDown, jDown, flowDirection, point2 % easting, point2 % northing)
             reachLength  = reachLength + Distance (point1, point2) !length (m)
             
		     IF ( CheckOutlet (row,col,iDown,jDown,flowDirection) ) THEN
			    outletSection = .TRUE.
			    endOfReach    = .TRUE.
		     ENDIF

             IF ( outletSection ) THEN
               current % i1 =  row
               current % j1 =  col
		       current % n_cells = n_cells ! unit: cells (integer)
		       current % length = reachLength ! !length (m)
		       current % area = flowAcc % mat(row,col) ! * &
		                        !flowAccumulation % cellsize**2
             
		       		  		       		       		       		      	        		          
		       ! if the reach length exceedes the maximum length =>
		       ! break the reach
             ELSEIF (maxReachLength > 0. .AND. &  !if maximum length = 0 I do not break the reach
		             reachLength  >= maxReachLength ) THEN
		       current % i1 = row
		       current % j1 = col
		       current % n_cells = n_cells
		       current % length = reachLength
		       !current % routingMethod = default
		       current % area = flowAcc % mat(row,col) ! * &
		                        !flowAccumulation % cellsize**2 !(m2)
		       !set elements for the following reach
		       ALLOCATE(current % next)
		       NULLIFY(current % next % next)
		       current => current % next
		       current % i0    = row
		       current % j0    = col
		       current % order = l
		       countReaches  = countReaches + 1
		       reachLength   = 0. 
		       n_cells       = 0
		       current % id  = countReaches
		       
		     END IF	
		     
		     IF ( .NOT. IsOutOfGrid(iDown,jDown,orders)) THEN
		         IF (orders%mat(iDown,jDown) > l) THEN
                   endOfReach = .TRUE.
                   current % i1 = iDown
                   current % j1 = jDown
                   IF (n_cells == 0) THEN !this case occurs when a reach has just been terminated because too long
                       current % n_cells = 1
                       CALL GetXY (row, col, flowDirection, point1 % easting, point1 % northing)
                       CALL GetXY (iDown, jDown, flowDirection, point2 % easting, point2 % northing)
		               current % length = Distance (point1, point2)
                   ELSE
		               current % n_cells = n_cells
		               current % length = reachLength 
		           END IF
		           current % area = flowAcc % mat(row,col) ! * &
		                            !flowAccumulation % cellsize**2
		         END IF  
		     END IF
		      		 			                  
             row = iDown
             col = jDown
                         
           END DO                     
        
        ENDIF
                
      ENDDO
    ENDDO

ENDDO


!set other parameters
!scan all list elements
current => list
DO 

  !calculate slope from digital elevation model and check negative slope
  current % slope = ( dem % mat(current%i0,current%j0) - &
                      dem % mat(current%i1,current%j1) ) / &
                      current % length
                      
  IF ( current % slope <= 0. ) THEN
    IF ( slopeCorrection > 0. ) current % slope  = slopeCorrection
  ENDIF
    
   !calculate spatial coordinate
   CALL GetXY (current % i0, current % j0, dem, current % x0, current % y0)
   CALL GetXY (current % i1, current % j1, dem, current % x1, current % y1)
  
  IF ( .NOT. ASSOCIATED (current % next) )THEN !found last element of list
    EXIT
  END IF
  
  current => current % next
  
END DO


!------------------------------------------------------------------------------
!(5.0) remove temporary variables
!------------------------------------------------------------------------------

CALL GridDestroy ( orders )

CALL GridDestroy ( orderBeginning )


!------------------------------------------------------------------------------
!(6.0) export files
!------------------------------------------------------------------------------

IF ( fileExport > 0 ) THEN
    CALL ExportReaches ( )
END IF

IF ( shpExport > 0 ) THEN
    CALL ExportShape ( dem, flowDirection )
END IF

!------------------------------------------------------------------------------
!(7.0) populate stream network 
!------------------------------------------------------------------------------
ALLOCATE ( streamNetwork % branch (countReaches) ) 

streamNetwork % nreach = countReaches

current => list
DO i = 1, countReaches
    
  streamNetwork % branch (i) % id = current % id
  streamNetwork % branch (i) % i0 = current % i0
  streamNetwork % branch (i) % j0 = current % j0
  streamNetwork % branch (i) % i1 = current % i1
  streamNetwork % branch (i) % j1 = current % j1
  
  streamNetwork % branch (i) % x0 = current % x0
  streamNetwork % branch (i) % y0 = current % y0
  streamNetwork % branch (i) % x1 = current % x1
  streamNetwork % branch (i) % y1 = current % y1
  
  streamNetwork % branch (i) % ncells = current % n_cells
  streamNetwork % branch (i) % slope = current % slope
  streamNetwork % branch (i) % length = current % length
  streamNetwork % branch (i) % area = current % area
  streamNetwork % branch (i) % order = current % order
    
  current => current % next

END DO
  
RETURN
END SUBROUTINE BuildReachNetwork

!==============================================================================
!| Description:
!   Export reaches on file
SUBROUTINE ExportReaches &
!
( )

IMPLICIT NONE

!local declaration:
INTEGER ( KIND = short ) :: fileunit

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

fileunit = GetUnit ()

OPEN (unit = fileunit, file = 'stream_network.txt' )

WRITE (fileunit, '(a)') 'surface_routing_reaches'
WRITE (fileunit, '(a)') '          ID          X0            Y0             &
                            X1             Y1                CELLS     STHRALER_ORDER   &
                            SLOPE_L/L        LENGTH_m       DRAINED_CELLS   '

!scan all list elements
current => list
DO 

    WRITE(fileunit, '(I,3X,E,E,E,E,I,I,10x,E,1X,E,1X,E,I,11X,E,1X,E,&
                        1X,E,1X,E,1X,E,5X,E,E,E)') &
                        current % id, current % x0, current % y0, &
                        current % x1, current % y1, current % n_cells, &
                        current % order, current % slope, current % length,  &
                        current % area 
                            
                            
    IF ( .NOT. ASSOCIATED (current % next) ) THEN !found last element of list
    EXIT
    END IF

    current => current % next

END DO

CLOSE (fileunit)



END SUBROUTINE ExportReaches


!==============================================================================
!| Description:
!   Define cells where the beginning of a reach of different order takes place.
SUBROUTINE MarkOrderBeginning &
!
(orders, flowDirection, orderBeginning)

IMPLICIT NONE

!Arguments with intent in:
TYPE(grid_integer),INTENT(in):: orders
TYPE(grid_integer),INTENT(in):: flowDirection

!Arguments with intent out
TYPE(grid_integer),INTENT(inout):: orderBeginning

LOGICAL                         :: outlet

INTEGER                         :: row,col !current cell
INTEGER                         :: iDown, jDown !downstream cell
INTEGER                         :: reachOrder
INTEGER                         :: maxPathOrder
INTEGER                         :: i,j
!- End of header --------------------------------------------------------------


 DO j=1,orderBeginning%jdim        
  DO i=1,orderBeginning%idim
    IF(CellIsSpring(i,j,flowDirection)) THEN !found a spring
	   orderBeginning%mat(i,j) = 1
       row         = i
       col         = j
       outlet      = .FALSE.
       reachOrder  = 1

       DO WHILE (.NOT. outlet) ! follow the reach till the outlet   

         CALL DownstreamCell(row,col,flowDirection%mat(row,col),iDown,jDown)          

         outlet = CheckOutlet (row,col,iDown,jDown,flowDirection)
         
         IF (.NOT. outlet) THEN                     
            IF(orders%mat(iDown,jDown) == reachOrder + 1) THEN
			   !found the beginning of a reach of increased order 
			   reachOrder = reachOrder + 1
               orderBeginning%mat(iDown,jDown) = reachOrder                                                                                                                                                                           
            ENDIF
         ENDIF
         row = iDown
         col = jDown
       END DO 
	                  
    ENDIF
  ENDDO
ENDDO  

!-----------remove duplicates------------

DO j=1,orderBeginning%jdim        ! scan entire grid
  DO i=1,orderBeginning%idim

    IF(CellIsSpring(i,j,flowDirection)) THEN !found a spring
       row          = i
       col          = j
       outlet       = .FALSE.
       maxPathOrder = 1

       DO WHILE (.NOT. outlet) ! follow the reach till the outlet    

         CALL DownstreamCell(row,col,flowDirection%mat(row,col),iDown,jDown)          

         outlet = CheckOutlet (row,col,iDown,jDown,flowDirection)
         
         IF (.NOT. outlet) THEN                     
            IF(orderBeginning%mat(iDown,jDown) <= maxPathOrder) THEN 
            !only one point where a N order reach begins is possible

               orderBeginning%mat(iDown,jDown) = 0       
                                                         
            ELSEIF (orderBeginning%mat(iDown,jDown) > maxPathOrder) THEN
               maxPathOrder = orderBeginning%mat(iDown,jDown)                                                                                                                                                                                  
            ENDIF
         ENDIF

         row = iDown
         col = jDown

       END DO                
    ENDIF

  ENDDO
ENDDO  

RETURN
END SUBROUTINE MarkOrderBeginning

!==============================================================================
!| Description:
!   export shape file of river network
SUBROUTINE ExportShape &
!
( dem, flowDirection )



IMPLICIT NONE

!arguments with intent (in):
TYPE (grid_real), INTENT (IN) :: dem
TYPE (grid_integer), INTENT (IN) :: flowDirection


!local declarations
INTEGER                           :: K, jin,iin,ivalle,jvalle
INTEGER                           :: dot
REAL                              :: X,Y
LOGICAL(4)                        :: res
CHARACTER(300)                    :: command
CHARACTER(300)                    :: filename
INTEGER                           :: system
!------------------------------------------------------------------------------

 
!------------------------------------------------------------------------------
!(1.1) create "generate file"
!------------------------------------------------------------------------------
OPEN(unit=876,file='tratti.gen')
OPEN(unit=877,file='tratti.csv')
WRITE(877,'(a)')'#ID;X0;Y0;X1;Y1;CELLS;STHRALER;SLOPE_L/L;LENGTH_m;&
                    DRAINED_CELLS'

!scan all list elements
current => list
DO 
       
    WRITE(876,*) current % id
    jin = current % j0
    iin = current % i0

    WRITE(877,'(I8,";",E14.7,";",E14.7,";",E14.7,";",E14.7,";",I8,";",I8,";",E14.7,";",&
                E14.7,";",E14.7,";",I8,";",E14.7,";",E14.7,";",E14.7,";",E14.7,";",E14.7)') &
	        current % id, current % x0, current % y0, &
            current % x1, current % y1, current % n_cells, &
            current % order, current % slope, current % length,  &
            current % area 
    	
    DO WHILE (.not.((jin.EQ.current%j1)  .AND. &
	                (iin.EQ.current%i1)))

	    X = DEM%xllcorner + jin *  DEM%cellsize - DEM%cellsize / 2.
	    Y = DEM%yllcorner + (DEM%idim - (iin-1)) * DEM%cellsize - DEM%cellsize / 2.

        WRITE(876,'(f20.7,a1,f20.7)') X,',',Y

        CALL DownstreamCell(iin,jin,flowDirection%mat(iin,jin),ivalle,jvalle)
    	 
	    iin = ivalle
	    jin = jvalle  
    	   													 
    END DO  !end single reach
       
        X = DEM%xllcorner + jin *  DEM%cellsize - DEM%cellsize / 2.
	    Y = DEM%yllcorner + (DEM%idim - (iin-1)) * DEM%cellsize - DEM%cellsize / 2.

        WRITE(876,'(f20.7,a1,f20.7)') X,',',Y	

        WRITE(876,'(a3)') 'end'
       
    IF ( .NOT. ASSOCIATED (current % next) )THEN !found last element of list
        EXIT
    END IF
       
    current => current % next
       
END DO	!end cycle on all reaches

WRITE(876,'(a3)') 'end'

CLOSE(876)
CLOSE(877)

!!------------------------------------------------------------------------------
!!(1.2) create shape file
!!------------------------------------------------------------------------------
filename = 'stream_network'
!dot = SCAN (filename,'.')
!IF (dot /= 0) filename = filename(1:dot-1)
command = 'gen2shp ' // filename(1:len_trim(filename)) // ' lines < tratti.gen'
res = SYSTEM(command) 

!create dbf
command='txt2dbf -C10 -R20.7 -R20.7 -R20.7 -R20.7 -I10 -I10 -R15.7 -R15.2  &
        -R15.2 -I10 -d ; tratti.csv ' // filename(1:len_trim(filename)) // '.dbf'
res = SYSTEM(command)

!------------------------------------------------------------------------------
!(1.3) delete intermediate files
!------------------------------------------------------------------------------

res = SYSTEM('del tratti.gen')
res = SYSTEM('del tratti.csv')




RETURN
END SUBROUTINE ExportShape



         

END MODULE RiverDrainage