!! Module to deal with river and basin morphology
!|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 |
! |----------|--------------|----------|
! | 0.1      |  17/Aug/2009 |  Original code |
! | 0.2      |  23/Oct/2010 |  Compute flow accumulation grid |
! | 0.3      |  04/Dec/2012 |  Compute longest flow length |
! | 0.4      |  03/Apr/2013 |  Derive aspect subroutine |
! | 0.5      |  10/Apr/2013 |  Derive curvature subroutine |
! | 0.6      |  13/Apr/2013 |  Fixed DeriveSlope to match ArcGis algorithm |
! | 0.7      |  16/Apr/2013 |  DrainageDensity rooutine |
! | 0.8      |  03/May/2013 |  TopographicIndex routine |
! | 0.9      |  10/Jun/2013 |  CurvatureMicroMet routine |
! | 1.0      |  09/Oct/2013 |  LongestFlowPathSlope |
! | 1.1      |  03/Feb/2021 |  Compute Centroid  |
! | 1.2      |  22/Apr/2024 |  ESRI and GRASS flow direction can be used  |
!
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
! This file is part of 
!
!   MOSAICO -- MOdular library for raSter bAsed hydrologIcal appliCatiOn.
! 
!   Copyright (C) 2011 Giovanni Ravazzani
!
!### Code Description 
!   Language:           Fortran 90. 
!
!   Software Standards: "European Standards for Writing and  
!   Documenting Exchangeable Fortran 90 Code". 
! 
!### Module Description
!   library to deal with river and basin morphology
MODULE Morphology         
			

! Modules used: 

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

USE GridLib, ONLY: &
    !Imported type definitions:
    grid_integer, &
    grid_real, &
    ! Imported routines:
    NewGrid, &
    GridDestroy, &
    ExportGrid, &
    !Imported parameters:
    ESRI_ASCII, &
    ESRI_BINARY, &
    NET_CDF

USE GridOperations, ONLY: &
    !Imported routines:
    IsOutOfGrid, &
    GetIJ, &
    CellArea, &
    GetXY, &
    GridResample, &
    ExtractBorder, &
    CellLength, &
    ASSIGNMENT(=)

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


USE LogLib, ONLY: &
    !Imported routines:
    Catch

USE Units, ONLY: &
    !Imported parameters:
    pi

USE StringManipulation, ONLY: &
    !Imported routines:
    ToString, &
    StringToLower


IMPLICIT NONE 

PUBLIC :: HortonOrders
PUBLIC :: CheckOutlet
PUBLIC :: DownstreamCell
PUBLIC :: FlowAccumulation
PUBLIC :: DeriveSlope
PUBLIC :: DeriveAspect
PUBLIC :: DeriveCurvature
PUBLIC :: CurvatureMicroMet
PUBLIC :: BasinDelineate
PUBLIC :: CellIsSpring
PUBLIC :: LongestFlowLength
PUBLIC :: LongestFlowPathSlope
PUBLIC :: Perimeter
PUBLIC :: Centroid
PUBLIC :: SetFlowDirectionConvention

PRIVATE :: ConfluenceIsAround
PRIVATE :: BasinMask
PRIVATE :: BasinArea
PRIVATE :: CentroidGridReal
PRIVATE :: CentroidGridInteger

INTERFACE Centroid
   MODULE PROCEDURE CentroidGridReal
   MODULE PROCEDURE CentroidGridInteger
END INTERFACE

!flow direction
INTEGER (KIND = short) :: NW !! North-West
INTEGER (KIND = short) :: N  !! North 
INTEGER (KIND = short) :: NE !! North-East
INTEGER (KIND = short) :: W  !! West 
INTEGER (KIND = short) :: E  !! East 
INTEGER (KIND = short) :: SW !! South-West
INTEGER (KIND = short) :: S  !! South 
INTEGER (KIND = short) :: SE !! South-East

!flow direction conventions
INTEGER (KIND = short), PARAMETER :: ESRI = 1, GRASS = 2
INTEGER (KIND = short) :: flowDirectionConvention

!=======         
CONTAINS
!======= 
! Define procedures contained in this module. 
    
!==============================================================================
!| Description:
!   Set flow direction convention:
!
!  -  ESRI: NW=32, N=64, NE=128, W=16, E=1, SW=8, S=4, SE=2
!  -  GRASS: NW=3, N=2, NE=1, W=4, E=8, SW=5, S=6, SE=7
!
SUBROUTINE SetFlowDirectionConvention &
!
(string)


IMPLICIT NONE

!Arguments with intent in:
CHARACTER (LEN = *), INTENT (IN) :: string

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

SELECT CASE ( TRIM (StringToLower(string) ) )
CASE ('esri')
    flowDirectionConvention = ESRI
    NW = 32; N = 64; NE = 128; W = 16; E = 1; SW = 8; S = 4; SE = 2
CASE ('grass')
    flowDirectionConvention = GRASS
    NW = 3; N = 2; NE = 1; W = 4; E = 8; SW = 5; S = 6; SE = 7
CASE DEFAULT
  CALL Catch ('error', 'Morphology', 'flow direction convention not found: ', &
             argument = TRIM (string) )
END SELECT

RETURN
END SUBROUTINE SetFlowDirectionConvention
    
    
    
!==============================================================================
!| Description:
!   returns the coordinate of the centroid given a `grid_integer` mask,
!   in the same coordinate reference system of `grid_integer`
SUBROUTINE CentroidGridInteger &
!
(grid, ccoord)


IMPLICIT NONE

!Arguments with intent in:
TYPE (grid_integer), INTENT (IN) :: grid  !!input grid


!Arguments with intent out or inout
TYPE (Coordinate), INTENT (OUT) :: ccoord  !! coordinate of computed centroid


!local declarations:
INTEGER (KIND = short) :: cj, ci !!coordinate 
INTEGER (KIND = short) :: i, j
INTEGER (KIND = short) :: countcell

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

!set coordinate reference system of centroid
ccoord % system = grid % grid_mapping

cj = 0
ci = 0
countcell = 0

DO i = 1, grid % idim
    DO j = 1, grid % jdim
        IF (grid % mat (i,j) /= grid % nodata) THEN
            countcell = countcell + 1
            cj = cj + j
            ci = ci + i
        END IF
    END DO
END DO

IF ( countcell > 0) THEN
    cj = cj / countcell
    ci = ci / countcell
ELSE
    CALL Catch ('error', 'Morphology', 'area zero found while computing centroid coordinates')
END IF

CALL GetXY (ci, cj, grid, ccoord % easting, ccoord % northing )

RETURN
END SUBROUTINE CentroidGridInteger

!==============================================================================
!| Description:
!   returns the coordinate of the centroid given a `grid_integer` mask,
!   in the same coordinate reference system of `grid_integer`
SUBROUTINE CentroidGridReal &
!
(grid, ccoord)


IMPLICIT NONE

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


!Arguments with intent out or inout
TYPE (Coordinate), INTENT (OUT) :: ccoord  !! coordinate of computed centroid


!local declarations:
INTEGER (KIND = short) :: cj, ci !!coordinate 
INTEGER (KIND = short) :: i, j
INTEGER (KIND = short) :: countcell

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

!set coordinate reference system to centroid
ccoord % system = grid % grid_mapping

cj = 0
ci = 0
countcell = 0

DO i = 1, grid % idim
    DO j = 1, grid % jdim
        IF (grid % mat (i,j) /= grid % nodata) THEN
            countcell = countcell + 1
            cj = cj + j
            ci = ci + i
        END IF
    END DO
END DO

IF ( countcell > 0) THEN
    cj = cj / countcell
    ci = ci / countcell
ELSE
    CALL Catch ('error', 'Morphology', 'area zero found while computing centroid coordinates')
END IF

CALL GetXY (ci, cj, grid, ccoord % easting, ccoord % northing )

RETURN
END SUBROUTINE CentroidGridReal

!==============================================================================
!| Description:
!   returns a grid_integer containing Horton orders. Horton orders are 
!   computed on the entire space-filled basin.
SUBROUTINE HortonOrders &
!
(flowDirection,orders,basinOrder)

