!! Contains subroutines used by Meteo module
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.0 - 19th October 2017   
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 19/Oct/2017 | Original code |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module containing subroutines used by Meteo module
!
MODULE MeteoUtilities      

! Modules used: 

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


USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE GridLib, ONLY : &
  !Imported types:
  grid_real, grid_integer, &
  !Imported parameters:
  NET_CDF, &
  !Imported routines:
  NewGrid, GridDestroy

USE Chronos, ONLY : &
  !Imported types:
  DateTime, &
  !Imported operations:
  OPERATOR (+), &
  ASSIGNMENT (=)

USE GridOperations, ONLY: &
  !Imported routines
  GridConvert, GridResample, &
  GetIJ

USE ObservationalNetworks, ONLY: &
  !Imported definitions:
  ObservationalNetwork

USE SpatialInterpolation, ONLY: &
  !imported routines:
  NearestNeighbor, IDW
  
USE Kriging, ONLY: &
    !imported routines:
    Krige

IMPLICIT NONE 

! Global (i.e. public) Declarations: 


!Local (i.e. private) declarations 


!Public routines
PUBLIC :: ReadField
PUBLIC :: ShiftMeteoWithLapse
PUBLIC :: ShiftBackWithLapse
PUBLIC :: Offset
PUBLIC :: Scale
PUBLIC :: Interpolate

!Local routines

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

!==============================================================================
!| Description:
!   read a time varying field from a netcdf file
SUBROUTINE ReadField &
!
(filename, time, dtAggr, dtGrid, aggrType, field, varName, &
 stdName, cellsize, dem, demHiRes, lapse )

USE StringManipulation, ONLY : &
! Imported routines:
StringCompact

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *), INTENT (IN) :: filename !!name of netcdf file
TYPE (DateTime),     INTENT (IN) :: time  !!time of the variable to read
INTEGER,             INTENT (IN) :: dtAggr !!aggregation time interval
INTEGER,             INTENT (IN) :: dtGrid !!time interval of grid in netcdf file
CHARACTER (LEN = *), INTENT (IN) :: aggrType !!aggregation type. 'M' = mean, 'C' = cumulated, 'X' = maximum, 'N' = minimum

CHARACTER (LEN = *), OPTIONAL, INTENT (IN) :: varName !!name of the variable to read
CHARACTER (LEN = *), OPTIONAL, INTENT (IN) :: stdName !!name of the variable to read
REAL,                OPTIONAL, INTENT (IN) :: cellsize
TYPE (grid_real),    OPTIONAL, INTENT (IN) :: dem
TYPE (grid_real),    OPTIONAL, INTENT (IN) :: demHiRes
TYPE (grid_real),    OPTIONAL, INTENT (IN) :: lapse

!Arguments with intent (out):
TYPE (grid_real), INTENT (INOUT) :: field

!Local declarations:
TYPE (grid_real) :: gridTemp !!temporary grid 
TYPE (grid_real) :: gridtemp2 !!temporary grid
TYPE (grid_real) :: gridTemp3 !!temporary grid
INTEGER :: nGrid !! number of grid to be read in netcdf file
TYPE (DateTime) :: gridTime
INTEGER :: i,j,k
CHARACTER (LEN = 100) :: variableName
CHARACTER (LEN = 100) :: standardName
REAL :: size

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

variableName = ''
standardName = ''
IF (PRESENT (cellsize)) THEN
  size  = cellsize
ELSE
  size = 0.
END IF

!read grid in netcdf file
IF (PRESENT(stdName)) THEN
  CALL NewGrid (GridTemp, filename, NET_CDF, stdName=stdName, time = time)
  standardName = stdName
ELSE IF (PRESENT(varName)) THEN
  CALL NewGrid (GridTemp, filename, NET_CDF, variable=varName, time = time)
  variableName = varName
ELSE !read info in field
  IF (TRIM(StringCompact(field % standard_name)) /= '') THEN
    CALL NewGrid (GridTemp, filename, NET_CDF, stdName=field % standard_name, time = time)
    standardName = field % standard_name
  ELSE IF  (TRIM(StringCompact(field % var_name)) /= '') THEN
    CALL NewGrid (GridTemp, filename, NET_CDF, variable=field % var_name, time = time)
    variableName = field % var_name
  ELSE
    CALL Catch ('error', 'MeteoUtilities',   &
			    'missing standard or variable name in grid while calling ReadField' )
  END IF
END IF

!compute number of grid to be read in netcdf file
nGrid = INT (dtAggr / dtGrid)

