!! compute grid statistics
!|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 - 24th Jan 2020  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 19/Dec/2016 | routines moved from GridOperations module |
! | 1.1      | 04/Sep/2018 | GetCV to compute coefficient of variation of a grid |
! | 1.2      | 24/Jan/2020 | Added subroutine UniqueValues to retrieve a list of distinct values from a grid_integer |
!
!
!### 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
!
!### Module Description: 
!   library to compute grids statistics
MODULE GridStatistics         
			

! Modules used: 

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

USE LogLib, ONLY: &
  ! Imported routines:
  Catch
 
USE GridLib, ONLY: &
  !Imported type definitions:
  grid_integer, grid_real

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

!USE GeoLib, ONLY: &
!Imported operators:
!OPERATOR (==)


IMPLICIT NONE 

!global (public declarations)

INTERFACE GetMean
   MODULE PROCEDURE GetMeanOfGridFloat
   MODULE PROCEDURE GetMeanOfGridInteger
END INTERFACE

INTERFACE GetStDev
   MODULE PROCEDURE GetStDevOfGridFloat
   MODULE PROCEDURE GetStDevOfGridInteger
END INTERFACE

INTERFACE GetMin
   MODULE PROCEDURE GetMinOfGridFloat
   MODULE PROCEDURE GetMinOfGridInteger
END INTERFACE

INTERFACE GetCV
   MODULE PROCEDURE GetCVofGridFloat
   !MODULE PROCEDURE GetCVofGridInteger
END INTERFACE

INTERFACE GetMax
   MODULE PROCEDURE GetMaxOfGridFloat
   MODULE PROCEDURE GetMaxOfGridInteger
END INTERFACE

INTERFACE GetSum
   MODULE PROCEDURE GetSumOfGridFloat
   MODULE PROCEDURE GetSumOfGridInteger
END INTERFACE

INTERFACE CountCells
   MODULE PROCEDURE CountCellsFloat
   MODULE PROCEDURE CountCellsInteger
END INTERFACE

INTERFACE MinMaxNormalization
   MODULE PROCEDURE MinMaxNormalizationFloat
   !MODULE PROCEDURE MinMaxNormalizationInteger !seems not useful
END INTERFACE

INTERFACE GetBias
   MODULE PROCEDURE GetBiasFloat
   !MODULE PROCEDURE GetBiasInteger !TO BE IMPLEMENTED IF REQUIRED
END INTERFACE

INTERFACE GetRMSE
   MODULE PROCEDURE GetRMSEfloat
   !MODULE PROCEDURE GetRMSEinteger !TO BE IMPLEMENTED IF REQUIRED
END INTERFACE

INTERFACE GetNSE
   MODULE PROCEDURE GetNSEfloat
   !MODULE PROCEDURE GetNSEinteger !TO BE IMPLEMENTED IF REQUIRED
END INTERFACE

INTERFACE GetR2
   MODULE PROCEDURE GetR2float
   !MODULE PROCEDURE GetR2integer !TO BE IMPLEMENTED IF REQUIRED
END INTERFACE

INTERFACE UniqueValues
   MODULE PROCEDURE UniqueValuesOfGridInteger
END INTERFACE

!global procedures:
PUBLIC :: GetMean
PUBLIC :: GetStDev
PUBLIC :: GetMin
PUBLIC :: GetMax
PUBLIC :: GetSum
PUBLIC :: CountCells
PUBLIC :: MinMaxNormalization
PUBLIC :: GetRMSE
PUBLIC :: GetBias
PUBLIC :: GetNSE
PUBLIC :: GetR2
PUBLIC :: GetCV
PUBLIC :: UniqueValues


! Local (i.e. private) Declarations: 
! Local Procedures:

PRIVATE :: GetMeanOfGridFloat
PRIVATE :: GetMeanOfGridInteger
PRIVATE :: GetStDevOfGridFloat
PRIVATE :: GetStDevOfGridInteger
PRIVATE :: GetMinOfGridFloat
PRIVATE :: GetMinOfGridInteger
PRIVATE :: GetMaxOfGridFloat
PRIVATE :: GetMaxOfGridInteger
PRIVATE :: GetSumOfGridFloat
PRIVATE :: GetSumOfGridInteger
PRIVATE :: CountCellsFloat
PRIVATE :: CountCellsInteger
PRIVATE :: MinMaxNormalizationFloat
PRIVATE :: GetBiasFloat
PRIVATE :: GetRMSEfloat
PRIVATE :: GetR2float

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





