!! Compute sediment erosion and deposition 
!|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.0 - 5th October 2011  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 05/Oct/2011 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
!  
!   compute sediment erosion and deposition
!   The interril and channelerosion are computed according
!   to Sun et al. (2002), which combines the use of USLE
!   data source with storm-based, distributed hydrological model
!
! References:
!
!  Sun, H., P. S. Cornish, and T. M. Daniell, Contour-based digital elevation 
!      modeling of watershed erosion and sedimentation:
!      Erosion and sedimentation estimation tool (EROSET), 
!      Water Resour. Res., 38(11), 1233, doi:10.1029/2001WR000960, 2002.
!
MODULE Sediment         

! Modules used: 

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

USE GridLib, ONLY: &
! Imported type definitions:
grid_real, grid_integer, &
!Imported routines:
NewGrid, &
!debug
exportgrid, esri_ascii

USE LogLib, ONLY: &
!Imported routines:
Catch

USE DomainProperties, ONLY: &
  !imported variables:
  domain  => mask

USE MorphologicalProperties, ONLY: &
    !imported variables:
    dtm => dem

!USE CommonDataIO, ONLY: &
!!Imported variables:
!domain  => bacino_bilancio, &
!dtm

USE Morphology, ONLY: &
!Imported routines:
DeriveSlope, DownstreamCell

USE Units, ONLY: &
!Imported parameters:
gravityAccel

USE HydroNetwork, ONLY : &
  !Imported definitions:
  NetworkSediment, &
  !Imported routines:
  ReadHydroNetwork

!USE propagazione_superficiale, ONLY: &
!!Imported variables:
!gridC1, gridC2, gridC3

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

!USE modelli_propagazione, ONLY: &
!!Imported routines:
!muskingum_cunge

USE RoutingModels, ONLY: &
  !Imported routines:
  MuskingumCungeTodini

IMPLICIT NONE 
         
! Global declarations:
TYPE (grid_real), PUBLIC :: interrillErosion !!interrill sediment detachment rate (kg/m2/s)
TYPE (grid_real), PUBLIC :: channelErosion !!channel sediment detachment rate (kg/m2/s)
TYPE (grid_real), PUBLIC :: QoutSed !!output sediment discharge
TYPE (grid_real), PUBLIC :: deltaSed !!variation of sediment storage (10^6 kg)
LOGICAL,          PUBLIC :: routeSediment !!flag to route sediment
TYPE (grid_integer), PUBLIC :: sedFlowDirection !!map of flow direction to route sediment
TYPE (NetworkSediment), PUBLIC :: sedReach !!drainage network to route sediment



!Global Procedures:
PUBLIC :: SedimentInit
PUBLIC :: InterrillDetachmentRate
PUBLIC :: SedimentRouting

!Local Procedures:
PRIVATE :: ComputeSlopeFactor
PRIVATE :: ChannelDetachmentRate 
PRIVATE :: Yang1973
PRIVATE :: Yang1979
PRIVATE :: FallVelocity
PRIVATE :: ShearVelocity
PRIVATE :: CriticalVelocity

!Local variables:
TYPE (grid_real), PRIVATE :: slopeFactor !!slope dactor computed by routine ComputeSlopeFactor
TYPE (grid_real), PRIVATE :: slope !!local slope (radians)
TYPE (grid_real), PRIVATE :: rusleK !!soil erodibility factor ( ( tons h) / ( MJ mm))
TYPE (grid_real), PRIVATE :: rusleC !!crop and management factor (-)
TYPE (grid_real), PRIVATE :: QinSS !!input discharge of suspended sediment at time t
TYPE (grid_real), PRIVATE :: QoutSS !!output discharge of suspended sediment at time t
TYPE (grid_real), PRIVATE :: PinSS !!input discharge of suspended sediment at time t-1
TYPE (grid_real), PRIVATE :: PoutSS !!output discharge of suspended sediment at time t-1
TYPE (grid_real), PRIVATE :: QinBL !!input discharge of bed load sediment at time t
TYPE (grid_real), PRIVATE :: QoutBL !!output discharge of bed load at time t
TYPE (grid_real), PRIVATE :: PinBL !!input discharge of bed load at time t-1
TYPE (grid_real), PRIVATE :: PoutBL !!output discharge of bed load at time t-1


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