gridTime = time
!read other grids in netcdf file
DO k = 1, nGrid - 1 
  gridTime = gridTime + dtGrid
   
  IF (PRESENT(stdName)) THEN
      CALL NewGrid (gridTemp3, filename, NET_CDF, stdName=stdName, time = gridTime)
  ELSE IF (PRESENT(varName)) THEN
      CALL NewGrid (gridTemp3, filename, NET_CDF, variable=varName, time = gridTime)
  ELSE !read info in field
      IF (TRIM(StringCompact(field % standard_name)) /= '') THEN
        CALL NewGrid (gridTemp3, filename, NET_CDF, stdName=field % standard_name, time = gridTime)
      ELSE IF  (TRIM(StringCompact(field % var_name)) /= '') THEN
        CALL NewGrid (gridTemp3, filename, NET_CDF, variable=field % var_name, time = gridTime)
      ELSE
        CALL Catch ('error', 'MeteoUtilities',   &
			        'missing standard or variable name in grid while calling ReadField' )
      END IF
  END IF
  
  DO i = 1, gridTemp % idim
    DO j = 1, gridTemp % jdim
      IF (gridTemp % mat (i,j) /= gridTemp % nodata ) THEN
         SELECT CASE (aggrType)  
           CASE ('M','C') !mean or cumulated
              gridTemp % mat(i,j) = gridTemp % mat(i,j) + gridTemp3 % mat(i,j)
           CASE ('N') !minimum
              gridTemp % mat (i,j) = MIN (gridTemp % mat (i,j), gridTemp3 % mat(i,j))
           CASE ('X') !maximum
              gridTemp % mat (i,j) = MAX (gridTemp % mat (i,j), gridTemp3 % mat(i,j))
         END SELECT
      END IF
    END DO
  END DO
  CALL GridDestroy (gridTemp3) 
END DO

SELECT CASE (aggrType)
  CASE ('M') !mean
    DO i = 1, gridTemp % idim
      DO j = 1, gridTemp % jdim
        IF (gridTemp % mat (i,j) /= gridTemp % nodata ) THEN
          gridTemp % mat (i,j) = gridTemp % mat (i,j) / REAL (nGrid)
        END IF
      END DO
    END DO
END SELECT

!update attribute
  field % next_time = gridTemp % next_time
  field % reference_time = gridTemp % reference_time
  field % current_time = gridTemp % current_time 
  field % time_index = gridTemp % time_index
  field % time_unit = gridTemp % time_unit
  field % var_name = gridTemp % var_name
  field % standard_name = gridTemp % standard_name
  field % file_name = gridTemp % file_name
  
!coordinate conversion
IF ( .NOT. gridTemp % grid_mapping == field % grid_mapping) THEN
  !set Coordinate reference system of temporary grid
  gridTemp2 % grid_mapping = field % grid_mapping
  !convert coordinate
  IF (size > 0.) THEN
      CALL GridConvert (gridTemp, gridTemp2, cellsize = size)
  ELSE
       CALL GridConvert (gridTemp, gridTemp2)
  END IF

  !apply lapse rate 
  IF (PRESENT (dem) .AND. PRESENT (demHiRes) .AND. PRESENT (lapse) ) THEN
    DO i = 1, gridTemp2 % idim
     DO j = 1, gridTemp2 % jdim
        IF (gridTemp2 % mat(i,j) /= gridTemp2 % nodata) THEN
              gridTemp2 % mat (i,j) = gridTemp2 % mat (i,j) + &
                      ( demHiRes % mat(i,j) - dem % mat (i,j) ) * &
                      lapse % mat(i,j) 
              
        END IF
     END DO
    END DO
  END IF


  !resample grid
  CALL GridResample (gridTemp2, field)
  CALL GridDestroy (gridTemp)
  CALL GridDestroy (gridTemp2)
ELSE
  !resample grid
  CALL GridResample (gridTemp, field)
  CALL GridDestroy (gridTemp) 
END IF

RETURN


END SUBROUTINE ReadField

!==============================================================================
!| Description:
!   apply offset to a grid_real
SUBROUTINE Offset &
!
(grid, off)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: off

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

!Local declarations
INTEGER :: i,j
!------------end of declaration------------------------------------------------

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

END SUBROUTINE Offset

!==============================================================================
!| Description:
!   apply scale factor to a grid_real
SUBROUTINE Scale &
!
(grid, sc)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: sc

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

!Local declarations
INTEGER :: i,j
!------------end of declaration------------------------------------------------

DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
      grid % mat (i,j) = grid % mat (i,j) * sc
    END IF
  END DO
END DO

END SUBROUTINE Scale


!==============================================================================
!| Description:
!   shift meteo observations to reference elevation applying lapse rate
SUBROUTINE ShiftMeteoWithLapse &
!
(input, lapse, refelev, output, dt)