!==============================================================================
!| Description:
!   compute mean of grid_real eventually constrained to a mask
FUNCTION GetMeanOfGridFloat &
!
(grid, maskReal, maskInteger) &
!
RESULT (mean)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: mean
REAL (KIND = float) :: countCells
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
mean = 0.
countCells = 0.
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                IF (grid % mat (i,j) /= grid % nodata) THEN
                   CALL Catch ('warning', 'GridStatistics',  &
                             'calculate mean: ', argument = &
                              'found nodata within mask' ) 
                   mean = mean + grid % mat (i,j)
                   countCells = countCells + 1.
                END IF
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 IF (grid % mat (i,j) /= grid % nodata) THEN
                     CALL Catch ('warning', 'GridStatistics',  &
                             'calculate mean: ', argument = &
                              'found nodata within mask' ) 
                     mean = mean + grid % mat (i,j)
                     countCells = countCells + 1
                 END IF
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                mean = mean + grid % mat (i,j)
                countCells = countCells + 1.
            END IF
        END DO
    END DO
END IF

mean = mean / countCells

RETURN

END FUNCTION GetMeanOfGridFloat


!==============================================================================
!| Description:
!   compute mean of grid_integer eventually constrained to a mask
FUNCTION GetMeanOfGridInteger &
!
(grid, maskReal, maskInteger) &
!
RESULT (mean)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger

!Local declarations:
REAL (KIND = float) :: mean
REAL (KIND = float) :: countCells
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
mean = 0.
countCells = 0.
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                mean = mean + grid % mat (i,j)
                countCells = countCells + 1.
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                mean = mean + grid % mat (i,j)
                countCells = countCells + 1.
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                mean = mean + grid % mat (i,j)
                countCells = countCells + 1.
            END IF
        END DO
    END DO
END IF

mean = mean / countCells

RETURN

END FUNCTION GetMeanOfGridInteger


!==============================================================================
!| Description:
!   compute unbiased standard deviation of grid_real optionally 
!   constrained to a mask
FUNCTION GetStDevOfGridFloat &
!
(grid, maskReal, maskInteger) &
!
RESULT (stDev)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: mean
REAL (KIND = float) :: stDev
REAL (KIND = float) :: countCells
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
stDev = 0.

!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate standard deviation: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean
    mean = GetMean (grid, maskReal = maskReal)

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                stDev = stDev + ( grid % mat (i,j) - mean) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate standard deviation: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean
    mean = GetMean (grid, maskInteger = maskInteger)

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                stDev = stDev + ( grid % mat (i,j) - mean) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO

ELSE

    !get mean
    mean = GetMean (grid)
    
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                stDev = stDev + ( grid % mat (i,j) - mean) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
END IF

stDev = (stDev / (countCells - 1.)) ** 0.5

RETURN

END FUNCTION GetStDevOfGridFloat


!==============================================================================
!| Description:
!   compute unbiased standard deviation of grid_integer optionally 
!   constrained to a mask
FUNCTION GetStDevOfGridInteger &
!
(grid, maskReal, maskInteger) &
!
RESULT (stDev)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: mean
REAL (KIND = float) :: stDev
REAL (KIND = float) :: countCells
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
stDev = 0.

!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate standard deviation: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean
    mean = GetMean (grid, maskReal = maskReal)

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                stDev = stDev + ( grid % mat (i,j) - mean) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate standard deviation: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean
    mean = GetMean (grid, maskInteger = maskInteger)

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                stDev = stDev + ( grid % mat (i,j) - mean) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO

ELSE

    !get mean
    mean = GetMean (grid)
    
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                stDev = stDev + ( grid % mat (i,j) - mean) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
END IF

stDev = (stDev / (countCells - 1.)) ** 0.5

RETURN

END FUNCTION GetStDevOfGridInteger



!==============================================================================
!| Description:
!   compute coefficient of variation of grid_real eventually constrained to a mask
FUNCTION GetCVofGridFloat &
!
(grid, maskReal, maskInteger) &
!
RESULT (cv)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: cv
REAL (KIND = float) :: mean
REAL (KIND = float) :: std

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