IMPLICIT NONE

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


!Arguments with intent out or inout
TYPE (grid_integer), INTENT (INOUT) :: orders
INTEGER, OPTIONAL, INTENT (OUT) :: basinOrder !!the maximum order of the basin

!local declarations:
LOGICAL  :: confluence !!true if confluence
LOGICAL  :: outlet  !!true if basin outlet

INTEGER  :: row, col !!current cell
INTEGER  :: iDown, jDown !!downstream cell
INTEGER  :: numConf !!number of confluences
INTEGER  :: order !! Horton order
INTEGER  :: cellsCount
INTEGER  :: i, j

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




order = 1
numConf = 1

DO WHILE (numConf > 0) ! se non trovo confluenze 
                                 ! di classe order
								 ! l'operazione è terminata

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

numConf = 0

!-----follow the reach till a confluence or a basin outlet------

DO j = 1,orders % jdim
  DO i = 1,orders % idim

    IF(CellIsSpring(i,j,flowDirection)) THEN !found a spring
           row                = i
           col                = j
           outlet             = .FALSE.
           confluence         = .FALSE.
           cellsCount         = 0
           orders % mat(i,j)  = 1
          
       DO WHILE (.NOT. outlet) ! follow the reach till the basin outlet 
	                                                            
          IF (orders % mat(row,col) == order ) THEN
              cellsCount = cellsCount + 1 			  
          ENDIF

          CALL DownstreamCell(row, col, &
							  flowDirection%mat(row,col), &
                              iDown, jDown)       
          
          IF (cellsCount >= 1 ) THEN  !I am in the reach of that order
          
          !check if downstream cell is a confluence to increment horton order 
          !Downstream the confluence, till the basin outlet, as temptative value,
          !order is increased by 1 (order + 1)
             IF ( .NOT. confluence  ) THEN
                CALL ConfluenceIsAround (iDown, jDown, row, col, &
								    flowDirection,confluence,orders,order)
                IF(confluence) numConf = numConf + 1 
             ENDIF

			 outlet = CheckOutlet (row,col,iDown,jDown,flowDirection)
			 
             IF (.NOT. outlet) THEN
                IF (.NOT. confluence) THEN
                   orders % mat(iDown,jDown) = order
                ELSE
                   orders % mat(iDown,jDown) = order + 1        
                ENDIF
             ENDIF

          ENDIF ! cellsCount >= 1 

          outlet = CheckOutlet(row,col,iDown,jDown,flowDirection)
          
          !loop
          row = iDown
          col = jDown

       END DO
                  
    ENDIF

  ENDDO
ENDDO  !ciclo sulla matrice ordini
!------------------------------------------------------------------------------
order = order + 1

ENDDO 

IF ( PRESENT (basinOrder) ) THEN
    basinOrder = order - 1
END IF

END SUBROUTINE HortonOrders


!==============================================================================
!| Description:
!   Scan the eigth cells surrounding the center cell (row,col) (neglecting 
!   the upstream cell (i,j) ) to find if a confluence is present.
!   The confluence cell must be of the same order or not yet defined.
SUBROUTINE ConfluenceIsAround &
!
(row,col,i,j,flowDir,confluence,orders,Norder)


IMPLICIT NONE

TYPE (grid_integer), INTENT(IN):: flowDir
TYPE (grid_integer), INTENT(IN):: orders
INTEGER, INTENT(IN):: i,j,row,col,Norder
LOGICAL, INTENT(OUT):: confluence
!-----------------------------------------------------------


IF(.NOT. IsOutOfGrid(row,col+1,flowDir) ) THEN
    IF(flowDir%mat(row,col+1) == W .AND.&
       (row /= i .OR. col+1 /= j).AND. &
       (orders%mat(row,col+1) == orders%nodata.OR. &
       orders%mat(row,col+1) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row+1,col+1,flowDir) ) THEN
    IF(flowDir%mat(row+1,col+1) == NW .AND.&
       (row+1 /=i .OR. col+1 /= j).AND. &
       (orders%mat(row+1,col+1) == orders%nodata.OR. &
       orders%mat(row+1,col+1) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row+1,col,flowDir) ) THEN
    IF(flowDir%mat(row+1,col) == N .AND.&
       (row+1 /=i .OR. col /= j).AND. &
       (orders%mat(row+1,col) == orders%nodata.OR. &
       orders%mat(row+1,col) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row+1,col-1,flowDir) ) THEN
    IF(flowDir%mat(row+1,col-1) == NE .AND.&
       (row+1 /= i .OR. col-1 /= j).AND. &
       (orders%mat(row+1,col-1) == orders%nodata.OR. &
       orders%mat(row+1,col-1) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row,col-1,flowDir) ) THEN
    IF(flowDir%mat(row,col-1) == E .AND.&
       (row /= i .OR. col-1 /= j).AND. &
       (orders%mat(row,col-1) == orders%nodata.OR. &
       orders%mat(row,col-1) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row-1,col-1,flowDir) ) THEN
    IF(flowDir%mat(row-1,col-1) == SE .AND.&
       (row-1 /= i .OR. col-1 /= j).AND. &
       (orders%mat(row-1,col-1) == orders%nodata.OR. &
       orders%mat(row-1,col-1) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row-1,col,flowDir) ) THEN
    IF(flowDir%mat(row-1,col) == S .AND.&
       (row-1 /= i .OR. col /= j).AND. &
       (orders%mat(row-1,col) == orders%nodata.OR. &
       orders%mat(row-1,col) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row-1,col+1,flowDir) ) THEN
    IF(flowDir%mat(row-1,col+1) == SW .AND.&
       (row-1 /= i .OR. col+1 /= j).AND. &
       (orders%mat(row-1,col+1) == orders%nodata.OR. &
       orders%mat(row-1,col+1) == Norder) ) THEN
       confluence = .TRUE.
    ENDIF
ENDIF


RETURN
END SUBROUTINE ConfluenceIsAround

!==============================================================================
!| Description:
!   if the downstream cell is a nodata value or out of grid space, 
!   or points toward the current cell, the current cell
!   is the basin outlet. 
LOGICAL FUNCTION CheckOutlet &
!
(icurrent, jcurrent, iDown, jDown, flowDirection)


IMPLICIT NONE
INTEGER, INTENT(in)           :: icurrent
INTEGER, INTENT(in)           :: jcurrent
INTEGER, INTENT(in)           :: iDown !riga cella di valle
INTEGER, INTENT(in)           :: jDown !colonna cella di valle
TYPE(grid_integer), INTENT(in):: flowDirection

!local declarations:
INTEGER :: iback, jback
!------------------------------end of declaration -----------------------------

CheckOutlet = .FALSE.

!first case
IF ( IsOutOfGrid (iDown,jDown,flowDirection) ) THEN
	CheckOutlet = .TRUE.
    RETURN
END IF

!second case
IF ( flowDirection % mat(iDown,jDown) == flowDirection % nodata)  THEN
	CheckOutlet = .TRUE.
    RETURN
END IF


!third case
CALL DownstreamCell (iDown,jDown,flowDirection%mat(iDown,jDown),iback,jback)
IF ( iback == icurrent .AND. jback == jcurrent) THEN
    CheckOutlet = .TRUE.
    RETURN
END IF


RETURN
END FUNCTION CheckOutlet



!==============================================================================
!| Description:
!   returns the position (is,js) of the downstream cell and, optionally,
!   the flow path length, considering cardinal and diagonal direction
SUBROUTINE DownstreamCell &
!
(iin, jin, dir, is, js, dx, grid)

IMPLICIT NONE

!Arguments with intent in:
INTEGER (KIND = short), INTENT (IN) :: iin, jin !!current cell
INTEGER (KIND = long), INTENT (IN) :: dir !!flow direction 
TYPE(grid_integer), INTENT (IN),OPTIONAL:: grid !!used to define coordinate reference system