!==============================================================================
!| Description:
!  Initialize sediment related variables
SUBROUTINE SedimentInit &
!
(iniFile, dtRoute, fileOutSedimentRouting)

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

USE GridOperations, ONLY: &
!Imported routines
GridByIni, CRSisEqual, &
GridByIni

USE StringManipulation, ONLY: &
!Imported routines:
ToString

IMPLICIT NONE

!Argument with intent in:
CHARACTER (LEN = 300), INTENT(IN) :: iniFile !!file containing configuration information 

!Argument with intent out:
INTEGER (KIND = short), INTENT(OUT) :: dtRoute !!time step for sediment routing
CHARACTER (LEN = 300), INTENT(OUT)  :: fileOutSedimentRouting


!local declarations:
TYPE(IniList) :: iniDB !!store configuration info
CHARACTER (LEN = 300) :: filename
INTEGER (KIND = short) :: k, iin, jin, is, js

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

CALL Catch ('info', 'Sediment', 'initialize sediment module ')

!--------------------------------------------
!  open and read configuration file
!--------------------------------------------

CALL IniOpen (iniFile, iniDB)

!-------------------------------------------
!  load parameters and options
!-------------------------------------------

!soil erodibility factor
IF (SectionIsPresent('soil-erodibility', iniDB)) THEN
    CALL GridByIni (iniDB, rusleK, section = 'soil-erodibility')
    IF  ( .NOT. CRSisEqual (mask = domain, grid = rusleK, checkCells = .TRUE.) ) THEN
       CALL Catch ('error', 'Sediment',   &
		        'wrong spatial reference in soil erodibility factor' )
    END IF
ELSE
     CALL Catch ('error', 'Sediment: ',   &
		        'missing soil-erodibility section in configuration file' )
END IF

!crop and management factor
IF (SectionIsPresent('crop-factor', iniDB)) THEN
    CALL GridByIni (iniDB, rusleC, section = 'crop-factor')
    IF  ( .NOT. CRSisEqual (mask = domain, grid = rusleC, checkCells = .TRUE.) ) THEN
       CALL Catch ('error', 'Sediment',   &
		        'wrong spatial reference in crop and management factor' )
    END IF
ELSE
     CALL Catch ('error', 'Sediment: ',   &
		        'missing crop-factor section in configuration file' )
END IF

