Sediment.f90 Source File
Compute sediment erosion and deposition
Source Code
!! 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