!Arguments with intent out:
INTEGER (KIND = short), INTENT (OUT) :: is,js !!downstream cell
REAL (KIND = float), INTENT (OUT) ,OPTIONAL :: dx !!flow path length [m]

!Local declarations
REAL (KIND = float) :: ddx,x,y,xs,ys

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

!east
IF ( dir == E ) THEN
    js = jin + 1
    is = iin
END IF

!south east
IF ( dir == SE ) THEN
    js = jin + 1
    is = iin + 1
END IF

!south
IF ( dir == S ) THEN
    js = jin
    is = iin + 1
END IF

!south west
IF ( dir == SW ) THEN
    js = jin - 1
    is = iin + 1
END IF

!west
IF ( dir == W ) THEN
    js = jin - 1
    is = iin
END IF

!north west
IF ( dir == NW ) THEN
    js = jin - 1
    is = iin - 1
END IF

!north
IF ( dir == N ) THEN
    js = jin
    is = iin - 1
END IF

!north easth
IF ( dir == NE ) THEN
    js = jin + 1
    is = iin - 1
END IF

IF ((PRESENT(dx)).and.(PRESENT(grid))) THEN
    CALL GetXY (iin,jin,grid,x,y)
    CALL GetXY (is,js,grid,xs,ys)
    
    point1 % system = grid % grid_mapping
    point2 % system = grid % grid_mapping

    point1 % northing = y  
    point1 % easting = x
    point2 % northing = ys  
    point2 % easting = xs

    ddx = distance(point1,point2)
    dx = ddx
    RETURN

ELSE IF ((PRESENT(dx)) .AND. .NOT. (PRESENT(grid))) THEN
   CALL Catch ('error', 'Morphology', 'missing grid while calculating downstream distance')
END IF

END SUBROUTINE DownstreamCell



!==============================================================================
!| Description:
!   returns the position (is,js) of the upstream cell and, optionally,
!   the flow path length, considering cardinal and diagonal direction
SUBROUTINE UpstreamCell &
!
(i, j, flowdir, flowacc, is, js, found, dx)

IMPLICIT NONE

!Arguments with intent in:
INTEGER (KIND = short), INTENT (IN) :: i, j !!current cell
TYPE (grid_integer), INTENT (IN) :: flowdir  !!flow direction map
TYPE (grid_real), INTENT (IN) :: flowacc !! flow accumulation map

!Arguments with intent out:
INTEGER (KIND = short), INTENT (OUT) :: is,js !!upstream cell
LOGICAL, INTENT(OUT) :: found  !! upsteam cell found
REAL (KIND = float), INTENT (OUT) ,OPTIONAL :: dx !!flow path length [m]

!Local declarations
REAL (KIND = float) :: x,y,xs,ys
INTEGER (KIND = short) :: row, col
REAL (KIND = float) :: minDeltaArea, deltaArea


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


minDeltaArea = 1000000000.
is = 0
js = 0
found = .FALSE.

!east cell
row = i
col = j + 1
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == W ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF

!south-east cell
row = i + 1
col = j + 1
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == NW ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF

!south cell
row = i + 1
col = j 
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == N ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF


!south-west cell
row = i + 1
col = j - 1
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == NE ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF

!west cell
row = i
col = j - 1
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == E ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF


!north-west cell
row = i - 1
col = j - 1
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == SE ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF


!north cell
row = i - 1
col = j 
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == S ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF

!north-east cell
row = i - 1
col = j + 1
IF(.NOT. IsOutOfGrid(row,col,flowDir) ) THEN
    IF(flowDir%mat(row,col) == SW ) THEN
        !check delta area
        deltaArea = flowacc % mat (i, j) - flowacc % mat (row,col)
        
        IF ( deltaArea > 0. .AND. deltaArea < minDeltaArea ) THEN
            found = .TRUE.
            minDeltaArea = deltaArea
            is = row
            js = col
        END IF

    ENDIF
ENDIF

IF (found) THEN
    IF (PRESENT (dx) ) THEN
        CALL GetXY (i,j,flowDir,x,y)
        CALL GetXY (is,js,flowDir,xs,ys)
        
        point1 % system = flowDir % grid_mapping
        point2 % system = flowDir % grid_mapping
        
        point1 % northing = y  
        point1 % easting = x
        point2 % northing = ys  
        point2 % easting = xs

       dx = distance(point1,point2)

    END IF
END IF

END SUBROUTINE UpstreamCell



!==============================================================================
!| Description:
!   For each cell, the routine calculates the maximum rate of change in value 
!   from that cell to its neighbors. Basically, the maximum change in 
!   elevation over the distance between the cell and its eight neighbors 
!   identifies the steepest downhill descent from the cell.
!   Conceptually, the routine fits a plane to the z-values of a 3 x 3 cell 
!   neighborhood around the processing or center cell. The slope value of this 
!   plane is calculated using the average maximum technique. The direction 
!   the plane faces is the aspect for the processing cell.
!   If there is a cell location in the neighborhood with a NoData z-value, 
!   the z-value of the center cell will be assigned to the location.
!
! References:
!
!   Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical 
!       Information Systems (Oxford University Press, New York), 190 pp.
SUBROUTINE DeriveSlope &
!
(dem, slope)


IMPLICIT NONE

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

!arguments with intent out
TYPE (grid_real), INTENT (out) :: slope ![rad]

!local variables
INTEGER :: ii,jj
REAL (KIND = float) :: dZdX !!rate of change in the x direction
REAL (KIND = float) :: dZdY !!rate of change in the y direction
REAL (KIND = float) :: a !!elevation of N-W cell
REAL (KIND = float) :: b !!elevation of N cell
REAL (KIND = float) :: c !!elevation of N-E cell
REAL (KIND = float) :: d !!elevation of W cell
REAL (KIND = float) :: f !!elevation of E cell
REAL (KIND = float) :: g !!elevation of S-W cell
REAL (KIND = float) :: h !!elevation of S cell
REAL (KIND = float) :: i !!elevation of S-E cell
REAL (KIND = float) :: L !!length

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

!allocate slope grid
CALL NewGrid (slope,dem)

!loop through dem grid
idim: DO ii = 1, dem % idim
  jdim: DO jj = 1, dem % jdim
    IF (dem % mat(ii,jj) /= dem % nodata) THEN
    
        !set cell size
        L = CellArea (dem,ii,jj) ** 0.5
        
        !north-west cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj-1,dem) ) THEN
          IF ( dem % mat (ii-1,jj-1) /= dem % nodata) THEN
            a = dem % mat (ii-1,jj-1)
          ELSE
            a = dem % mat (ii,jj)
          END IF
        ELSE  
          a = dem % mat (ii,jj)  
        END IF
        
        !north cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj,dem) ) THEN
          IF ( dem % mat (ii-1,jj) /= dem % nodata) THEN
            b = dem % mat (ii-1,jj)
          ELSE
            b = dem % mat (ii,jj)
          END IF
        ELSE  
          b = dem % mat (ii,jj)  
        END IF
        
        !north-east cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj+1,dem) ) THEN
          IF ( dem % mat (ii-1,jj+1) /= dem % nodata) THEN
            c = dem % mat (ii-1,jj+1)
          ELSE
            c = dem % mat (ii,jj)
          END IF
        ELSE  
          c = dem % mat (ii,jj)  
        END IF
        
        !west cell
        IF ( .NOT. IsOutOfGrid (ii,jj-1,dem) ) THEN
          IF ( dem % mat (ii,jj-1) /= dem % nodata) THEN
            d = dem % mat (ii,jj-1)
          ELSE
            d = dem % mat (ii,jj)
          END IF
        ELSE  
          d = dem % mat (ii,jj)  
        END IF
        
        !east cell
        IF ( .NOT. IsOutOfGrid (ii,jj+1,dem) ) THEN
          IF ( dem % mat (ii,jj+1) /= dem % nodata) THEN
            f = dem % mat (ii,jj+1)
          ELSE
            f = dem % mat (ii,jj)
          END IF
        ELSE  
          f = dem % mat (ii,jj)  
        END IF
        
        !south-west cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj-1,dem) ) THEN
          IF ( dem % mat (ii+1,jj-1) /= dem % nodata) THEN
            g = dem % mat (ii+1,jj-1)
          ELSE
            g = dem % mat (ii,jj)
          END IF
        ELSE  
          g = dem % mat (ii,jj)  
        END IF
        
        !south cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj,dem) ) THEN
          IF ( dem % mat (ii+1,jj) /= dem % nodata) THEN
            h = dem % mat (ii+1,jj)
          ELSE
            h = dem % mat (ii,jj)
          END IF
        ELSE  
          h = dem % mat (ii,jj)  
        END IF
        
        !south-east cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj+1,dem) ) THEN
          IF ( dem % mat (ii+1,jj+1) /= dem % nodata) THEN
            i = dem % mat (ii+1,jj+1)
          ELSE
            i = dem % mat (ii,jj)
          END IF
        ELSE  
          i = dem % mat (ii,jj)  
        END IF
        
        !compute the rate of change in the x direction
        dZdX = ((c + 2.*f + i) - (a + 2.*d + g)) / (8. * L)
        
        !compute the rate of change in the y direction
        
        dZdY = ((g + 2.*h + i) - (a + 2.*b + c))  / (8. * L)
        
        !compute slope
        slope % mat (ii,jj) = ATAN ( (dZdX ** 2. + dZdY ** 2.) ** 0.5 )
        
    
    END IF
  END DO jdim