!sediment routing
IF (SectionIsPresent('route-sediment', iniDB)) THEN
    dtRoute = IniReadInt ('time-step', iniDB, section = 'route-sediment')
    IF (dtRoute > 0) THEN
        
        routeSediment = .TRUE.
        !read file containing cross sections to be included in output file 
        fileOutSedimentRouting =  IniReadString ('xs-output', iniDB, section = 'route-sediment')
        
        !read flow direction
        CALL GridByIni (iniDB, sedFlowDirection, section = 'route-sediment', &
                        subsection = 'flow-direction' ) 
        !read drainage network
        filename = IniReadString ('sediment-reach', iniDB, section = 'route-sediment')
        CALL ReadHydroNetwork (filename=filename, &
                        domain = sedFlowDirection, network = sedReach)
        !check consistency of drainage network
        DO k = 1, sedReach % nreach
		        iin = sedReach % branch(k) % i0
		        jin = sedReach % branch(k) % j0
	            DO WHILE ( .NOT.((jin == sedReach % branch(k) % j1) .AND. &
			                      (iin == sedReach % branch(k) %i1)) )
		            IF(domain%mat(iin,jin) == domain%nodata) THEN
			          CALL Catch ('error', 'Sediment',   &
			                        'error in checking river drainage: ' ,  &
			                        argument = 'reach out of the basin row ' // &
			                        ToString(iin) // ' col ' // ToString(jin))
		            END IF
			        CALL DownstreamCell (iin, jin, sedFlowDirection % mat(iin,jin), is, js)
				          jin = JS
				          iin = IS
		        END DO
	            !last cell of last reach
	            IF (K == sedReach%nreach) THEN
		            IF(domain%mat(iin,Jin) == domain%nodata) THEN
		                CALL Catch ('error', 'Sediment',   &
	                          'error in checking drainage network: ' ,  &
	                        argument = 'reach out of the basin row ' // &
	                        ToString(iin) // ' col ' // ToString(jin)) 
		            END IF
	            END IF
        END DO
	    !initialize sediment routing grids
	    CALL NewGrid (QinSS, domain)
	    CALL NewGrid (QoutSS, domain)
	    CALL NewGrid (PinSS, domain)
	    CALL NewGrid (PoutSS, domain)
	    CALL NewGrid (QinBL, domain)
	    CALL NewGrid (QoutBL, domain)
	    CALL NewGrid (PinBL, domain)
	    CALL NewGrid (PoutBL, domain)
	    CALL NewGrid (QoutSed, domain) 
    ELSE
        routeSediment = .FALSE.
    END IF
    
ELSE
     routeSediment = .FALSE.
END IF

!-------------------------------------------------------
!             compute slope factor
!-------------------------------------------------------
!derive slope

CALL DeriveSlope (dtm, slope)

CALL ComputeSlopeFactor (slope)

!-------------------------------------------------------
!            initialize detachment rate grids
!-------------------------------------------------------

CALL NewGrid (interrillErosion, domain)

!-------------------------------------------------------
!    initialize !variation of sediment storage  grid
!-------------------------------------------------------
!assume initial value = 0.
CALL NewGrid (deltaSed, domain, 0.) 

!----------------------------------------------------
!  Configuration terminated. Deallocate ini database
!----------------------------------------------------

CALL IniClose (iniDB)

RETURN

END SUBROUTINE SedimentInit

!==============================================================================
!| Description:
!  compute sediment detachment rate (kg/s). Interril area detachment
!  is considered to be due to the raindrop impact, which breaks the
!  cohesive bonds between particles, making them available for transport.
!
! References:
!
!  Foster, G. R., Modeling the erosion process, in Hydrologic Modeling of
!    Small Watersheds, edited by C. T. Haan, H. P. Johnson, and D. L.
!    Brakensiek, pp. 297–383, Am. Soc. of Agric. Eng., St. Joseph, Minn., 1982   
SUBROUTINE InterrillDetachmentRate &
!
(pRate)

IMPLICIT NONE

!Argument with intent(in):
TYPE (grid_real), INTENT(IN) :: pRate !! precipitation rate (m/s)

!Local declarations:
INTEGER :: i,j

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

DO j = 1, interrillErosion % jdim
  DO i = 1, interrillErosion % idim
    IF (interrillErosion % mat (i,j) /= interrillErosion % nodata) THEN             
                                      !conversion m/s -> mm/h 
      interrillErosion % mat (i,j) = 0.0138 * (pRate % mat(i,j) * 3600000. )**2.0 * &
                            rusleK % mat(i,j) * slopeFactor % mat(i,j) * &
                            rusleC % mat(i,j) / 3600. * CellArea (interrillErosion, i,j)
                            
    END IF
  END DO
END DO

END SUBROUTINE InterrillDetachmentRate


!==============================================================================
!| Description:
!  compute sediment detachment rate per second in channel element (kg/s)
!  using simplified approach by Sun et al. (2002).
!
! References:
!
!  Sun, H., P. S. Cornish, and T. M. Daniell, Contour-based digital 
!  elevation modeling of watershed erosion and sedimentation:
!  Erosion and sedimentation estimation tool (EROSET), 
!  Water Resour. Res., 38(11), 1233, doi:10.1029/2001WR000960, 2002.  
FUNCTION ChannelDetachmentRate &
!
(flow, s, k, c) &
!
RESULT (d)

!Arguments with intent in:
REAL (KIND = float), INTENT(IN) :: flow !!water discharge in channel (m3/s)
REAL (KIND = float), INTENT(IN) :: s !!channel slope (m/m)
REAL (KIND = float), INTENT(IN) :: k !!soil erodibility factor ((tons h) / (MJ mm))
REAL (KIND = float), INTENT(IN) :: c !!crop and management factor (-)

!local declarations:
REAL (KIND = float) :: d
REAL (KIND = float), PARAMETER :: alpha = 4300. 
REAL (KIND = float), PARAMETER :: beta  = 2.

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

d = alpha * flow ** beta * s * k * c

RETURN
END FUNCTION ChannelDetachmentRate


!============================================================================
!| Description:
!  compute sediment transport capacity using Yang (1973) approach (kg/s)
!  limitation: Yang1973 equation can be applied for noncohesive natural 
!  beds with particle sizes between 0.062 mm and 2 mm, a specific 
!  gravity of 2.65 g/cm3, and a shape factor of 0.7
!
! References:
!
!  Yang, C. T., Unit stream power and sediment transport, 
!    J. Hydraul. Div. Am. Soc. Civ. Eng., 98(HY10), 1805– 1826, 1972.
!
!  Yang, C. T., Incipient motion and sediment transport, 
!    J. Hydraul. Div. Am. Soc. Civ. Eng., 99(HY10), 1679– 1705, 1973.
FUNCTION Yang1973 &
!
(flow, s, v, d, dp) &
!
RESULT (tc)

IMPLICIT NONE

!Arguments with intent in:
REAL (KIND = float), INTENT(IN) :: flow !!water discharge in channel (m3/s)
REAL (KIND = float), INTENT(IN) :: s !!channel slope (m/m)
REAL (KIND = float), INTENT(IN) :: v !!water flow velocity in channel (m/s)
REAL (KIND = float), INTENT(IN) :: d !!particle size of bed material (mm)
REAL (KIND = float), INTENT(IN) :: dp !!water depth (m)

!local declarations:
REAL (KIND = float) :: tc !!computed transport capacity (kg/s)
REAL (KIND = float) :: vs !! channel Unit Stream Power (m/s)
REAL (KIND = float) :: vsCrit !! critical channel Unit Stream Power (m/s)
REAL (KIND = float) :: kVisc = 0.000001004 !!kinematic viscosity of water at 20 °C(m2/s)
REAL (KIND = float) :: I, J !! terms to compute transport capacity in Yang method
REAL (KIND = float) :: conc !!sediment concentration (mg/l)

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

!compute channel unit stream power
vs = v * s

!compute channel critical unit stream power
vsCrit = CriticalVelocity(d, dp, s) * s

!compute I and J

!I = 5.435 - 0.286 * LOG (FallVelocity(d) * (d/1000.) / kVisc) - &
I = 6.0 - 0.286 * LOG (FallVelocity(d) * (d/1000.) / kVisc) - &
    0.457 * LOG (ShearVelocity(dp,s)/FallVelocity(d)) 
    
!J = 1.799 - 0.409 * LOG (FallVelocity(d) * (d/1000.) / kVisc) - &
J = 1.5 - 0.409 * LOG (FallVelocity(d) * (d/1000.) / kVisc) - &
    0.314 * LOG (ShearVelocity(dp,s)/FallVelocity(d)) 
    
    
!if (isnan (I)) then
!  write(*,*) 'I isnan', FallVelocity(d), ShearVelocity(dp,s)
!   pause
!end if    
!
!if (isnan (J)) then
!  write(*,*) 'J isnan', FallVelocity(d), ShearVelocity(dp,s)
!   pause
!end if    
    
!compute sediment concentration (ppm or mg/l)
IF (vs > vsCrit) THEN
   conc = EXP ( I + J * LOG((vs - vsCrit)/FallVelocity(d)) )
ELSE
   conc = 0.
END IF

IF (conc < 0.) THEN
  conc = 0.
END IF

!compute transport capacity (kg/s)
tc = flow * conc / 1000.


END FUNCTION Yang1973 




!============================================================================
!| Description:
!  compute sediment transport capacity using Yang (1979) approach (kg/s)
!  limitation: Yang1979 equation can be applied when the dimensionless
!  unit stream power is relatively small with respect to the
!  prevailing value of unit stream power. Used for finer suspended
!  sediment trasnport in the Yellow River with median particle
!  diameter of 0.067 mm.
!
! References:
!
!  Yang, C. T., Unit stream power equations for total load.
!    J. Hydro., Amsterdam, The Netherlands, Vol. 40, 123–138.
!
!  Yang, C. T., Molinas, A., Wu, B., Sediment transport in the Yellow River, 
!    J. Hydraul. Eng., 122(5), 237–244, 1996.
FUNCTION Yang1979 &
!
(flow, s, v, d, dp) &
!
RESULT (tc)

IMPLICIT NONE

!Arguments with intent in:
REAL (KIND = float), INTENT(IN) :: flow !!water discharge in channel (m3/s)
REAL (KIND = float), INTENT(IN) :: s !!channel slope (m/m)
REAL (KIND = float), INTENT(IN) :: v !!water flow velocity in channel (m/s)
REAL (KIND = float), INTENT(IN) :: d !!particle size of bed material (mm)
REAL (KIND = float), INTENT(IN) :: dp !!water depth (m)

!local declarations:
REAL (KIND = float) :: tc !!computed transport capacity (kg/s)
REAL (KIND = float) :: vs !! channel Unit Stream Power (m/s)
REAL (KIND = float) :: vsCrit !! critical channel Unit Stream Power (m/s)
REAL (KIND = float) :: kVisc = 0.000001004 !!kinematic viscosity of water at 20 °C(m2/s)
REAL (KIND = float) :: I, J !! terms to compute transport capacity in Yang method
REAL (KIND = float) :: conc !!sediment concentration (mg/l)

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

!compute channel unit stream power
vs = v * s

!compute channel critical unit stream power
vsCrit = CriticalVelocity(d, dp, s) * s

!compute I and J

I = 5.165 - 0.153 * LOG (FallVelocity(d) * (d/1000.) / kVisc) - &
    0.297 * LOG (ShearVelocity(dp,s)/FallVelocity(d)) 
    
J = 1.780 - 0.360 * LOG (FallVelocity(d) * (d/1000.) / kVisc) - &
    0.480 * LOG (ShearVelocity(dp,s)/FallVelocity(d)) 
    
!compute sediment concentration (ppm or mg/l)
IF (vs > vsCrit) THEN
   conc = EXP ( I + J * LOG((vs - vsCrit)/FallVelocity(d)) )
ELSE
   conc = 0.
END IF

IF (conc < 0.) THEN
  conc = 0.
END IF

!compute transform capacity (kg/s)
tc = flow * conc / 1000.


END FUNCTION Yang1979


!============================================================================
!| Description:
!  compute fall velocity of sediment (m/s)
!
! References:
!
!  Wikipedia contributors, "Sediment transport," Wikipedia, 
!     The Free Encyclopedia, http://en.wikipedia.org/w/index.php?title=Sediment_transport&oldid=449390416 
!    (accessed October 29, 2011).
FUNCTION FallVelocity &
!
(d) &
!
RESULT (fv)

IMPLICIT NONE

REAL (KIND = float), INTENT(IN) :: d !!particle size of bed material (mm)

!local declarations:
REAL (KIND = float) :: fv !!fall velocity (m/s)

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

fv = (16.17 * (d/1000.)**2.0) / (0.000018 + (12.1275 * (d/1000.)**3.)**0.5)


END FUNCTION FallVelocity


!============================================================================
!| Description:
!  compute shear velocity (m/s)
FUNCTION ShearVelocity &
!
(dp, s) &
!
RESULT (sv)

IMPLICIT NONE

REAL (KIND = float), INTENT(IN) :: dp !!water depth (m)
REAL (KIND = float), INTENT(IN) :: s !!channel slope (m/m)

!local declarations:
REAL (KIND = float) :: sv !!shear velocity (m/s)

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

sv = (gravityAccel * dp * s)**0.5


END FUNCTION ShearVelocity


!============================================================================
!| Description:
!  compute critical velocity for incipient motion (m/s)
!
! References:
!
!  Yang, C.H., Stall, J.B., Unit stream power for sediment transport
!    in natural rivers, University of Illinois Water Resources
!    Center Research Report No. 88, 1974
!    url: http://web.extension.illinois.edu/iwrc/pdf/88.pdf
FUNCTION CriticalVelocity &
!
(d, dp, s) &
!
RESULT (cv)

IMPLICIT NONE

!Arguments with intent in:
REAL (KIND = float), INTENT(IN) :: d !!particle size of bed material (mm)
REAL (KIND = float), INTENT(IN) :: dp !!water depth (m)
REAL (KIND = float), INTENT(IN) :: s !!channel slope (m/m)

!local declarations:
REAL (KIND = float) :: cv !!critical velocity (m/s)
REAL (KIND = float) :: kVisc = 0.000001004 !!kinematic viscosity of water at 20 °C(m2/s)
REAL (KIND = float) :: Restar


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

Restar = ShearVelocity(dp,s) * (d/1000.) / kVisc

IF (Restar < 70.) THEN
   
   cv = FallVelocity(d) * ( 2.5 / (LOG (Restar) - 0.06) + 0.66 )

ELSE

   cv = FallVelocity(d) * 2.05

END IF



END FUNCTION CriticalVelocity




!==============================================================================
!! Description:
!   compute slope factor
SUBROUTINE ComputeSlopeFactor &
!
(slope)

IMPLICIT NONE

!Arguments with intent in
TYPE (grid_real), INTENT (IN) :: slope

!local declarations:
INTEGER :: i,j

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

CALL NewGrid (slopeFactor, slope)

DO j = 1, slope % jdim
  DO i = 1, slope % idim
    IF (slope % mat (i,j) /= slope % nodata) THEN
       slopeFactor % mat (i,j) = 1.05 - 0.85 * EXP( -4.0 * SIN(slope % mat(i,j)) )
       IF ( slopeFactor % mat (i,j) < 0.2) THEN
         !filter negative value due to negative slope in pits
         slopeFactor % mat (i,j) = 0.2
       END IF
    END IF
  END DO
END DO

RETURN
END SUBROUTINE ComputeSlopeFactor


!==============================================================================
!| Description:
!  compute sediment routing in channelized drainage network
!  using Muskingum-Cunge method.
!  Interrill erosion is supposed to act as suspended sediment so
!  deposition is not contemplated. Bedload sediment is divided in
!  different grainsize classes and routing is performed for 
!  each class separately. Deposition is contemplated comparing 
!  sediment discharge to transport capacity. 
!  Flow Routing method code remind:
!     0  --> hillslope, water flow rate is not defined
!     1,2,3,11,30  --> channel routing, water flow is computed according 
!                      to different schemes
!     5  --> Instantaneous mass transferring inside lakes (infinitive celerity)
!     1000:2000 --> reservoir
!
! References:
!
!  Sun, H., P. S. Cornish, and T. M. Daniell, Contour-based digital 
!  elevation modeling of watershed erosion and sedimentation:
!  Erosion and sedimentation estimation tool (EROSET), 
!  Water Resour. Res., 38(11), 1233, doi:10.1029/2001WR000960, 2002. 
!
!  Singh, V.P., Quiroga, C.A., A dam-breach erosion model: I. formulation,
!  Water resources management, 1, 177-197, 1987. 
SUBROUTINE SedimentRouting &
!
(dt, flow, v, dp)


!Arguments with intent in:
INTEGER (KIND = short) :: dt !!routing time step
TYPE (grid_real), INTENT(IN) :: flow
TYPE (grid_real), INTENT(IN) :: v !!water flow velocity in channel (m/s)
TYPE (grid_real), INTENT(IN) :: dp !!water depth (m)

!local declarations:
INTEGER (KIND = short) :: k, iin, jin, is, js
REAL (KIND = float) :: ddx, cappa
REAL (KIND = float) :: vel
REAL (KIND = float) :: transportCapacity
REAL (KIND = float) :: wd
REAL (KIND = float) :: width

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

QinSS = 0.

DO K=1,sedReach%nreach
  iin=sedReach%branch(k)%i0
  jin=sedReach%branch(k)%j0
  DO WHILE ( .NOT.((jin == sedReach%branch(k)%j1) .AND. &
	               (iin == sedReach%branch(k)%i1)) )
	               
     CALL DownstreamCell (iin, jin, sedFlowDirection % mat(iin,jin), &
                          is, js, ddx, sedFlowDirection)
     
     SELECT CASE (sedReach % branch (k) % routingMethod)
     
       CASE (0)
          !compute cumulated sediment rate on hillslope along drainage network 
          !first cell of order 1 reach
          IF (sedReach % branch (k) % order == 1 .AND. &
              sedReach % branch (k) % I0 == iin   .AND. &
              sedReach % branch (k) % J0 == jin) THEN
             QinSS % mat(iin,jin) = interrillErosion % mat(iin,jin)
             QinSS % mat(is,js) =  QinSS % mat(iin,jin) + interrillErosion % mat(is,js)
          ELSE
             QinSS % mat(is,js) = QinSS % mat(is,js) + QinSS % mat(iin,jin) +  interrillErosion % mat(is,js)
          END IF
		  !QinSS % mat(is,js) = QinSS % mat(is,js) + QinSS % mat(iin,jin)
		  
		  !update storage sediment variation. On hillslope variation is always negative
		  deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) - &
		                             interrillErosion % mat(iin,jin) * dt / 1000000.
		                             

       CASE (1,2,3,11,30)
          ! sediment routing in channel reach  
          
          !suspended solid
          !^^^^^^^^^^^^^^^
!          QoutSS % mat(iin,jin) = QinSS % mat(iin,jin)  * gridC1 % mat(iin,jin) + &
!                                  PinSS % mat(iin,jin)  * gridC2 % mat(iin,jin) + &
!                                  PoutSS % mat(iin,jin) * gridC3 % mat(iin,jin)
      
           !DA SISTEMARE - GR 12/10/2016
           !vel = flow % mat (iin,jin) / (dp % mat (iin,jin) * dp % mat (iin,jin) * sedReach % branch (k) % B0) 
          
           IF (vel > 0.) THEN
 
              cappa = ddx / vel
              !cappa = 3
              
              CALL MuskingumCungeTodini ( dt = dt, dx = ddx, &
                                 Qin =  QinSS%mat(iin,jin), &
                                 Pin =  PinSS%mat(iin,jin), &
                                 Pout =  PoutSS%mat(iin,jin), &
                                 Qlat = 0., Plat = 0., &
                                 B0 = 0., alpha = 0., slope = 0., n = 0., &
                                 Qout = QoutSS%mat(iin,jin), &
                                 depth = wd, width =  width, k = cappa, x = 0. )
       
  
           ELSE
              QoutSS % mat(iin,jin) = 0.
           END IF

          !filter negative discharge                                            
          IF ( QoutSS % mat(iin,jin) < 0.) THEN
             QoutSS % mat(iin,jin) = 0.
          END IF
          
           !limit suspended load to transport capacity
           !assume particle size coming from hillslope of 0.02 mm
          transportCapacity = Yang1979 (flow % mat(iin,jin),       &
                                        sedReach % branch(k) % slope, &
                                        vel, 0.02, &
                                        dp % mat(iin,jin)          )
          IF ( QoutSS % mat(iin,jin) > transportCapacity ) THEN
            
!             deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
!		                            (QoutSS % mat(iin,jin) - transportCapacity) * dt / 1000000.
             QoutSS % mat(iin,jin) = transportCapacity
             QinSS % mat(iin,jin) =  transportCapacity
             
          END IF
          
           !update storage sediment variation
          deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
		                            (QinSS % mat(iin,jin) - QoutSS % mat(iin,jin)) * dt / 1000000.
          
          
          !update downstream input discharge
          QinSS % mat(is,js) =  QinSS % mat(is,js) + QoutSS % mat(iin,jin)
          !store computed discharge for next time step
          PinSS % mat(iin,jin) = QinSS % mat(iin,jin)
          PoutSS % mat(iin,jin) = QoutSS % mat(iin,jin)
          
   
          !bed load transport
          !^^^^^^^^^^^^^^^^^^
          !compute source term
          QinBL % mat(iin,jin) = QinBL % mat(iin,jin) + ChannelDetachmentRate (flow % mat(iin,jin), &
                                                  sedReach%branch(k)%slope, &
                                                  rusleK % mat (iin,jin),    &
                                                  rusleC % mat (iin,jin)     )
                                                  
          !update storage sediment variation
!		  deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) - &
!		                             QinBL % mat(iin,jin) * dt / 1000000.
                                                                 
          !route downstream                                        