IMPLICIT NONE

!Arguments with intent(in):
TYPE (ObservationalNetwork), INTENT(IN) :: input  !!actual station network
TYPE (grid_real), INTENT(IN) :: lapse  !! lapse rate grid
REAL (KIND = float), INTENT(IN) :: refelev !!reference elevation
INTEGER, OPTIONAL, INTENT(IN) :: dt !! used when lapse rate is a flux

!Arguments with intent(inout):
TYPE (ObservationalNetwork), INTENT(INOUT) :: output !!station network at reference elevation

!local eclarations:
INTEGER :: i, j, r, c
REAL (KIND = float) :: x, y
LOGICAL            :: check
REAL (KIND = float) :: deltat

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

IF (PRESENT (dt)) THEN
   deltat = dt
ELSE
   deltat = 1.
END IF

!dato al livello di riferimento
	DO i = 1, input % countObs
		IF (input % obs(i) % value == input % nodata) THEN
			output % obs(i) % value = output % nodata
		ELSE
		    x = input % obs(i) % xyz % easting
		    y = input % obs(i) % xyz % northing
	        CALL GetIJ (x, y, lapse, r, c, check)
	        IF (check) THEN ! lapse rate is defined at location x,y
	            IF (lapse % mat (r,c) /= lapse % nodata) THEN 
			              output % obs(i) % value = input % obs(i) % value + &
                          (refelev - input % obs(i) % xyz % elevation) * &
                           lapse % mat(r,c) * deltat
			        ELSE
			            output % obs(i) % value = input % obs(i) % value
			        END IF
			    ELSE
			        output % obs(i) % value = input % obs(i) % value
			END IF
		END IF
	END DO


RETURN
END SUBROUTINE ShiftMeteoWithLapse


!==============================================================================
!| Description:
!   interpolate site data to space
SUBROUTINE Interpolate &
!
(network, method, near, idw_power, anisotropy, varmodel, lags, maxlag, grid, devst)

IMPLICIT NONE

!arguments with intent in:
TYPE (ObservationalNetwork), INTENT(IN) :: network
INTEGER (KIND = short),      INTENT(IN) :: method
INTEGER (KIND = short),      INTENT(IN) :: near
REAL (KIND = float),         INTENT(IN) :: idw_power
INTEGER (KIND = short),      INTENT(IN) :: anisotropy
INTEGER (KIND = short),      INTENT(IN) :: varmodel
INTEGER (KIND = short),      INTENT(IN) :: lags
REAL (KIND = float),         INTENT(IN) :: maxlag

!arguments with intent inout:
TYPE (grid_real), INTENT(INOUT) :: grid 
TYPE (grid_real), INTENT(INOUT) :: devst       

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


SELECT CASE (method)
    CASE (1) !thiessen
		    CALL NearestNeighbor (network = network, &
                    grid = grid )
    CASE (2) !IDW
        CALL IDW (network = network, &
                  grid = grid, &
                  method = 1, &
                  power = idw_power, &
                  near = near)
    CASE (3) !Kriging
        CALL Krige (sitedata = network, &
                        anisotropy = anisotropy, &
                        nlags = lags, &
                        maxlag = maxlag, &
                        neighbors = near, &
                        varModel = varmodel, &
                        grid = grid, &
                        gridvar = devst )
END SELECT

RETURN
END SUBROUTINE Interpolate


!==============================================================================
!| Description:
!   Shift back interpolated field to terrain surface
SUBROUTINE ShiftBackWithLapse &
!
(grid, dem, lapse, refelev, dt)

IMPLICIT NONE

!arguments with intent in:
TYPE (grid_real), INTENT(IN) :: dem
TYPE (grid_real), INTENT(IN) :: lapse
REAL (KIND = float), INTENT(IN) :: refelev !!ereference elevation
INTEGER (KIND = short), OPTIONAL, INTENT(IN) :: dt

!arguments with intent inout:
TYPE (grid_real), INTENT(INOUT) :: grid

!local declarations:
INTEGER (KIND = short) :: i,j
INTEGER (KIND = short) :: deltat

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

IF (PRESENT (dt)) THEN
   deltat = dt
ELSE
   deltat = 1.
END IF

DO i = 1, grid % idim
    DO j = 1, grid % jdim
			   
				IF (grid % mat(i,j) == grid % nodata) CYCLE
				 
				grid % mat(i,j) = grid % mat(i,j) - &
				              (refelev - dem % mat(i,j)) * &
				              lapse % mat(i,j) * deltat
    END DO
END DO

RETURN
END SUBROUTINE ShiftBackWithLapse


END MODULE MeteoUtilities