!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate coefficient of variation: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    mean = GetMean (grid, maskReal = maskReal)
    std = GetStDev (grid, maskReal = maskReal)

ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate coefficient of variation: ', argument = &
        'coordinate reference system of mask differs from input grid' )
        
    END IF
    
    mean = GetMean (grid, maskInteger = maskInteger)
    std = GetStDev (grid, maskInteger = maskInteger)

ELSE
    
    mean = GetMean (grid)
    std = GetStDev (grid)
    
END IF

cv = std / mean

RETURN

END FUNCTION GetCVofGridFloat



!==============================================================================
!| Description:
!   compute minimum of grid_real eventually constrained to a mask
FUNCTION GetMinOfGridFloat &
!
(grid, maskReal, maskInteger) &
!
RESULT (min)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: min
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
min = HUGE (min)
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate min: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                IF (grid % mat(i,j) < min) THEN
                   min = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate min: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                IF (grid % mat(i,j) < min) THEN
                   min = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                IF (grid % mat(i,j) < min) THEN
                   min = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
END IF

RETURN

END FUNCTION GetMinOfGridFloat


!==============================================================================
!| Description:
!   compute minimum of grid_integer eventually constrained to a mask
FUNCTION GetMinOfGridInteger &
!
(grid, maskReal, maskInteger) &
!
RESULT (min)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
INTEGER (KIND = long) :: min
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
min = HUGE (min)
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate min: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                IF (grid % mat(i,j) < min) THEN
                   min = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate min: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                IF (grid % mat(i,j) < min) THEN
                   min = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                IF (grid % mat(i,j) < min) THEN
                   min = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
END IF

RETURN

END FUNCTION GetMinOfGridInteger


!==============================================================================
!| Description:
!   compute maximum of grid_real eventually constrained to a mask
FUNCTION GetMaxOfGridFloat &
!
(grid, maskReal, maskInteger) &
!
RESULT (max)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: max
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
max = - HUGE (max)
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate max: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                IF (grid % mat(i,j) > max) THEN
                   max = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate max: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                IF (grid % mat(i,j) > max) THEN
                   max = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                IF (grid % mat(i,j) > max) THEN
                   max = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
END IF

RETURN

END FUNCTION GetMaxOfGridFloat


!==============================================================================
!| Description:
!   compute maximum of grid_integer eventually constrained to a mask
FUNCTION GetMaxOfGridInteger &
!
(grid, maskReal, maskInteger) &
!
RESULT (max)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
INTEGER (KIND = long) :: max
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
max = - HUGE (max)
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate max: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                IF (grid % mat(i,j) > max) THEN
                   max = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate max: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                IF (grid % mat(i,j) > max) THEN
                   max = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                IF (grid % mat(i,j) > max) THEN
                   max = grid % mat(i,j)
                END IF
            END IF
        END DO
    END DO
END IF

RETURN

END FUNCTION GetMaxOfGridInteger


!==============================================================================
!| Description:
!   compute Root Mean Square Error between two grids real 
!   optional argument:
!     `mask`: compute RMSE only on assigned mask
!     `nrmse`: compute Normalizes Root Mean Square Error
FUNCTION GetRMSEfloat &
!
(grid1, grid2, maskReal, maskInteger, nrmse) &
!
RESULT (rmse)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid1
TYPE (grid_real), INTENT(IN) :: grid2
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger
LOGICAL, OPTIONAL, INTENT(IN) :: nrmse


!Local declarations:
INTEGER (KIND = long) :: i, j
REAL (KIND = float) :: rmse, mean
INTEGER             :: countCells
!---------------------------end of declarations--------------------------------

rmse = 0.
countCells = 0

!check grid1 and grid2 have the same coordinate reference system
 IF ( .NOT. CRSisEqual(grid1,grid2) ) THEN
      CALL Catch ('error', 'GridStatistics',  &
      'calculate RMSE: ', argument = &
      'coordinate reference system of grid1 differs from grid2' )
END IF

IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate RMSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean
    !mean = GetMean (grid, maskReal = maskReal)

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                rmse = rmse + ( grid1 % mat (i,j) - grid2 % mat (i,j)) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
    
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate RMSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
  
    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 rmse = rmse + ( grid1 % mat (i,j) - grid2 % mat (i,j)) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO

ELSE
    
    DO j = 1, grid1 % jdim
        DO i = 1, grid1 % idim
            IF (grid1 % mat(i,j) /= grid1 % nodata) THEN
                rmse = rmse + ( grid1 % mat (i,j) - grid2 % mat (i,j)) **2.
                countCells = countCells + 1.
            END IF
        END DO
    END DO
END IF

rmse = ( rmse / countCells ) ** 0.5

IF (PRESENT(nrmse)) THEN
  IF (nrmse) THEN
    mean = GetMean (grid1)
    rmse = rmse / mean
  END IF
END IF

RETURN

END FUNCTION GetRMSEfloat


!==============================================================================
!| Description:
!   compute Nash and Sutcliffe efficiency index, equivalent to 
!   coefficient of determination.
!   optional argument:
!     `mask`: compute RMSE only on assigned mask
!     `nrmse`: compute normalized Nash efficiency
! ***
!  References: <br>
!    Nash, J. E. and J. V. Sutcliffe (1970), River flow forecasting through 
!      conceptual models part I — A discussion of principles, Journal of Hydrology, 
!      10 (3), 282–290. <br><br>
!    Moriasi, D. N.; Arnold, J. G.; Van Liew, M. W.; Bingner,R. L.; Harmel, R. D.; 
!      Veith, T. L. (2007), "Model Evaluation Guidelines for Systematic Quantification
!      of Accuracy in Watershed Simulations", Transactions of the ASABE, 50 (3), 885–900.
!      http://swat.tamu.edu/media/1312/moriasimodeleval.pdf <br><br>
!    Nossent, J., Bauwens, W., Application of a normalized Nash-Sutcliffe efficiency 
!      to improve the accuracy of the Sobol’ sensitivity analysis of a hydrological model
!      Geophysical Research Abstracts Vol. 14, EGU2012-237, 2012 EGU General Assembly 2012
FUNCTION GetNSEfloat &
!
(grid1, grid2, maskReal, maskInteger, normalized) &
!
RESULT (nash)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid1  !modeled
TYPE (grid_real), INTENT(IN) :: grid2  !observation
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger
LOGICAL, OPTIONAL, INTENT(IN) :: normalized


!Local declarations:
INTEGER (KIND = long) :: i, j
REAL (KIND = float) :: nash, meanobs, nashnum, nashden
!---------------------------end of declarations--------------------------------

nash = 0.
nashnum = 0.
nashden = 0.

!check grid1 and grid2 have the same coordinate reference system
 IF ( .NOT. CRSisEqual(grid1,grid2) ) THEN
      CALL Catch ('error', 'GridStatistics',  &
      'calculate NSE: ', argument = &
      'coordinate reference system of grid1 differs from grid2' )
END IF

IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate NSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean of observations
     meanobs = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                nashnum = nash + ( grid1 % mat (i,j) - grid2 % mat (i,j) ) **2.
                nashden = nashden + ( grid2 % mat (i,j) - meanobs ) **2.
            END IF
        END DO
    END DO
    
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate RMSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    !get mean of observations
     meanobs = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 nashnum = nashnum + ( grid1 % mat (i,j) - grid2 % mat (i,j) ) **2.
                 nashden = nashden + ( grid2 % mat (i,j) - meanobs ) **2.
            END IF
        END DO
    END DO

ELSE
   
    !get mean of observations
     meanobs = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, grid1 % jdim
        DO i = 1, grid1 % idim
            IF (grid1 % mat(i,j) /= grid1 % nodata) THEN
                 nashnum = nashnum + ( grid1 % mat (i,j) - grid2 % mat (i,j) ) **2.
                 nashden = nashden + ( grid2 % mat (i,j) - meanobs ) **2.
            END IF
        END DO
    END DO
END IF

IF (nashden /= 0.) THEN
   nash = 1. - nashnum / nashden
ELSE
   nash = -9999.9
END IF

IF (PRESENT(normalized)) THEN
  IF (normalized) THEN
    IF (nash == -9999.9) THEN
      nash = -9999.9
    ELSE
      nash = 1. / (2. - nash)
    END IF
  END IF
END IF

RETURN

END FUNCTION GetNSEfloat



