LongestFlowPathSlope Function

public function LongestFlowPathSlope(fdir, dem, x, y, outsection, headsection) result(slope)

compute slope of longest flow path (m/m)

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: fdir
type(grid_real), intent(in) :: dem
real(kind=float), intent(in) :: x
real(kind=float), intent(in) :: y
type(Coordinate), intent(out), optional :: outsection

coordinate of outlet section
type(Coordinate), intent(out), optional :: headsection

coordinate of head section

Return Value real(kind=float)


Variables

Type Visibility Attributes Name Initial
type(grid_integer), public :: basin
integer(kind=short), public :: col

current cell
integer(kind=short), public :: i
integer(kind=short), public :: iDown

downstream cell
integer(kind=short), public :: iHead

coordinate of the head of the longest flow path
integer(kind=short), public :: iOutlet

coordinate of the basin outlet
integer(kind=short), public :: iSpring

coordinate of the spring
integer(kind=short), public :: j
integer(kind=short), public :: jDown

downstream cell
integer(kind=short), public :: jHead

coordinate of the head of the longest flow path
integer(kind=short), public :: jOutlet

coordinate of the basin outlet
integer(kind=short), public :: jSpring

coordinate of the spring
real(kind=float), public :: length
real(kind=float), public :: lmax
logical, public :: outlet
integer(kind=short), public :: row

current cell
real(kind=float), public :: xd
real(kind=float), public :: xu
real(kind=float), public :: yd
real(kind=float), public :: yu

Source Code

FUNCTION LongestFlowPathSlope &
!
(fdir, dem, x, y, outsection, headsection) &
!
RESULT (slope)

IMPLICIT NONE

!Arguments with intent (in)
TYPE(grid_integer),INTENT(IN) :: fdir !flow direction
TYPE(grid_real),INTENT(IN) :: dem !digital elevation model
REAL (KIND = float), INTENT(in) :: x, y  

!Arguments with intent out
TYPE (Coordinate), OPTIONAL, INTENT (OUT) :: outsection  !!coordinate of outlet section
TYPE (Coordinate), OPTIONAL, INTENT (OUT) :: headsection  !!coordinate of head section

!local declarations
REAL (KIND = float) :: slope
REAL (KIND = float) :: lmax
REAL (KIND = float) :: length
REAL (KIND = float) :: xu, yu, xd, yd
TYPE(grid_integer) :: basin
INTEGER (KIND = short) :: row, col !!current cell
INTEGER (KIND = short) :: iDown, jDown !!downstream cell
INTEGER (KIND = short) :: iSpring, jSpring !!coordinate of the spring
INTEGER (KIND = short) :: iHead, jHead !!coordinate of the head of the longest flow path
INTEGER (KIND = short) :: iOutlet, jOutlet !!coordinate of the basin outlet
INTEGER (KIND = short) :: i, j
LOGICAL :: outlet

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

!delineate river basin
CALL BasinDelineate (fdir, x, y, basin)

!overlay flow direction map
CALL GridResample (fdir, basin)

!find local coordinate of basin outlet
CALL GetIJ (x, y, fdir, iOutlet, jOutlet)


point1 % system = basin % grid_mapping
point2 % system = basin % grid_mapping
lmax = 0.
iHead = 0
jHead = 0
!loop trough basin
DO j = 1,basin % jdim
  DO i = 1,basin % idim
    IF (basin % mat (i,j) /= basin % nodata) THEN
        
        IF(CellIsSpring(i,j,basin)) THEN !found a spring
             
              length = 0.
              iSpring = i
              jSpring = j

              !follow the reach till basin outlet
               row                = i
               col                = j
               outlet             = .FALSE.
              
           DO WHILE (.NOT. outlet) ! follow the reach till the basin outlet 
    	                                                            
              CALL DownstreamCell(row, col, basin%mat(row,col), iDown, jDown)                         
              
              CALL GetXY (row,col,basin,xu,yu)
              CALL GetXY (iDown,jDown,basin,xd,yd)
              point1 % northing = yu  
              point1 % easting = xu
              point2 % northing = yd  
              point2 % easting = xd

              length = length + Distance(point1,point2)
              
              outlet = CheckOutlet(row,col,iDown,jDown,basin)
              IF (outlet) THEN
                 IF (length > lmax) THEN
                    !update lmax
                    lmax = length
                    !update head 
                    iHead = iSpring
                    jHead = jSpring
                 END IF
              END IF
              
              !loop
              row = iDown
              col = jDown

           END DO
                      
        ENDIF
    END IF
  ENDDO
ENDDO 

!compute slope (m/m)
slope = ( dem % mat(iHead,jHead) - dem % mat (iOutlet,jOutlet) ) / lmax

IF (PRESENT(outsection)) THEN
    outsection % system = basin % grid_mapping
    CALL GetXY (iOutlet, jOutlet, basin, outsection % easting, outsection % northing)
END IF

IF (PRESENT(headsection)) THEN
    headsection % system = basin % grid_mapping
    CALL GetXY (iHead, jHead, basin, headsection % easting, headsection % northing)
END IF

RETURN
END FUNCTION LongestFlowPathSlope