END DO idim 


RETURN 
END SUBROUTINE DeriveSlope



!==============================================================================
!| Description:
! Aspect identifies the downslope direction of the maximum rate of change 
! in value from each cell to its neighbors. 
! Aspect can be thought of as the slope direction. The values of the output 
! raster will be the compass direction of the aspect.
! A moving window visits each cell in the input raster and for each cell 
! in the center of the window, an aspect value is calculated using an 
! algorithm that incorporates the values of the cell's eight neighbors. 
! If any neighborhood cells are NoData, they are first assigned the value 
! of the center cell, then the aspect is computed.
! The cells are identified as letters 'a' to 'i', with 'e' representing the 
! cell for which the aspect is being calculated.
! The rate of change in the x direction for cell 'e' is calculated with 
! the following algorithm:
!```
!    [dz/dx] = ((c + 2f + i) - (a + 2d + g)) / 8
!```
! The rate of change in the y direction for cell 'e' is calculated 
! with the following algorithm:
!```
!    [dz/dy] = ((g + 2h + i) - (a + 2b + c))  / 8
!```
! Taking the rate of change in both the x and y direction for 
! cell 'e', aspect is calculated using:
!```
!    aspect = 57.29578 * atan2 ([dz/dy], -[dz/dx])
!```
! The aspect value is then converted to compass direction values 
! (0–2pi rad), according to the following rule:
!```
!    if aspect < 0
!
!    cell = pi/2 - aspect
!
!    else if aspect > pi/2
!
!    cell = 2pi - aspect + pi/2
!
!    else
!
!    cell = pi/2 - aspect
!```
! Flat areas in the input raster—areas where the slope is zero are assigned 
! an aspect of -1.
!
! Reference:
!
! Burrough, P. A. and McDonell, R.A., 1998. Principles of Geographical 
!   Information Systems (Oxford University Press, New York), p. 190
SUBROUTINE DeriveAspect &
!
(dem, aspect)


IMPLICIT NONE

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

!arguments with intent out
TYPE (grid_real), INTENT (out) :: aspect !(rad)

!local variables
INTEGER :: ii,jj
REAL (KIND = float) :: dZdX !!rate of change in the x direction
REAL (KIND = float) :: dZdY !!rate of change in the y direction
REAL (KIND = float) :: aspectIJ
REAL (KIND = float) :: a !!elevation of N-W cell
REAL (KIND = float) :: b !!elevation of N cell
REAL (KIND = float) :: c !!elevation of N-E cell
REAL (KIND = float) :: d !!elevation of W cell
REAL (KIND = float) :: f !!elevation of E cell
REAL (KIND = float) :: g !!elevation of S-W cell
REAL (KIND = float) :: h !!elevation of S cell
REAL (KIND = float) :: i !!elevation of S-E cell

!------------------------------end of declaration -----------------------------
!allocate aspect grid
CALL NewGrid (aspect,dem)

!loop through dem grid
idim: DO ii = 1, dem % idim
  jdim: DO jj = 1, dem % jdim
    IF (dem % mat(ii,jj) /= dem % nodata) THEN
        
        !north-west cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj-1,dem) ) THEN
          IF ( dem % mat (ii-1,jj-1) /= dem % nodata) THEN
            a = dem % mat (ii-1,jj-1)
          ELSE
            a = dem % mat (ii,jj)
          END IF
        ELSE  
          a = dem % mat (ii,jj)  
        END IF
        
        !north cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj,dem) ) THEN
          IF ( dem % mat (ii-1,jj) /= dem % nodata) THEN
            b = dem % mat (ii-1,jj)
          ELSE
            b = dem % mat (ii,jj)
          END IF
        ELSE  
          b = dem % mat (ii,jj)  
        END IF
        
        !north-east cell
        IF ( .NOT. IsOutOfGrid (ii-1,jj+1,dem) ) THEN
          IF ( dem % mat (ii-1,jj+1) /= dem % nodata) THEN
            c = dem % mat (ii-1,jj+1)
          ELSE
            c = dem % mat (ii,jj)
          END IF
        ELSE  
          c = dem % mat (ii,jj)  
        END IF
        
        !west cell
        IF ( .NOT. IsOutOfGrid (ii,jj-1,dem) ) THEN
          IF ( dem % mat (ii,jj-1) /= dem % nodata) THEN
            d = dem % mat (ii,jj-1)
          ELSE
            d = dem % mat (ii,jj)
          END IF
        ELSE  
          d = dem % mat (ii,jj)  
        END IF
        
        !east cell
        IF ( .NOT. IsOutOfGrid (ii,jj+1,dem) ) THEN
          IF ( dem % mat (ii,jj+1) /= dem % nodata) THEN
            f = dem % mat (ii,jj+1)
          ELSE
            f = dem % mat (ii,jj)
          END IF
        ELSE  
          f = dem % mat (ii,jj)  
        END IF
        
        !south-west cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj-1,dem) ) THEN
          IF ( dem % mat (ii+1,jj-1) /= dem % nodata) THEN
            g = dem % mat (ii+1,jj-1)
          ELSE
            g = dem % mat (ii,jj)
          END IF
        ELSE  
          g = dem % mat (ii,jj)  
        END IF
        
        !south cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj,dem) ) THEN
          IF ( dem % mat (ii+1,jj) /= dem % nodata) THEN
            h = dem % mat (ii+1,jj)
          ELSE
            h = dem % mat (ii,jj)
          END IF
        ELSE  
          h = dem % mat (ii,jj)  
        END IF
        
        !south-east cell
        IF ( .NOT. IsOutOfGrid (ii+1,jj+1,dem) ) THEN
          IF ( dem % mat (ii+1,jj+1) /= dem % nodata) THEN
            i = dem % mat (ii+1,jj+1)
          ELSE
            i = dem % mat (ii,jj)
          END IF
        ELSE  
          i = dem % mat (ii,jj)  
        END IF
        
        !compute the rate of change in the x direction
        dZdX = ((c + 2.*f + i) - (a + 2.*d + g)) / 8.
        
        !compute the rate of change in the y direction
        
        dZdY = ((g + 2.*h + i) - (a + 2.*b + c))  / 8.
        
        !compute aspect
        IF (dZdX == 0. .AND. dZdY == 0.) THEN
            aspect % mat (ii,jj)  = -1.
            CYCLE
        ELSE
            aspectIJ  =  ATAN2 (dZdY, -dZdX)
        END IF
        
        
        !convert aspect value to compass direction values (0–2Pi radians)
        
        IF (aspectIJ < 0) THEN

            aspect % mat (ii,jj)  = pi/2. - aspectIJ

        ELSE IF (aspectIJ > pi /2.) THEN 

           aspect % mat (ii,jj) = 2.*pi - aspectIJ + pi/2.
           
        ELSE

            aspect % mat (ii,jj) = pi/2. - aspectIJ 
           
        END IF
       
    END IF
  END DO jdim