!==============================================================================
!| Description:
!   compute R2, square of the Pearson correlation coefficient 
!   optional argument:
!     `mask`: compute RMSE only on assigned mask
! ***
!  References: <br>
!    Karl Pearson (20 June 1895) "Notes on regression and inheritance in 
!      the case of two parents," Proceedings of the Royal Society of 
!      London, 58 : 240–242
FUNCTION GetR2float &
!
(grid1, grid2, maskReal, maskInteger) &
!
RESULT (r2)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid1  
TYPE (grid_real), INTENT(IN) :: grid2 
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
INTEGER (KIND = long) :: i, j
REAL (KIND = float) :: r2, mean1, mean2, num, den1, den2, den
!---------------------------end of declarations--------------------------------

num = 0.
den = 0.
den1 = 0.
den2 = 0.

!check grid1 and grid2 have the same coordinate reference system
 IF ( .NOT. CRSisEqual(grid1,grid2) ) THEN
      CALL Catch ('error', 'GridStatistics',  &
      'calculate R2: ', argument = &
      'coordinate reference system of grid1 differs from grid2' )
END IF

IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate R2: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
    
    !get mean of grid1 and grid2
     mean1 = GetMean (grid1, maskReal = maskReal)
     mean2 = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                num = num + ( grid1 % mat (i,j) - mean1 ) * &
                            ( grid2 % mat (i,j) - mean2 )
                den1 = den1 + ( grid1 % mat (i,j) - mean1 ) ** 2.
                den2 = den2 + ( grid2 % mat (i,j) - mean2 ) ** 2.
            END IF
        END DO
    END DO
    
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate R2: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

     !get mean of grid1 and grid2
      mean1 = GetMean (grid1, maskReal = maskReal)
      mean2 = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 num = num + ( grid1 % mat (i,j) - mean1 ) * &
                            ( grid2 % mat (i,j) - mean2 )
                 den1 = den1 + ( grid1 % mat (i,j) - mean1 ) ** 2.
                 den2 = den2 + ( grid2 % mat (i,j) - mean2 ) ** 2.
            END IF
        END DO
    END DO

ELSE !compute over entire grid
   
     !get mean of grid1 and grid2
      mean1 = GetMean (grid1, maskReal = maskReal)
      mean2 = GetMean (grid2, maskReal = maskReal)
    
    !compute numerator and denominator
    DO j = 1, grid1 % jdim
        DO i = 1, grid1 % idim
            IF (grid1 % mat(i,j) /= grid1 % nodata) THEN
                 num = num + ( grid1 % mat (i,j) - mean1 ) * &
                            ( grid2 % mat (i,j) - mean2 )
                 den1 = den1 + ( grid1 % mat (i,j) - mean1 ) ** 2.
                 den2 = den2 + ( grid2 % mat (i,j) - mean2 ) ** 2.
            END IF
        END DO
    END DO
END IF

den = (den1 * den2) ** 0.5
IF (den /= 0.) THEN
   r2 = ( num / den ) ** 2.
ELSE
   r2 = -9999.9
END IF


RETURN

END FUNCTION GetR2float


!==============================================================================
!| Description:
!   compute mean Bias between two grids real 
!   optional argument:
!     `mask`: compute Bias only on assigned mask
!     `rbias`: compute relative Bias
FUNCTION GetBiasFloat &
!
(grid1, grid2, maskReal, maskInteger, rbias) &
!
RESULT (bias)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid1
TYPE (grid_real), INTENT(IN) :: grid2
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger
LOGICAL, OPTIONAL, INTENT(IN) :: rbias


!Local declarations:
INTEGER (KIND = long) :: i, j
REAL (KIND = float) :: bias, mean
INTEGER             :: countCells
!---------------------------end of declarations--------------------------------

bias = 0.
mean = 0.
countCells = 0

!check grid1 and grid2 have the same coordinate reference system
 IF ( .NOT. CRSisEqual(grid1,grid2) ) THEN
      CALL Catch ('error', 'GridStatistics',  &
      'calculate Bias: ', argument = &
      'coordinate reference system of grid1 differs from grid2' )
END IF

IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate Bias: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
   
    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                bias = bias + ( grid1 % mat (i,j) - grid2 % mat (i,j))
                countCells = countCells + 1.
                mean = mean + grid2 % mat (i,j)
            END IF
        END DO
    END DO
    
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid1) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate RMSE: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF
  
    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                 bias = bias + ( grid1 % mat (i,j) - grid2 % mat (i,j))
                 countCells = countCells + 1.
                 mean = mean + grid2 % mat (i,j)
            END IF
        END DO
    END DO