!          QoutBL % mat(iin,jin) = QinBL % mat(iin,jin)  * gridC1 % mat(iin,jin) + &
!                              PinBL % mat(iin,jin)  * gridC2 % mat(iin,jin) + &
!                              PoutBL % mat(iin,jin) * gridC3 % mat(iin,jin)
                              
            IF (vel > 0.) THEN
              cappa = ddx / vel
              
              CALL MuskingumCungeTodini ( dt = dt, dx = ddx, &
                                 Qin =  QinBL%mat(iin,jin), &
                                 Pin =  PinBL%mat(iin,jin), &
                                 Pout =  PoutBL%mat(iin,jin), &
                                 Qlat = 0., Plat = 0., &
                                 B0 = 0., alpha = 0., slope = 0., n = 0., &
                                 Qout = QoutBL%mat(iin,jin), &
                                 depth = wd, width = width, k = cappa, x = 0. )
  
           ELSE
              QoutBL % mat(iin,jin) = 0.
           END IF                   
                              
                              
                              
          !filter negative discharge                                            
          IF ( QoutBL % mat(iin,jin) < 0.) THEN
             QoutBL % mat(iin,jin) = 0.
          END IF 
                                
          !limit bed load to transport capacity
          transportCapacity = Yang1973 (flow % mat(iin,jin),       &
                                        sedReach%branch(k)%slope, &
                                        vel, sedReach%branch(k)%d50, &
                                        dp % mat(iin,jin)          )
          IF ( QoutBL % mat(iin,jin) > transportCapacity ) THEN
             !update storage sediment variation if deposition is positive