END DO idim 

RETURN

END SUBROUTINE DeriveAspect


!==============================================================================
!| Description:
! Calculate curvature according to method implemented in MICROMET
!
! References:
!
!   Liston, G.E., Elder, K., A Meteorological Distribution System for 
!   High-Resolution Terrestrial Modeling, JOURNAL OF HYDROMETEOROLOGY,
!   7, 217-234, 2006.
SUBROUTINE CurvatureMicroMet &
!
(dem, curve_len_scale, curvature)

IMPLICIT NONE

!arguments with intent in
TYPE(grid_real),    INTENT(in) :: dem
REAL (KIND = float), INTENT(in) :: curve_len_scale

!arguments with intent out
TYPE (grid_real), INTENT (out) :: curvature

!local declarations:
 REAL (KIND = float) :: deltaxy
 REAL (KIND = float) :: z, zw, ze, zs, zn, zsw, zne, znw, zse
 REAL (KIND = float) :: curve_max
 INTEGER (KIND = short) :: inc
 INTEGER (KIND = short) :: i, j


!-----------------------------end of declarations------------------------------
!allocate grid
CALL NewGrid (curvature, dem)


!compute curvature
DO i = 1, dem % idim
  DO j = 1, dem % jdim
   IF (dem % mat (i,j) /= dem % nodata) THEN
   
     !average grid increment
     deltaxy = CellArea (dem,i,j) ** 0.5

     !Convert the length scale to an appropriate grid increment.
     inc = max(1,nint(curve_len_scale/deltaxy))
     
     !srt z
     z = dem % mat(i,j)
     
     !set zn
     IF ( .NOT. IsOutOfGrid (i-inc,j,dem) ) THEN
        IF ( dem % mat (i-inc,j) /= dem % nodata) THEN
           zn = dem % mat (i-inc,j)
        ELSE
           zn = dem % mat (i,j)
        END IF  
      ELSE 
           zn = dem % mat (i,j)
      END IF
      
      !set zne
      IF ( .NOT. IsOutOfGrid (i-inc,j+inc,dem) ) THEN
         IF ( dem % mat (i-inc,j+inc) /= dem % nodata) THEN
            zne = dem % mat (i-inc,j+inc)
         ELSE
            zne = dem % mat (i,j)
         END IF  
       ELSE 
            zne = dem % mat (i,j)
       END IF
       
       !set ze
       IF ( .NOT. IsOutOfGrid (i,j+inc,dem) ) THEN
          IF ( dem % mat (i,j+inc) /= dem % nodata) THEN
             ze = dem % mat (i,j+inc)
          ELSE
             ze = dem % mat (i,j)
          END IF  
        ELSE 
             ze = dem % mat (i,j)
        END IF
        
        !set zse
        IF ( .NOT. IsOutOfGrid (i+inc,j+inc,dem) ) THEN
           IF ( dem % mat (i+inc,j+inc) /= dem % nodata) THEN
              zse = dem % mat (i+inc,j+inc)
           ELSE
              zse = dem % mat (i,j)
           END IF  
         ELSE 
              zse = dem % mat (i,j)
         END IF
         
         !set zs
         IF ( .NOT. IsOutOfGrid (i+inc,j,dem) ) THEN
            IF ( dem % mat (i+inc,j) /= dem % nodata) THEN
               zs = dem % mat (i+inc,j)
            ELSE
               zs = dem % mat (i,j)
            END IF  
          ELSE 
               zs = dem % mat (i,j)
          END IF
          
          !set zsw
          IF ( .NOT. IsOutOfGrid (i+inc,j-inc,dem) ) THEN
             IF ( dem % mat (i+inc,j-inc) /= dem % nodata) THEN
                zsw = dem % mat (i+inc,j-inc)
             ELSE
                zsw = dem % mat (i,j)
             END IF  
           ELSE 
                zsw = dem % mat (i,j)
           END IF
           
           !set zw
           IF ( .NOT. IsOutOfGrid (i,j-inc,dem) ) THEN
              IF ( dem % mat (i,j-inc) /= dem % nodata) THEN
                 zw = dem % mat (i,j-inc)
              ELSE
                 zw = dem % mat (i,j)
              END IF  
            ELSE 
                 zw = dem % mat (i,j)
            END IF
            
            !set znw
            IF ( .NOT. IsOutOfGrid (i-inc,j-inc,dem) ) THEN
               IF ( dem % mat (i-inc,j-inc) /= dem % nodata) THEN
                  znw = dem % mat (i-inc,j-inc)
               ELSE
                  znw = dem % mat (i,j)
               END IF  
             ELSE 
                  znw = dem % mat (i,j)
             END IF
             
             
             !curvature
              curvature % mat(i,j) = (4.0 * z - znw - zse - zsw - zne ) / &
                   (sqrt(2.0) * 16.0 * real(inc) * deltaxy) + &
                   (4.0 * z - zs - zn - ze - zw) / &
                   (16.0 * real(inc) * deltaxy)
      
   END IF
  END DO
END DO


! Scale the curvature such that the max abs(curvature) has a value
!   of abs(0.5).  Include a 1 mm curvature in curve_max to prevent
!   divisions by zero in flat terrain where the curvature is zero.

curve_max = 0.0 + 0.001

DO i = 1, dem % idim
  DO j = 1, dem % jdim
    IF (dem % mat (i,j) /= dem % nodata) THEN 
       curve_max = MAX(curve_max,ABS(curvature%mat(i,j)))
    END IF
  END DO
END DO

DO i = 1, dem % idim
  DO j = 1, dem % jdim
    IF (dem % mat (i,j) /= dem % nodata) THEN 
       curvature%mat(i,j) = curvature%mat(i,j) / (2.0 * curve_max)
    END IF
  END DO
END DO
       

RETURN

END SUBROUTINE CurvatureMicroMet

!==============================================================================
!| Description:
! Calculates the curvature [1/m] of a raster surface, optionally including 
! profile and plan curvature. If any neighborhood cells are NoData, they 
! are first assigned the value of the center cell.
! It calculates the second derivative value of the input surface on a 
! cell-by-cell basis. For each cell, a fourth-order polynomial of the form:
!```
! Z = Ax²y² + Bx²y + Cxy² + Dx² + Ey² + Fxy + Gx + Hy + I
!```
! is fit to a surface composed of a 3x3 window. 
!```
! A = [(Z1 + Z3 + Z7 + Z9) / 4  - (Z2 + Z4 + Z6 + Z8) / 2 + Z5] / L4
! B = [(Z1 + Z3 - Z7 - Z9) /4 - (Z2 - Z8) /2] / L3
! C = [(-Z1 + Z3 - Z7 + Z9) /4 + (Z4 - Z6)] /2] / L3
! D = [(Z4 + Z6) /2 - Z5] / L2
! E = [(Z2 + Z8) /2 - Z5] / L2
! F = (-Z1 + Z3 + Z7 - Z9) / 4L2
! G = (-Z4 + Z6) / 2L
! H = (Z2 - Z8) / 2L
! I = Z5
! 
!   curvature = -2 (E + D)
!
!   profileCurvature = -2 (D G^2 + E H^2 + F G H) / (G^2 + H^2)
!         
!   planCurvature  = -2 (D H^2 + E G^2 - F G H) / (G^2 + H^2)
!```
! Profile curvature is parallel to the direction of the maximum slope. 
! A negative value indicates that the surface is upwardly convex at 
! that cell. A positive profile indicates that the surface is upwardly 
! concave at that cell. A value of zero indicates that the surface 
! is linear. Planform curvature (commonly called plan curvature) is 
! perpendicular to the direction of the maximum slope. A positive value 
! indicates the surface is sidewardly convex at that cell. A negative plan
! indicates the surface is sidewardly concave at that cell. A value of 
! zero indicates the surface is linear. Plan curvature relates to the 
! convergence and divergence of flow across a surface.
!
! References:
!
! Moore, I. D., R. B. Grayson, and A. R. Landson. 1991. Digital Terrain 
!    Modelling: A Review of Hydrological, Geomorphological, and Biological 
!    Applications. Hydrological Processes 5: 3–30.
!
! Zeverbergen, L. W., and C. R. Thorne. 1987. Quantitative Analysis of 
!    Land Surface Topography. Earth Surface Processes and Landforms 12: 47–56.
SUBROUTINE DeriveCurvature &
!
(dem, curvature, profileCurvature, planCurvature)