ELSE
    
    DO j = 1, grid1 % jdim
        DO i = 1, grid1 % idim
            IF (grid1 % mat(i,j) /= grid1 % nodata) THEN
                bias = bias + ( grid1 % mat (i,j) - grid2 % mat (i,j))
                countCells = countCells + 1.
                mean = mean + grid2 % mat (i,j)
            END IF
        END DO
    END DO
END IF

bias = bias / countCells

IF (PRESENT(rbias)) THEN
  IF (rbias) THEN
    mean = mean / countCells
    bias = bias / mean
  END IF
END IF

RETURN

END FUNCTION GetBiasFloat


!==============================================================================
!| Description:
!   compute sum of grid_real eventually constrained to a mask
FUNCTION GetSumOfGridFloat &
!
(grid, maskReal, maskInteger) &
!
RESULT (sum)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
REAL (KIND = float) :: sum
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
sum = 0.
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                sum = sum + grid % mat (i,j)
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate sum: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata .AND. &
                grid % mat(i,j) /= grid % nodata ) THEN
                sum = sum + grid % mat (i,j)
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                sum = sum + grid % mat (i,j)
            END IF
        END DO
    END DO
END IF

RETURN

END FUNCTION GetSumOfGridFloat



!==============================================================================
!| Description:
!   compute sum of grid_integer eventually constrained to a mask
FUNCTION GetSumOfGridInteger &
!
(grid, maskReal, maskInteger) &
!
RESULT (sum)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_integer), INTENT(IN) :: grid
TYPE (grid_real), OPTIONAL,  INTENT(IN) :: maskReal
TYPE (grid_integer), OPTIONAL,  INTENT(IN) :: maskInteger


!Local declarations:
INTEGER (KIND = long) :: sum
INTEGER (KIND = long) :: i, j
!---------------------------end of declarations--------------------------------
sum = 0.
!check that grid and mask have the same coordinate reference system
IF (PRESENT (maskReal)) THEN
    IF ( .NOT. CRSisEqual(maskReal,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskReal % jdim
        DO i = 1, maskReal % idim
            IF (maskReal % mat(i,j) /= maskReal % nodata) THEN
                sum = sum + grid % mat (i,j)
            END IF
        END DO
    END DO
ELSE IF (PRESENT (maskInteger)) THEN
    IF ( .NOT. CRSisEqual(maskInteger,grid) ) THEN
        CALL Catch ('error', 'GridStatistics',  &
        'calculate mean: ', argument = &
        'coordinate reference system of mask differs from input grid' )
    END IF

    DO j = 1, maskInteger % jdim
        DO i = 1, maskInteger % idim
            IF (maskInteger % mat(i,j) /= maskInteger % nodata) THEN
                sum = sum + grid % mat (i,j)
            END IF
        END DO
    END DO

ELSE
    DO j = 1, grid % jdim
        DO i = 1, grid % idim
            IF (grid % mat(i,j) /= grid % nodata) THEN
                sum = sum + grid % mat (i,j)
            END IF
        END DO
    END DO
END IF

RETURN

END FUNCTION GetSumOfGridInteger


!==============================================================================
!| Description:
!   count number of cells different from nodata
FUNCTION CountCellsFloat &
!
(grid) &
!
RESULT (count)

IMPLICIT NONE

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

!Local declaration
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: i, j

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

count = 0
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat(i,j) /= grid % nodata) THEN
      count = count + 1
    END IF
  END DO
END DO

RETURN
END FUNCTION CountCellsFloat


!==============================================================================
!| Description:
!   count number of cells different from nodata
FUNCTION CountCellsInteger &
!
(grid) &
!
RESULT (count)

IMPLICIT NONE

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

!Local declaration
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: i, j

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

count = 0
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat(i,j) /= grid % nodata) THEN
      count = count + 1
    END IF
  END DO
END DO

RETURN
END FUNCTION CountCellsInteger


!==============================================================================
!| Description:
!   performs a linear transformation on the original data values.
!   Suppose that `minX` and `maxX` are the minimum and maximum of feature X.
!   We would like to map interval `[minX, maxX]` into a new interval
!   `[new_minX, new_maxX]`. Consequently, every value `v` from the original
!   interval will be mapped into value `new_v` following formula:
! 
!  ` new_v = (v - minX) / (maxX - minX) * (new_maxX - new_minX) + new_minX`
!
!  If `new_minX` and `new_maxX` are not given values are mapped 
!  to `[0,1]` interval
! ***
!  References:<br>
!  Hann, J., Kamber, M. (2000). Data Mining: Concepts and Techniques. 
!    Morgan Kaufman Publishers. 
SUBROUTINE MinMaxNormalizationFloat &
!
(gridIn, gridOut, min, max) 

IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real), INTENT(IN) :: gridIn
REAL (KIND = float), OPTIONAL, INTENT(IN) :: min
REAL (KIND = float), OPTIONAL, INTENT(IN) :: max

!Arguments with intent(inout):
TYPE (grid_real), INTENT(INOUT) :: gridOut

!Local declaration
INTEGER (KIND = short) :: i, j
REAL (KIND = float) :: old_min, old_max
REAL (KIND = float) :: new_min, new_max

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

!set new_min and new_max
IF (PRESENT(min)) THEN
  new_min = min
ELSE
  new_min = 0.
END IF

IF (PRESENT(max)) THEN
  new_max = max
ELSE
  new_max = 0.
END IF

!get min and max of gridIn
old_min = GetMin (gridIn)
old_max = GetMax (gridIn)

!normalize grid. gridout is supposed to be initialized outside the subroutine
DO i = 1, gridIn % idim
  DO j = 1, gridOut % jdim
    IF (gridIn % mat (i,j) /= gridIn % nodata) THEN
       gridOut % mat (i,j) = (gridIn % mat (i,j) -  old_min) /  &
                             (old_max - old_min) * (new_max - new_min) + new_min
    END IF
  END DO

END DO

RETURN
END SUBROUTINE MinMaxNormalizationFloat


!==============================================================================
!| Description:
!   also called zero-mean normalization. The values of attribute `X` are 
!   normalized using the mean and standard deviation of `X`. A new value 
!   `new_v` is obtained using the following expression:
!
!   `new_v = (v - mX) / sX`
!
!   where `mX` and `sX` are the mean and standard deviation of attribute `X`, 
!   respectively. After zero-mean normalizing each feature will have a mean
!   value of 0. Also, the unit of each value will be the number of (estimated) 
!   standard deviations away from the (estimated) mean.
SUBROUTINE ZscoreNormalizationFloat &
!
(gridIn, gridOut) 

IMPLICIT NONE

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

!Arguments with intent(inout):
TYPE (grid_real), INTENT(INOUT) :: gridOut

!Local declaration
INTEGER (KIND = short) :: i, j
REAL (KIND = float) :: mean, std

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

!get mean
mean = GetMean (gridIn)

!get standard deviation
std = getStDev (gridIn)

!normalize grid. gridout is supposed to be initialized outside the subroutine
DO i = 1, gridIn % idim
  DO j = 1, gridIn % jdim
    IF (gridIn % mat (i,j) /= gridIn % nodata) THEN
      gridOut % mat (i,j) = (gridIn % mat (i,j) - mean) / std 
    END IF
  END DO

END DO

RETURN

END SUBROUTINE ZscoreNormalizationFloat


!------------------------------------------------------------------------------
!| Description
!   return a vector of distinct values from a grid_integer
SUBROUTINE UniqueValuesOfGridInteger &
!
(grid, values)

IMPLICIT NONE

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

!arguments with intent (out):
INTEGER (KIND = short), INTENT (OUT), ALLOCATABLE :: values (:)

!Local declarations:
INTEGER (KIND = short), ALLOCATABLE :: vector (:)
INTEGER (KIND = short), ALLOCATABLE :: unique (:)
INTEGER (KIND = short) :: i = 0
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: min_val, max_val
!------------------------------------end of declarations-----------------------

! vectorize grid
CALL Grid2Vector (grid, vector)

!allocate temporary unique array
ALLOCATE ( unique ( SIZE(vector) ) )

!elaborate
min_val = MINVAL (vector) - 1
max_val = MAXVAL (vector)
DO WHILE (min_val<max_val)
    i = i + 1
    min_val = MINVAL (vector, mask = vector > min_val)
    unique(i) = min_val
END DO
ALLOCATE ( values (i), source = unique(1:i) )   !<-- Or, just use unique(1:i) 

!free memory
DEALLOCATE (vector)
DEALLOCATE (unique)


RETURN
END SUBROUTINE UniqueValuesOfGridInteger

END MODULE GridStatistics