!             deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
!                     (QoutBL % mat(iin,jin) - transportCapacity) * dt / 1000000.
             QoutBL % mat(iin,jin) = transportCapacity
             QinBL % mat(iin,jin) = transportCapacity
             
          END IF
          !update storage sediment variation 
           deltaSed % mat(iin,jin) =  deltaSed % mat(iin,jin) + &
		                            (QinBL % mat(iin,jin) - QoutBL % mat(iin,jin)) * dt / 1000000.
          
                                                            
          !update downstream input discharge
          QinBL % mat(is,js) =  QinBL % mat(is,js) + QoutBL % mat(iin,jin)
          !store computed discharge for next time step
          PinBL % mat(iin,jin) = QinBL % mat(iin,jin)
          PoutBL % mat(iin,jin) = QoutBL % mat(iin,jin)
          
          !build QoutSed
          !^^^^^^^^^^^^^
          QoutSed  % mat(iin,jin) = QoutSS % mat(iin,jin)  + QoutBL % mat(iin,jin)
          !QoutSed  % mat(iin,jin) = QoutBL % mat(iin,jin)
          
          
     END SELECT
     
     !go downstream 
     jin = js
     iin = is
  END DO
  !last cell of last reach
  IF (K == sedReach%nreach) THEN
      
  END IF
END DO


END SUBROUTINE SedimentRouting






END MODULE Sediment