IMPLICIT NONE

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

!arguments with intent out
TYPE (grid_real), INTENT (out) :: curvature
TYPE (grid_real), OPTIONAL, INTENT (out) :: profileCurvature
TYPE (grid_real), OPTIONAL, INTENT (out) :: planCurvature

!local variables
INTEGER :: i,j
REAL (KIND = float) :: D, E, F, G, H
REAL (KIND = float) :: z1, z2, z3, z4, z5, z6, z7, z8, z9
REAL (KIND = float) :: L

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

!allocate grids
CALL NewGrid (curvature, dem)

IF (PRESENT(profileCurvature)) THEN
  CALL NewGrid (profileCurvature, dem)
END IF

IF (PRESENT(planCurvature)) THEN
  CALL NewGrid (planCurvature, dem)
END IF


!loop through dem grid
idim: DO i = 1, dem % idim
  jdim: DO j = 1, dem % jdim
    IF (dem % mat(i,j) /= dem % nodata) THEN
    
       !set length
       L = CellArea (dem,i,j) ** 0.5
       
       !set z5
       z5 = dem % mat (i,j)
    
       !set z2
        IF ( .NOT. IsOutOfGrid (i-1,j,dem) ) THEN
          IF ( dem % mat (i-1,j) /= dem % nodata) THEN
             z2 = dem % mat (i-1,j)
          ELSE
             z2 = dem % mat (i,j)
          END IF  
        ELSE 
           z2 = dem % mat (i,j)
        END IF
        
        !set z3
        IF ( .NOT. IsOutOfGrid (i-1,j+1,dem) ) THEN
          IF ( dem % mat (i-1,j+1) /= dem % nodata) THEN
             z3 = dem % mat (i-1,j+1)
          ELSE
             z3 = dem % mat (i,j)
          END IF  
        ELSE 
           z3 = dem % mat (i,j)
        END IF
        
        !set z6
        IF ( .NOT. IsOutOfGrid (i,j+1,dem) ) THEN
          IF ( dem % mat (i,j+1) /= dem % nodata) THEN
             z6 = dem % mat (i,j+1)
          ELSE
             z6 = dem % mat (i,j)
          END IF  
        ELSE 
           z6 = dem % mat (i,j)
        END IF
        
        !set z9
        IF ( .NOT. IsOutOfGrid (i+1,j+1,dem) ) THEN
          IF ( dem % mat (i+1,j+1) /= dem % nodata) THEN
             z9 = dem % mat (i+1,j+1)
          ELSE
             z9 = dem % mat (i,j)
          END IF  
        ELSE 
           z9 = dem % mat (i,j)
        END IF
        
        !set z8
        IF ( .NOT. IsOutOfGrid (i+1,j,dem) ) THEN
          IF ( dem % mat (i+1,j) /= dem % nodata) THEN
             z8 = dem % mat (i+1,j)
          ELSE
             z8 = dem % mat (i,j)
          END IF  
        ELSE 
           z8 = dem % mat (i,j)
        END IF
        
        !set z7
        IF ( .NOT. IsOutOfGrid (i+1,j-1,dem) ) THEN
          IF ( dem % mat (i+1,j-1) /= dem % nodata) THEN
             z7 = dem % mat (i+1,j-1)
          ELSE
             z7 = dem % mat (i,j)
          END IF  
        ELSE 
           z7 = dem % mat (i,j)
        END IF
        
        !set z4
        IF ( .NOT. IsOutOfGrid (i,j-1,dem) ) THEN
          IF ( dem % mat (i,j-1) /= dem % nodata) THEN
             z4 = dem % mat (i,j-1)
          ELSE
             z4 = dem % mat (i,j)
          END IF  
        ELSE 
           z4 = dem % mat (i,j)
        END IF
        
         !set z1
        IF ( .NOT. IsOutOfGrid (i-1,j-1,dem) ) THEN
          IF ( dem % mat (i-1,j-1) /= dem % nodata) THEN
             z1 = dem % mat (i-1,j-1)
          ELSE
             z1 = dem % mat (i,j)
          END IF  
        ELSE 
           z1 = dem % mat (i,j)
        END IF
        
        !compute D, E, F, G, H
         D = ((z4 + z6) /2. - z5) / L**2
         E = ((z2 + z8) /2. - z5) / L**2
         F = (-z1 + z3 + z7 - z9) / (4. * L**2)
         G = (-z4 + z6) / (2. * L)
         H = (z2 - z8) / (2. * L)
         
         curvature % mat (i,j) = -2. * (E + D)
         
         IF (PRESENT(profileCurvature)) THEN
           IF ( (G**2. + H**2.) > 0.) THEN
               profileCurvature % mat (i,j) = -2. * (D*G**2. + E*H**2. + F*G*H) / (G**2. + H**2.)
           ELSE
               profileCurvature % mat (i,j) = 0.
           END IF
         END IF
         
         IF (PRESENT(planCurvature)) THEN
            IF ( (G**2. + H**2.) > 0.) THEN
               planCurvature % mat (i,j) = -2. * (D*H**2. + E*G**2. - F*G*H) / (G**2. + H**2.)
            ELSE
                planCurvature % mat (i,j) = 0.
            END IF 
         END IF
       
    
    END IF
  END DO jdim
END DO idim


RETURN

END SUBROUTINE DeriveCurvature



!==============================================================================
!| Description:
!   compute map of flow accumulation (m2)
!   Input grid: flow direction
SUBROUTINE FlowAccumulation &
!
(fdir, facc)

IMPLICIT NONE

!arguments with intent in
TYPE(grid_integer), INTENT(in):: fdir

!arguments with intent out
TYPE (grid_real), INTENT (out) :: facc

!local variables
INTEGER :: i,j

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

!allocate new flow accumulation grid using flow direction grid as template
CALL NewGrid (facc, fdir)

!compute grid containing area of each cell
DO i = 1, fdir % idim
  DO j = 1, fdir % jdim
    IF (fdir % mat (i,j) /= fdir % nodata) THEN
      facc % mat(i,j) = CellArea (facc, i, j)
    ELSE
      facc % mat(i,j) = facc % nodata
    END IF
   END DO
END DO

DO i = 1, fdir % idim
  DO j = 1, fdir % jdim
    IF (fdir % mat (i,j) /= fdir % nodata) THEN
      CALL BasinArea (i, j, fdir, facc % mat(i,j))
    ELSE
      facc % mat(i,j) = facc % nodata
    END IF
  END DO
END DO


END SUBROUTINE FlowAccumulation



!==============================================================================
!| Description:
!   compute basin area (m2)
RECURSIVE SUBROUTINE BasinArea &
!
(r, c, fdir, area)

IMPLICIT NONE

TYPE(grid_integer),INTENT(IN) :: fdir
INTEGER, INTENt(in) :: r,c 
REAL (KIND = float), INTENT(inout) :: area
!------------------------------end of declaration -----------------------------

IF ( .NOT. IsOutOfGrid(r,c+1,fdir) ) THEN
    IF(fdir%mat(r,c+1) == W) THEN
       area = area + CellArea (fdir, r, c+1)
       CALL BasinArea (r,c+1,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r+1,c+1,fdir) ) THEN
    IF(fdir%mat(r+1,c+1) == NW ) THEN
       area = area + CellArea (fdir, r+1, c+1)
       CALL BasinArea (r+1,c+1,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r+1,c,fdir) ) THEN
    IF(fdir%mat(r+1,c) == N ) THEN
       area = area + CellArea (fdir, r+1, c)
       CALL BasinArea (r+1,c,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r+1,c-1,fdir) ) THEN
    IF(fdir%mat(r+1,c-1) == NE ) THEN
       area = area + CellArea (fdir, r+1, c-1)
       CALL BasinArea (r+1,c-1,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r,c-1,fdir) ) THEN
    IF(fdir%mat(r,c-1) == E) THEN
        area = area + CellArea (fdir, r, c-1)
        CALL BasinArea (r,c-1,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r-1,c-1,fdir) ) THEN
    IF(fdir%mat(r-1,c-1) == SE ) THEN
       area = area + CellArea (fdir, r-1, c-1)
       CALL BasinArea (r-1,c-1,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r-1,c,fdir) ) THEN
    IF(fdir%mat(r-1,c) == S ) THEN
       area = area + CellArea (fdir, r-1, c)
       CALL BasinArea (r-1,c,fdir,area)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r-1,c+1,fdir) ) THEN
    IF(fdir%mat(r-1,c+1) == SW ) THEN
       area = area + CellArea (fdir, r-1, c+1)
       CALL BasinArea (r-1,c+1,fdir,area) 
    END IF
END IF

END SUBROUTINE BasinArea



!==============================================================================
!| Description:
!   find the cells that are springs, defined as those cells that have not
!   any other cells upstream
FUNCTION CellIsSpring &
!
(row, col, flowDir) &
!
RESULT (spring)

IMPLICIT NONE

! Arguments with intent(in):
INTEGER, INTENT(in) :: row, col
TYPE (grid_integer), INTENT(in) :: flowDir

! Local variables:
LOGICAL :: spring

!------------end of declaration------------------------------------------------
spring = .TRUE.

IF (flowDir % mat(row,col) == flowDir % nodata) THEN
  spring = .FALSE.
  RETURN
END IF

IF(.NOT. IsOutOfGrid(row,col+1,flowDir) ) THEN
    IF(flowDir%mat(row,col+1) == W ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row+1,col+1,flowDir) ) THEN
    IF(flowDir%mat(row+1,col+1) == NW  ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row+1,col,flowDir) ) THEN
    IF(flowDir%mat(row+1,col) == N  ) THEN
      spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row+1,col-1,flowDir) ) THEN
    IF(flowDir%mat(row+1,col-1) == NE  ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row,col-1,flowDir) ) THEN
    IF(flowDir%mat(row,col-1) == E  ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row-1,col-1,flowDir) ) THEN
    IF(flowDir%mat(row-1,col-1) == SE  ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row-1,col,flowDir) ) THEN
    IF(flowDir%mat(row-1,col) == S  ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF

IF(.NOT. IsOutOfGrid(row-1,col+1,flowDir) ) THEN
    IF(flowDir%mat(row-1,col+1) == SW  ) THEN
       spring = .FALSE.
       RETURN
    ENDIF
ENDIF


END FUNCTION CellIsSpring


!==============================================================================
!| Description:
!   compute mask of river basin given map of flow direction and the coordinate
!   of the outlet point
SUBROUTINE BasinDelineate &
!
(fdir,x,y, mask)

IMPLICIT NONE

!arguments with intent in
TYPE(grid_integer), INTENT(in):: fdir
REAL (KIND = float), INTENT(in) :: x, y !!coordinate of outlet

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

!local variables
INTEGER :: i,j
LOGICAL :: check

!------------------------------end of declaration -----------------------------
                                         
IF (.NOT. ALLOCATED (mask % mat)) THEN
  !allocate new grid
  CALL NewGrid (mask, fdir)
END IF

mask % mat = mask % nodata

!compute row and column of outlet
CALL GetIJ (x, y, fdir, i, j, check)
	
!compute basin mask
CALL BasinMask(mask,fdir,i,j)

END SUBROUTINE BasinDelineate
	


!==============================================================================
!| Description:
!   search for cells included in the river basin
RECURSIVE SUBROUTINE BasinMask &
!
(basin, fdir, r, c)

IMPLICIT NONE

TYPE(grid_integer),INTENT(IN) :: fdir
TYPE(grid_integer),INTENT(INOUT) :: basin
INTEGER, INTENT(in) :: r,c  

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

IF ( .NOT. IsOutOfGrid(r,c+1,fdir) ) THEN
    IF(fdir%mat(r,c+1)==W.AND. basin%mat(r,c+1)/=1) THEN
       basin%mat(r,c+1) = 1
       CALL BasinMask(basin,fdir,r,c+1) 
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r+1,c+1,fdir) ) THEN
    IF(fdir%mat(r+1,c+1)==NW .AND. basin%mat(r+1,c+1)/=1) THEN
       basin%mat(r+1,c+1) = 1
       CALL BasinMask(basin,fdir,r+1,c+1) 
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r+1,c,fdir) ) THEN
    IF(fdir%mat(r+1,c)==N .AND. basin%mat(r+1,c)/=1) THEN
       basin%mat(r+1,c) = 1
       CALL BasinMask(basin,fdir,r+1,c) 
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r+1,c-1,fdir) ) THEN
    IF(fdir%mat(r+1,c-1)==NE .AND. basin%mat(r+1,c-1)/=1) THEN
       basin%mat(r+1,c-1) = 1
       CALL BasinMask(basin,fdir,r+1,c-1)
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r,c-1,fdir) ) THEN
    IF(fdir%mat(r,c-1)==E .AND. basin%mat(r,c-1)/=1) THEN
       basin%mat(r,c-1) = 1
       CALL BasinMask(basin,fdir,r,c-1) 
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r-1,c-1,fdir) ) THEN
    IF(fdir%mat(r-1,c-1)==SE .AND. basin%mat(r-1,c-1)/=1) THEN
       basin%mat(r-1,c-1) = 1
       CALL BasinMask(basin,fdir,r-1,c-1) 
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r-1,c,fdir) ) THEN
    IF(fdir%mat(r-1,c)==S .AND. basin%mat(r-1,c)/=1) THEN
       basin%mat(r-1,c) = 1
       CALL BasinMask(basin,fdir,r-1,c) 
    END IF
END IF

IF ( .NOT. IsOutOfGrid(r-1,c+1,fdir) ) THEN
    IF(fdir%mat(r-1,c+1)==SW .AND. basin%mat(r-1,c+1)/=1) THEN
       basin%mat(r-1,c+1) = 1
       CALL BasinMask(basin,fdir,r-1,c+1) 
    END IF
END IF

END SUBROUTINE BasinMask


!==============================================================================
!| Description:
!   compute longest flow length (m)
FUNCTION LongestFlowLength &
!
(fdir, x, y) &
!
RESULT (lmax)

IMPLICIT NONE

!Arguments with intent (in)
TYPE(grid_integer),INTENT(IN) :: fdir
REAL (KIND = float), INTENT(in) :: x, y  

!local declarations
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) :: i, j
LOGICAL :: outlet

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

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

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


point1 % system = basin % grid_mapping
point2 % system = basin % grid_mapping
lmax = 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.
       
              !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
                    lmax = length
                 END IF
              END IF
              
              !loop
              row = iDown
              col = jDown

           END DO
                      
        ENDIF
    END IF
  ENDDO
ENDDO 


RETURN
END FUNCTION LongestFlowLength


!==============================================================================
!| Description:
!   compute slope of longest flow path (m/m)
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




!==============================================================================
!| Description:
!   Drainage Density (1/m) of a basin is the total line length of all streams 
!   divided by basin area. If flow accumulation is not given, it is computed
!   from flow direction. If coordinate of closing section is not given
!   drainage density is computed for the entire grid.
FUNCTION DrainageDensity &
!
(channel, fdir, x, y, flowAcc) &
!
RESULT (dd)

IMPLICIT NONE

!Arguments with intent (in)
TYPE(grid_integer),INTENT(IN) :: channel !!mask of channel cells. NoData for non channel cells
TYPE(grid_integer),INTENT(IN) :: fdir !!flow direction
REAL (KIND = float), INTENT(IN),OPTIONAL :: x, y !!coordinate of basin closing section 
TYPE (grid_real),INTENT(IN),OPTIONAL :: flowAcc !!flow accumulation

!local declarations
REAL (KIND = float) :: dd
TYPE (grid_real) :: facc
TYPE (grid_integer) :: mask
INTEGER (KIND = short) :: i,j
INTEGER (KIND = short) :: is, js !!coordinate of downstream cell in local reference system
REAL (KIND = float) :: length !!total river length
REAL (KIND = float) :: cLength !!cell length
LOGICAL :: check

!------------------------------end of declaration -----------------------------
!set flowaccumulation
IF (PRESENT (flowAcc)) THEN
   CALL NewGrid (facc, fdir)
   facc = flowAcc

ELSE
  CALL FlowAccumulation (fdir, facc)
END IF

IF (PRESENT(x) .AND. PRESENT(y)) THEN !delineate watershed mask
   CALL BasinDelineate (fdir,x,y, mask)
ELSE  !drainage density is computed on entire grid
   CALL NewGrid (mask, fdir)
END IF

!compute total channel length [m]
length = 0.
DO i = 1, mask % idim
  DO j = 1, mask % jdim
    IF (mask % mat (i,j) /= mask % nodata) THEN
       IF (channel % mat (i,j) /= channel % nodata) THEN
          CALL DownstreamCell (i, j, fdir % mat(i,j), is, js, dx = cLength, grid = fdir)
          length = length + cLength
       END IF
    END IF
  END DO
END DO

!compute drainage density [1/m]
CALL GetIJ (x, y, fdir, i, j, check)

IF (.NOT. check) THEN
  CALL Catch ('error', 'Morphology', 'Drainage density computation: ', &
             argument = 'closing section outside grid dimension')
END IF

IF (channel % mat(i,j) == channel % nodata) THEN
  CALL Catch ('error', 'Morphology', 'Drainage density computation: ', &
             argument = 'closing section on hillslope')
END IF

dd = length / facc % mat(i,j) 

!deallocate
CALL GridDestroy (facc)
CALL GridDestroy (mask)

RETURN

END FUNCTION DrainageDensity


!==============================================================================
!| Description:
!   Topographic index `TI` takes into account both a local slope geometry 
!   and site location in the landscape combining data on slope steepness `S` 
!   and specific catchment area `As`:
!```  
!   TI = ln (As/S)
!```  
!  For grid structures the specific catchment area  can be obtained by 
!  dividing the contributing area of the cell by the effective contour length. 
!  This is the length of contour line within the grid cell over which flow 
!  can pass. It is calculated as the length of the line through the grid 
!  cell centre and perpendicular to the aspect direction.
!
!  References:
!
!  P. J. J. DESMET & G. GOVERS (1996): Comparison of routing algorithms 
!     for digital elevation models and their implications for predicting 
!     ephemeral gullies, International Journal of Geographical 
!     Information Systems, 10:3, 311-331
SUBROUTINE TopographicIndex &
!
(dem, fdir, facc, xcoord, ycoord, tivalue, tigrid)

IMPLICIT NONE

!Argumentd with intent(in):
TYPE (grid_real), INTENT(IN) :: dem !!digital elevation model
TYPE (grid_integer), INTENT(IN) :: fdir !!flow direction
TYPE (grid_real), OPTIONAL, INTENT(IN) :: facc !!flow accumulation in cell unit
REAL (KIND = float) , OPTIONAL, INTENT(IN) :: xcoord !! x coordinate of cell for computing index
REAL (KIND = float) , OPTIONAL, INTENT(IN) :: ycoord !! y coordinate of cell for computing index

!Arguments with intent(out):
REAL (KIND = float) , OPTIONAL, INTENT(OUT) :: tivalue !! index of cell with coordinates (x,y)
TYPE (grid_real), OPTIONAL, INTENT(OUT) :: tigrid !!topographic index

!local declarations
TYPE (grid_real) :: area !! contributing area (m^2)
TYPE (grid_real) :: aspect !!aspect (rad)
TYPE (grid_real) :: slope !!slope (rad)
REAL (KIND = float) :: as !! specific catchment area (m^2 / m)
REAL (KIND = float) :: x !!factor converting cellsize along direction...
                         !!...perpendicular to the aspect direction
REAL (KIND = float) :: alpha !! local aspect (rad)
REAL (KIND = float) :: cellsize !!cell size (m)
INTEGER (KIND = short) :: i, j
!------------------------------end of declaration -----------------------------

!allocate topographic index grid
CALL NewGrid (tigrid, dem)

IF ( .NOT. PRESENT(facc) ) THEN
   !compute contributing area
   CALL FlowAccumulation (fdir, area)
ELSE
    CALL  NewGrid (area, fdir)
    DO i = 1, area % idim
      DO j = 1, area % jdim
          IF (area % mat (i,j) /= area % nodata ) THEN
              area % mat (i,j) = facc % mat (i,j) ! (* CellArea (facc, i, j))
          END IF
      END DO
    END DO
END IF

!compute aspect
CALL DeriveAspect (dem, aspect)

!compute slope
CALL DeriveSlope (dem, slope)

!if present tivalue compute topographic index only for cell with coordinate xcoord, ycoord
IF (PRESENT (tivalue) ) THEN
    !check wether coordinates of cell is given
    IF ( PRESENT (xcoord) .AND. PRESENT (ycoord) ) THEN
        CALL GetIJ (xcoord, ycoord, fdir, i, j)
        
        alpha = aspect % mat (i,j)
        x = 1. / MAX(ABS(SIN(alpha)),ABS(COS(alpha)))
        cellsize = ( CellArea (fdir, i, j) ) ** 0.5 ![m]
        as = area % mat (i,j) / (cellsize * x)
        IF ( TAN (slope % mat (i,j)) == 0.) THEN
            tivalue = LOG (as / 0.0001 )
        ELSE
            tivalue = LOG (as / TAN (slope % mat (i,j)))
        END IF
 
    ELSE
       CALL Catch ('error', 'Morphology', 'missing coordinate for compuing topographic index of given cell')
    END IF
    
END IF

    
!If present tigrid, compute topographic index for every cell of the grid
IF (PRESENT (tigrid) ) THEN

        DO i = 1, tigrid % idim
          DO j = 1, tigrid %  jdim
            IF (tigrid % mat (i,j) /= tigrid % nodata) THEN
              alpha = aspect % mat (i,j)
              x = 1. / MAX(ABS(SIN(alpha)),ABS(COS(alpha)))
              cellsize = ( CellArea (tigrid, i, j) ) ** 0.5 ![m]
              as = area % mat (i,j) / (cellsize * x)
              IF ( TAN (slope % mat (i,j)) == 0.) THEN
                 tigrid % mat (i,j) = LOG (as / 0.0001)
              ELSE
                 tigrid % mat (i,j) = LOG (as / TAN (slope % mat (i,j)))
              END IF
     
            END IF 
          END DO
        END DO
        
END IF


!purge
CALL GridDestroy (area)
CALL GridDestroy (aspect)
CALL GridDestroy (slope)

RETURN

END SUBROUTINE TopographicIndex


!==============================================================================
!| Description:
!   find the cells that are springs, defined as those cells that have not
!   any other cells upstream
FUNCTION Perimeter &
!
(grid) &
!
RESULT (p)

IMPLICIT NONE

! Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid

! Local variables:
REAL (KIND = float) :: p
TYPE (grid_integer) :: border
INTEGER :: i,j

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


CALL ExtractBorder (grid, border, cardinal = .TRUE.)

p = 0.

DO i = 1, border % idim
    DO j = 1, border % jdim
        IF (border % mat (i,j) /= border % nodata ) THEN
            p = p + CellLength (border, i,j)
        END IF
    END DO
END DO

!increase a default value of 9% to consider diagonal length

p = p * 1.09

CALL GridDestroy (border)


RETURN
END FUNCTION Perimeter



END MODULE Morphology