!! Interpolate sparse point measurements to regular grid  
!|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.5 - 18th December 2018  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 03/Jun/2011 | Original code |
! | 1.1      | 01/apr/2013 | added Barnes Objective Analysis interpolation |
! | 1.2      | 20/Feb/2016 | added Kriging interpolation |
! | 1.3      | 09/Mar/2016 | WindMicromet added |
! | 1.4      | 01/Mar/2018 | WindGonzalezLongatt added |
! | 1.5      | 18/Dec/2018 | Kriging was moved to a specific module |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
!   set of generic routines to convert sparse point measurements to regular grid  
!   Methods implemented:
!
! * Inverse Distance Weighted
! * Nearest Neighbor (Thiessen polygons)
! * Barnes Objective Analysis
! * Kriging with automatic fitting of linear semivariogram (moved to a specific module with several semivariogram models)
! * Wind with topographic correction Micromet
! * Wind with orographic correction (Gonzales-Longatt etal., 2015)
!
! References:
!
!  González-Longatt, F., Medina, H., Serrano González, J., Spatial interpolation 
!     and orographic correction to estimate wind energy resource in Venezuela.
!     Renewable and Sustainable Energy Reviews, 48, 1-16, 2015.
!
MODULE SpatialInterpolation


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

USE LogLib, ONLY: &
  !Imported routines:
  Catch

USE ObservationalNetworks, ONLY: &
  !Imported definitions:
  ObservationalNetwork, &
  !Imported routines:
  ActualObservations, CopyNetwork, &
  DestroyNetwork

USE GridLib, ONLY: &
  !Imported definitions:
  grid_real, grid_integer, &
  !Imported routines:
  NewGrid, GridDestroy, &
  !Imported operators:
  ASSIGNMENT( = )

USE GridOperations, ONLY: &
  !Imported routines:
  GetXY, CellArea, CellIsValid


USE GridStatistics, ONLY: &
  !Imported routines:
  CountCells
  
USE GeoLib, ONLY: &
  !Imported variables:
  point1, point2, &
  !Imported routines:
  Distance, &
  !Imported operators:
  ASSIGNMENT( = )


IMPLICIT NONE

! Global (i.e. public) Declarations: 

! Global Procedures:

PUBLIC :: IDW
PUBLIC :: NearestNeighbor
PUBLIC :: BarnesOI
PUBLIC :: WindMicromet
PUBLIC :: WindGonzalezLongatt
!PUBLIC :: Zonal TODO


!Private procedures
PRIVATE :: Sort
PRIVATE :: Inverse

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

!==============================================================================
!| Description:
!   Inverse Distance Weighted interpolation. Accept as optional argument
!   the power parameter and the number of near observations to 
!   included in interpolation. Can use Shepard's method (Shepard 1968) 
!   or Franke & Nielson, 1980
!
! References:
!
!   Shepard, D. (1968) A two-dimensional interpolation function for 
!   irregularly-spaced data, Proc. 23rd National Conference ACM, ACM, 517-524. 
!
!   Franke, R. and G. Nielson (1980), Smooth interpolation of large sets
!   of scattered data. International Journal of Numerical Methods 
!   in Engineering, 15, 1691-1704.
SUBROUTINE IDW &
!
(network, grid, method, power, near)


IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network
INTEGER (KIND = short), INTENT (IN) :: method
REAL (KIND = float), OPTIONAL, INTENT (IN) :: power
INTEGER (KIND = short), OPTIONAL, INTENT (IN) :: near

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

!Local declarations
TYPE (ObservationalNetwork) :: activeNetwork
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: i, j, k, s, t
INTEGER (KIND = short) :: nearPoints
REAL (KIND = float)    :: actPower
REAL (KIND = float),POINTER   :: dist(:,:)  
REAL (KIND = float)   :: distIJ
REAL (KIND = float)   :: denom 
REAL (KIND = float)   :: weight 
!-----------------------end of declarations------------------------------------

!check for available measurements
CALL ActualObservations (network, count, activeNetwork)

!allocate distances vector
ALLOCATE (dist (activeNetwork % countObs + 1,2) )

!set points
point1 % system = grid % grid_mapping
point2 % system = grid % grid_mapping

!set near
IF (PRESENT (near)) THEN
  IF (activeNetwork % countObs <= near) THEN 
    nearPoints = activeNetwork % countObs
  ELSE
    nearPoints = near
  END IF
ELSE
  nearPoints = activeNetwork % countObs !use all data
END IF

!set power
IF (PRESENT (power)) THEN
  actPower = power
ELSE
  actPower = 2. !default to square
END IF

DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
      dist = -9999.99
      !compute distance from cell center to observations
      CALL GetXY (i,j,grid, point1 % easting, point1 % northing)
      DO k = 1, activeNetwork % countObs        
         point2 % northing = activeNetwork % obs (k) % xyz % northing  
         point2 % easting = activeNetwork % obs (k) % xyz % easting
         distIJ = Distance (point1,point2)
                  
		 !search for neighbours
		 DO s = 1, nearPoints
			IF ( dist(s,1) == -9999.99 .OR. &
			     dist(s,1) > distIJ       ) THEN
				DO t = nearPoints, s, -1
					dist(t+1,1) = dist(t,1)
					dist(t+1,2) = dist(t,2)
				END DO
				dist(s,1) = distIJ
				dist(s,2) = activeNetwork % obs(k) % value
				EXIT
			END IF
		 END DO
	  END DO	 
	  !compute denominator
	  denom = 0.0
	  IF (method == 1) THEN !Shepard's method
	      DO s = 1, nearPoints
	 	      IF (dist(s,1) <= 1.) THEN
			     !observation in cell, set minimum distance to avoid dividing by zero
			     dist(s,1) = 1.
		      END IF
		      denom = denom + dist(s,1) ** (-actPower)
	      END DO
	      ! weighted value
	      grid % mat(i,j) = 0
	      DO s = 1, nearPoints
	        weight = dist(s,1) ** (-actPower) / denom
	        grid % mat(i,j) = grid % mat(i,j) + weight * dist(s,2)
	      END DO
	   ELSE !Franke & Nielson method
	      DO s = 1, nearPoints
	 	      IF (dist(s,1) <= 1.) THEN
			     !observation in cell, set minimum distance to avoid dividing by zero
			     dist(s,1) = 1.
		      END IF
		      denom = denom + ( (dist(nearPoints,1) - dist(s,1)) / (dist(nearPoints,1) * dist(s,1)) ) ** actPower
	      END DO
	      ! weighted value
	      grid % mat(i,j) = 0
	      DO s = 1, nearPoints
	        weight = ( (dist(nearPoints,1) - dist(s,1)) / (dist(nearPoints,1) * dist(s,1)) ) ** actPower / denom
	        grid % mat(i,j) = grid % mat(i,j) + weight * dist(s,2)
	      END DO
	   END IF		 
    END IF
  END DO
END DO

DEALLOCATE(dist)
DEALLOCATE(activeNetwork % obs)

RETURN

END SUBROUTINE IDW



!==============================================================================
!| Description:
!   This subroutine implements the method used in the MICROMET program (see
!   reference). Wind speed is interpolated accounting for wind direction
!   and an empirical weigth that considers slope and curvature (topographic
!   effect)
!
! References:
!
!  Liston, G. E., Elder, K.,  2006. A Meteorological Distribution System 
!  for High-Resolution Terrestrial Modeling (MicroMet). 
!  Journal of Hudrometeorology, 7, 217-234. 
SUBROUTINE WindMicromet &
!
(speed, dir, slope, curvature, slope_az, slopewt, curvewt, grid, winddir)

USE Units, ONLY: &
!Imported parameters
degToRad, radToDeg, pi

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: speed !windspeed observations at stations 
TYPE (ObservationalNetwork), INTENT(IN) :: dir !wind direction observations at stations 
TYPE (grid_real), INTENT(IN) :: slope !terrain slope grid
TYPE (grid_real), INTENT(IN) :: curvature !terrain curvature grid
TYPE (grid_real), INTENT(IN) :: slope_az !terrain slope azimuth
REAL (KIND = float),  INTENT(IN) :: slopewt !slope weighting factor
REAL (KIND = float),  INTENT(IN) :: curvewt !curvature weighting factor

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid !wind speed grid [m/s]
TYPE (grid_real), OPTIONAL, INTENT (INOUT) :: winddir !wind direction [deg]


!Local declarations
TYPE (ObservationalNetwork) :: temp_speed, temp_dir, active_speed, active_dir
TYPE (ObservationalNetwork) :: uwind, vwind
INTEGER :: count_speed, count_dir, i, j
TYPE (grid_real) :: uwind_grid !zonal wind [m/s]
TYPE (grid_real) :: vwind_grid !meridional wind [m/s]
TYPE (grid_real) :: winddir_grid !wind direction [rad]
TYPE (grid_real) :: wind_slope !slope in the direction of the wind [rad]
REAL (KIND = float) :: wslope_max ![rad]
REAL (KIND = float) :: windwt !wind weighting factor
REAL (KIND = float) :: thetad !diverting factor
REAL (KIND = float),  PARAMETER :: windspd_min = 0.00001 ![m/s]


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

!check for available measurements. Both speed and direction must exist
!first check nodata. If only speed or direction are missing, set both to missing 
! assume dir and speed sations are in memory in the same order
CALL CopyNetwork (speed,temp_speed)
CALL CopyNetwork (dir,temp_dir)

DO i = 1, temp_speed %countObs
  IF (temp_speed % obs (i) % value == temp_speed % nodata) THEN
      !set dir to nodata
      temp_dir % obs (i) % value = temp_dir % nodata
  ELSE IF (temp_dir % obs (i) % value == temp_dir % nodata) THEN
      !set speed to nodata
      temp_speed % obs (i) % value = temp_speed % nodata
  END IF
END DO

!now remove nodata
CALL ActualObservations (temp_speed, count_speed, active_speed)
CALL ActualObservations (temp_dir, count_dir, active_dir)

!initialize zonal and meridional component network
CALL CopyNetwork (active_speed, uwind)
CALL CopyNetwork (active_speed, vwind)

!initialize grid
CALL NewGrid (uwind_grid, grid)
CALL NewGrid (vwind_grid, grid)
CALL NewGrid (winddir_grid, grid)
CALL NewGrid (wind_slope, grid)

!compute u and v at stations
DO i = 1, active_speed % countObs
  uwind % obs (i) % value = - active_speed % obs (i) % value * &
                          SIN ( active_dir % obs (i) % value * degToRad)
  vwind % obs (i) % value = - active_speed % obs (i) % value * &
                          COS ( active_dir % obs (i) % value * degToRad)
END DO

!interpolate u and v grid 
CALL IDW (network=uwind, grid=uwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=uwind, grid=uwind_grid, gammapar = 0.6)
     
CALL IDW (network=vwind, grid=vwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=vwind, grid=vwind_grid, gammapar = 0.6)
   
!compute wind direction and windspeed grid 
DO i = 1, grid % idim
  DO j = 1, grid % jdim
      IF (grid % mat (i,j) /= grid % nodata) THEN
           
        ! Some compilers do not allow both u and v to be 0.0 in
        !   the atan2 computation.
        IF (ABS(uwind_grid % mat (i,j)) < 1e-10) &
              uwind_grid % mat (i,j) = 1e-10
           
        winddir_grid %mat (i,j) = &
            ATAN2(uwind_grid % mat (i,j),vwind_grid % mat (i,j))
  
        IF (winddir_grid %mat (i,j) >= pi) THEN
          winddir_grid % mat (i,j) = winddir_grid % mat (i,j) - pi
        ELSE
          winddir_grid % mat (i,j) = winddir_grid % mat (i,j) + pi
        END IF
           
        !windspeed
        grid % mat (i,j) = &
            SQRT(uwind_grid % mat (i,j)**2 + vwind_grid % mat (i,j)**2)
            
      END IF
  END DO
END DO

! Compute the slope in the direction of the wind.
DO i = 1, grid % idim
  DO j = 1, grid % jdim
      IF (grid % mat (i,j) /= grid % nodata) THEN
     
      wind_slope % mat (i,j) = slope % mat (i,j) * &
            COS ( winddir_grid % mat (i,j) - slope_az % mat (i,j) )
         
      END IF
  END DO
END DO

! Scale the wind slope such that the max abs(wind slope) has a value
! of abs(0.5).  Include a 1 mm slope in slope_max to prevent
! divisions by zero in flat terrain where the slope is zero.
wslope_max = 0.0 + 0.001
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
        wslope_max = MAX (wslope_max,ABS(wind_slope % mat (i,j)))
    END IF
  END DO
END DO
      
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
        wind_slope % mat (i,j) = wind_slope % mat (i,j) / (2.0 * wslope_max)
    END IF
  END DO
END DO

! Calculate the wind speed and direction adjustments.  The
! curvature and wind_slope values range between -0.5 and +0.5.
! Valid slopewt and curvewt values are between 0 and 1, with
! values of 0.5 giving approximately equal weight to slope and
! curvature.  I suggest that slopewt and curvewt be set such
! that slopewt + curvewt = 1.0.  This will limit the total
! wind weight to between 0.5 and 1.5 (but this is not required).
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
            
      ! Compute the wind weighting factor.
      windwt = 1.0 + slopewt * wind_slope % mat (i,j) + &
      curvewt * curvature % mat (i,j)
      !if (windwt>0) write(*,*) windwt,  wind_slope % mat (i,j), curvature % mat (i,j)
      !pause
            
      ! Generate the terrain-modified wind speed.
      grid % mat (i,j) = windwt * grid % mat (i,j)
      
      !compute diverting factor
      thetad = - 0.5 * wind_slope % mat (i,j) * &
          SIN ( 2. * (slope_az % mat (i,j) - winddir_grid % mat (i,j)) )
      
      IF (PRESENT (winddir)) THEN
      
         winddir % mat (i,j) = ( winddir_grid % mat (i,j) + thetad ) * radToDeg
         
      END IF   
              
    END IF
  END DO
END DO

!deallocate memory
CALL DestroyNetwork (temp_speed)
CALL DestroyNetwork (temp_dir)
CALL DestroyNetwork (active_speed)
CALL DestroyNetwork (active_dir)
CALL DestroyNetwork (uwind)
CALL DestroyNetwork (vwind)

CALL GridDestroy (uwind_grid)
CALL GridDestroy (vwind_grid)
CALL GridDestroy (winddir_grid)
CALL GridDestroy (wind_slope)

RETURN
END SUBROUTINE WindMicromet


!==============================================================================
!| Description:
!   This subroutine implements the method presented by Gonzalez-Longatt et al.
!   (2015). Zonal and meridional components are computed and then orographic 
!   correction is applied. The two components are then re-composed to
!   provide final result.
!
! References:
!
!  González-Longatt, F., Medina, H., Serrano González, J., Spatial interpolation 
!     and orographic correction to estimate wind energy resource in Venezuela.
!     Renewable and Sustainable Energy Reviews, 48, 1-16, 2015.
!
SUBROUTINE WindGonzalezLongatt &
!
(speed, dir, dem, grid, winddir)

USE Units, ONLY: &
!Imported parameters
degToRad, radToDeg, pi

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: speed !windspeed observations at stations 
TYPE (ObservationalNetwork), INTENT(IN) :: dir !wind direction observations at stations 
TYPE (grid_real), INTENT(IN) :: dem !digital elevation model

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid !wind speed grid [m/s]
TYPE (grid_real), OPTIONAL, INTENT (INOUT) :: winddir !wind direction [deg]


!Local declarations
TYPE (ObservationalNetwork) :: temp_speed, temp_dir, active_speed, active_dir
TYPE (ObservationalNetwork) :: uwind, vwind
INTEGER :: count_speed, count_dir, i, j
TYPE (grid_real) :: uwind_grid !zonal wind [m/s]
TYPE (grid_real) :: vwind_grid !meridional wind [m/s] !OK
REAL (KIND = float) :: vxcalc !zonal wind with orograpcic correction applied [m/s]
REAL (KIND = float) :: vycalc !meridional wind with orograpcic correction applied [m/s]
REAL (KIND = float) :: corr1, corr2, cellsize
LOGICAL :: validNorth, validSouth, validEast, validWest


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

!check for available measurements. Both speed and direction must exist
!first check nodata. If only speed or direction are missing, set both to missing 
! assume dir and speed stations are in memory in the same order
CALL CopyNetwork (speed,temp_speed)
CALL CopyNetwork (dir,temp_dir)

DO i = 1, temp_speed %countObs
  IF (temp_speed % obs (i) % value == temp_speed % nodata) THEN
      !set dir to nodata
      temp_dir % obs (i) % value = temp_dir % nodata
  ELSE IF (temp_dir % obs (i) % value == temp_dir % nodata) THEN
      !set speed to nodata
      temp_speed % obs (i) % value = temp_speed % nodata
  END IF
END DO

!now remove nodata
CALL ActualObservations (temp_speed, count_speed, active_speed)
CALL ActualObservations (temp_dir, count_dir, active_dir)

!initialize zonal and meridional component network
CALL CopyNetwork (active_speed, uwind)
CALL CopyNetwork (active_speed, vwind)

!initialize grid
CALL NewGrid (uwind_grid, grid)
CALL NewGrid (vwind_grid, grid)

!compute u and v at stations
DO i = 1, active_speed % countObs
  uwind % obs (i) % value = - active_speed % obs (i) % value * &
                          SIN ( active_dir % obs (i) % value * degToRad)
  vwind % obs (i) % value = - active_speed % obs (i) % value * &
                          COS ( active_dir % obs (i) % value * degToRad)
END DO

!interpolate u and v grid 
CALL IDW (network=uwind, grid=uwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=uwind, grid=uwind_grid, gammapar = 0.6)
     
CALL IDW (network=vwind, grid=vwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=vwind, grid=vwind_grid, gammapar = 0.6)


!apply orographic correction
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN

       !set cellsize
       cellsize = (CellArea (grid,i,j)  ) ** 0.5 ![m]
       
       !check for neighbour cells. Nodata and out
       !of grid cells must not be considered
       !for orographic correction     
       validNorth = CellIsValid (i-1, j, grid)
       validSouth = CellIsValid (i+1, j, grid)
       validEast  = CellIsValid (i, j+1, grid)
       validWest  = CellIsValid (i, j-1, grid)

       !correct zonal (along x) component
       IF (validEast .AND. validWest)  THEN !apply correction where east and west data are defined

          !west correction factor component
          corr1 = (dem % mat (i,j)  -  dem % mat (i,j-1)) * &
                (uwind_grid % mat (i,j-1) - uwind_grid % mat (i,j)  ) / cellsize
          
          !east correction factor component
          corr2 = (dem % mat (i,j+1)  -  dem % mat (i,j)) * &
                (uwind_grid % mat (i,j) - uwind_grid % mat (i,j+1)  ) / cellsize
          vxcalc = uwind_grid % mat(i,j) + 0.5 * ( corr1 + corr2)
       ELSE !do not apply correction
          vxcalc = uwind_grid % mat(i,j)
       END IF

       !correct meridional (along y) component
       IF (validNorth .AND. validSouth)  THEN  !apply correction north and south data are defined

          !south correction factor component
          corr1 = (dem % mat (i,j)  -  dem % mat (i+1,j)) * &
                (vwind_grid % mat (i+1,j) - vwind_grid % mat (i,j)  ) / cellsize
          
          !north correction factor component
          corr2 = (dem % mat (i-1,j)  -  dem % mat (i,j)) * &
                (vwind_grid % mat (i,j) - vwind_grid % mat (i-1,j)  ) / cellsize
          vycalc = vwind_grid % mat(i,j) + 0.5 * ( corr1 + corr2)

       ELSE !do not apply correction
          vycalc = vwind_grid % mat(i,j)
       END IF  
       
       !compute wind direction if required
       IF (PRESENT (winddir)) THEN
           ! Some compilers do not allow both u and v to be 0.0 in
           ! the atan2 computation.
           IF ( ABS(vxcalc ) < 1e-10) vxcalc = 1e-10
           
           winddir % mat (i,j) = ATAN2(vxcalc,vycalc)
  
           IF (winddir % mat (i,j) >= pi) THEN
              winddir % mat (i,j) = winddir % mat (i,j) - pi
           ELSE
              winddir % mat (i,j) = winddir % mat (i,j) + pi
           END IF
           !convert radians to deg
           winddir % mat (i,j) = winddir % mat (i,j) * radToDeg
       ENDIF
           
        !windspeed
        grid % mat (i,j) = SQRT(vxcalc**2. + vycalc**2.)
        
    END IF
  END DO
END DO


!deallocate memory
CALL DestroyNetwork (temp_speed)
CALL DestroyNetwork (temp_dir)
CALL DestroyNetwork (active_speed)
CALL DestroyNetwork (active_dir)
CALL DestroyNetwork (uwind)
CALL DestroyNetwork (vwind)

CALL GridDestroy (uwind_grid)
CALL GridDestroy (vwind_grid)

RETURN
END SUBROUTINE WindGonzalezLongatt




!==============================================================================
!| Description:
!   The nearest neighbor algorithm selects the value of the nearest point 
!   and does not consider the values of neighboring points at all, 
!   yielding a piecewise-constant interpolant
!
! References:
!
!  A.H. Thiessen. 1911. Precipitation averages for large areas. 
!  Monthly Weather Review, 39(7): 1082-1084. 
SUBROUTINE NearestNeighbor &
!
(network, grid, weights, gridPolygons)

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network

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

!Arguments with intent (out):
REAL, OPTIONAL, INTENT(OUT), POINTER      :: weights(:)
TYPE (grid_integer), OPTIONAL, INTENT(OUT) :: gridPolygons

!Local declarations
TYPE (ObservationalNetwork) :: activeNetwork
INTEGER (KIND = short) :: count
INTEGER(KIND = short), ALLOCATABLE  :: polygons(:,:)
REAL (KIND = float) :: dist, dist_min
INTEGER (KIND = short) :: i, j, k

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

!check for available measurements
CALL ActualObservations (network, count, activeNetwork)

!set points
point1 % system = grid % grid_mapping
point2 % system = grid % grid_mapping

!allocate polygons matrix
ALLOCATE (polygons (grid % idim, grid % jdim) )
polygons = -1


DO i = 1, grid % idim
	col:DO j = 1, grid % jdim
		IF (grid % mat(i,j) == grid % nodata) THEN
			CYCLE COL
		ELSE
			!initialize
			grid % mat(i,j) = 0.0
		END IF
		dist = -9999.99
		dist_min = -9999.99
		
		!compute distance from cell center to observations
		CALL GetXY (i,j,grid,point1 % easting,point1 % northing)
        DO k = 1, activeNetwork % countObs 
           point2 % northing = activeNetwork % obs (k) % xyz % northing  
           point2 % easting = activeNetwork % obs (k) % xyz % easting
           dist = Distance(point1,point2)
		   !search for nearest points
			IF ( dist_min == -9999.99 .OR. dist_min > dist ) THEN
					polygons(i,j) = k
				    dist_min = dist
			END IF
		END DO
	END DO col
END DO

!finalize interpolation
DO i = 1, grid % idim
	DO j = 1, grid % jdim
		IF (grid % mat(i,j) == grid % nodata) CYCLE
		k = polygons (i,j)
		grid % mat(i,j) = activeNetwork % obs (k) % value
	END DO
END DO

IF (PRESENT (gridPolygons)) THEN
    CALL NewGrid (gridPolygons,grid)
	gridPolygons % mat = polygons
END IF

IF (PRESENT(weights)) THEN
    ALLOCATE(weights(activeNetwork % countObs))
    weights = 0.
	    DO i = 1, grid % idim
		    DO j = 1, grid % jdim
			    IF (grid % mat(i,j) == grid % nodata) CYCLE
			    weights(polygons(i,j)) = weights(polygons(i,j)) + 1
		    END DO
	    END DO
    weights = weights / CountCells(grid)
END IF

RETURN
END SUBROUTINE NearestNeighbor



!==============================================================================
!| Description:
!   The method constructs a grid of size determined by the distribution 
!   of the two dimensional data points. Using this grid, the function 
!   values are calculated at each grid point. To do this the method 
!   utilises a series of Gaussian functions, given a distance weighting 
!   in order to determine the relative importance of any given measurement 
!   on the determination of the function values. Correction passes are 
!   then made to optimise the function values, by accounting for the 
!   spectral response of the interpolated points.
!
! References:
!
!  Barnes, S. L., 1964: A technique for maximizing details in numerical 
!     map analysis. Journal of Applied Meteorology, 3, 395–409.
!
!  Koch, S., M. desJardins,and P. Kocin, 1983: An Interactive Barnes 
!    Objective Map Analysis Scheme for Use with Satellite and 
!    Convectional Data. Journal of Appl. Meteor., 22, 1487-1503
!
!  Maddox, R., 1980 : An objective technique for seaprating 
!    macroscale and mesoscale features in meteorological data. 
!    Mon. Wea. Rev., 108, 1108-1121
!
SUBROUTINE BarnesOI &
!
(network, grid, gammapar)

USE Units, ONLY: &
!Imported parameters
pi

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: network

REAL (KIND = float), OPTIONAL, INTENT(IN) :: gammapar

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

!local declarations:
REAL (KIND = float) :: dn         !!data spacing used in the analysis
REAL (KIND = float) :: gamma      !! smoothing parameter [0.2 - 1.0]
REAL (KIND = float) :: xkappa_1   !! Gauss constant for first pass
REAL (KIND = float) :: xkappa_2   !! Gauss constant for second pass
REAL (KIND = float) :: rmax_1     !! maximum scanning radii, for first
REAL (KIND = float) :: rmax_2     !! and second passes
REAL (KIND = float) :: anum_1     !! numerator, beyond scanning radius,
REAL (KIND = float) :: anum_2     !! for first and second passes
REAL (KIND = float) :: xa,ya      !! x, y coords of current station
REAL (KIND = float) :: xb,yb      !! x, y coords of current station
REAL (KIND = float) :: w1,w2      !! weights for Gauss-weighted average
REAL (KIND = float) :: wtot1,wtot2 !!sum of weights
REAL (KIND = float) :: ftot1,ftot2 !!accumulators for values, corrections
REAL (KIND = float) :: dsq !!square of distance between two points
REAL (KIND = float) :: dnMax, dnMin
REAL (KIND = float) :: nc, nr !! number of columns and rows in the grid
REAL (KIND = float) :: cellsize ![m]
REAL (KIND = float), ALLOCATABLE :: dvar (:) !!estimated error
TYPE (ObservationalNetwork) :: activeNetwork
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: i, j, k !! loop indexes

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

!check for available measurements
CALL ActualObservations (network, count, activeNetwork)

!Allocate dvar
ALLOCATE (dvar(activeNetwork % countObs))

!set points
point1 % system = grid % grid_mapping
point2 % system = grid % grid_mapping

!set gamma to default or user specified value
IF (PRESENT(gammapar)) THEN
  gamma = gammapar
ELSE
  gamma = 0.2 !default value
END IF

!get dn
!compute average cellsize.
cellsize = ( CellArea (grid, grid % idim / 2, grid % jdim / 2) ) ** 0.5 ![m]
nc = grid % jdim
nr = grid % idim
dnMax = ( (nr * cellsize) * (nc * cellsize) ) ** 0.5 * &
        ((1.0 + activeNetwork % countObs ** 0.5) / (activeNetwork % countObs - 1.0))

dnMin = ( (nr * cellsize) * (nc * cellsize) / activeNetwork % countObs ) ** 0.5

dn = 0.5 * (dnMin + dnMax)
!debug
dn = 6. * cellsize

! Compute the first and second pass values of the scaling parameter
! and the maximum scanning radius used by the Barnes scheme.
! Values above this maximum will use a 1/r**2 weighting.

! First-round values
	xkappa_1 = 5.052 * (2. * dn / Pi) ** 2.0

! Define the maximum scanning radius to have weight defined by
! wt = 1.0 x 10**(-30) = exp(-rmax_1/xkappa_1)
! Also scale the 1/r**2 wts so that when dsq = rmax, the wts match.
	rmax_1 = xkappa_1 * 30.0 * log(10.0)
	anum_1 = 1.E-30 * rmax_1
	
! Second-round values
	xkappa_2 = gamma * xkappa_1
	rmax_2 = rmax_1 * gamma
	anum_2 = 1.E-30 * rmax_2
	
! Scan each input data point and construct estimated error, dvar, at that point
outer_loop: DO i = 1, activeNetwork % countObs

        point1 % northing = activeNetwork % obs (i) % xyz % northing  
        point1 % easting  = activeNetwork % obs (i) % xyz % easting
		wtot1 = 0.0
		ftot1 = 0.0
		
   inner_loop: DO j = 1, activeNetwork % countObs
   
            point2 % northing = activeNetwork % obs (j) % xyz % northing  
            point2 % easting  = activeNetwork % obs (j) % xyz % easting
            dsq = Distance (point1,point2) ** 2.
            
            IF (dsq <= rmax_1) THEN
               w1 = exp((- dsq)/xkappa_1)
            ELSE !Assume a 1/r**2 weight
               w1 = anum_1/dsq
            END IF

			wtot1 = wtot1 + w1
			ftot1 = ftot1 + w1 * activeNetwork % obs (j) % value
    END DO inner_loop

!
!		if(wtot1==0.0) printf("stn wt totals zero\n");

        dvar(i) = activeNetwork % obs (i) % value - ftot1/wtot1

END DO outer_loop

! Grid-prediction loop.  Generate the estimate using first set of
! weights, and correct using error estimates, dvar, and second
! set of weights.

row_loop: DO i = 1, grid % idim
  col_loop: DO j = 1, grid % jdim
  
  IF (grid % mat (i,j) /= grid % nodata) THEN
      CALL GetXY (i,j,grid, point1 % easting, point1 % northing)

      ! Scan each input data point.
      ftot1 = 0.0
      wtot1 = 0.0
      ftot2 = 0.0
      wtot2 = 0.0
      
      stations_loop: DO k = 1, activeNetwork % countObs
          point2 % northing = activeNetwork % obs (k) % xyz % northing  
          point2 % easting  = activeNetwork % obs (k) % xyz % easting
          
          dsq = Distance (point1,point2) ** 2.

          IF (dsq<=rmax_2) THEN
              w1 = exp((- dsq)/xkappa_1)
              w2 = exp((- dsq)/xkappa_2)
          ELSE IF (dsq<=rmax_1) THEN
              w1 = exp((- dsq)/xkappa_1)
              w2 = anum_2/dsq;
          ELSE  
              !Assume a 1/r**2 weight.
              w1 = anum_1/dsq
              !With anum_2/dsq.
              w2 = gamma * w1
          END IF

          wtot1 = wtot1 + w1
          wtot2 = wtot2 + w2
          ftot1 = ftot1 + w1 * activeNetwork % obs (k) % value
          ftot2 = ftot2 + w2 * dvar (k)
 
      END DO stations_loop
!
!			if (wtot1==0.0 || wtot2==0.0) printf("wts total zero\n");
!
			grid % mat (i,j) =  ftot1/wtot1 + ftot2/wtot2

  END IF
  END DO col_loop
END DO row_loop



DEALLOCATE(activeNetwork % obs)
DEALLOCATE (dvar)


RETURN
END SUBROUTINE BarnesOI


!============================================================
!| Description:
! Inverse matrix
! Method: Based on Doolittle LU factorization for Ax=b
!-----------------------------------------------------------
! input:
! gamma(nest,nest) - array of covariance coefficients for matrix gamma
! nest      - dimension
! outpu:
! inv_gamma(nest,nest) - inverse matrix of A
! 
! IMPORTANT: the original matrix gamma(nest,nest) is destroyed 
! during the calculation. Define A for avoiding this.
!
 SUBROUTINE inverse(gamma,nest,inv_gamma)

  IMPLICIT NONE 
! Calling variables  
  INTEGER, INTENT(IN) :: nest
  DOUBLE PRECISION, DIMENSION(nest,nest), INTENT(IN) :: gamma 
  DOUBLE PRECISION, DIMENSION(nest,nest), INTENT(OUT) :: inv_gamma
! Local variables  
  DOUBLE PRECISION, DIMENSION(nest,nest)  :: A,L,U
  DOUBLE PRECISION, DIMENSION(nest)  ::  b,d,x
  DOUBLE PRECISION ::  coeff
  INTEGER :: i, j, k
  
! step 0: initialization for matrices 
  A=gamma
  inv_gamma=0.
  L=0.
  U=0.
  b=0.
  d=0.
  x=0.
  coeff=0.
  

! step 1: forward elimination
  DO k=1, nest-1
     DO i=k+1,nest
        coeff=A(i,k)/A(k,k)
        
        L(i,k) = coeff
        DO j=k+1,nest
           A(i,j) = A(i,j)-coeff*A(k,j)
        ENDDO
     ENDDO
  ENDDO
  

! Step 2: prepare L and U matrices 
! L matrix is A matrix of the elimination coefficient
! + the diagonal elements are 1.0
  DO i=1,nest
     L(i,i) = 1.
  ENDDO
! U matrix is the upper triangular part of A
  DO j=1,nest
     DO i=1,j
        U(i,j) = A(i,j)
     ENDDO
  ENDDO

! Step 3: compute columns of the inverse matrix C
  DO k=1,nest
     b(k)=1.0
     d(1) = b(1)
! Step 3a: Solve Ld=b using the forward substitution
     DO i=2,nest
        d(i)=b(i)
        DO j=1,i-1
           d(i) = d(i) - L(i,j)*d(j)
        ENDDO
     ENDDO
! Step 3b: Solve Ux=d using the back substitution
     x(nest)=d(nest)/U(nest,nest)
     DO i = nest-1,1,-1
        x(i) = d(i)
        DO j=nest,i+1,-1
           x(i)=x(i)-U(i,j)*x(j)
        ENDDO
        x(i) = x(i)/u(i,i)
     ENDDO
! Step 3c: fill the solutions x(nest) into column k of C
     DO i=1,nest
        inv_gamma(i,k) = x(i)
     ENDDO
     b(k)=0.0
  ENDDO
 END SUBROUTINE inverse

 
!==============================================================================
!| Description:
! Sort distances in increasing order
!
!    Sort an array and make the same interchanges in
!            an auxiliary array.  
!
!   Description of Parameters
!      distance - array of values to be sorted   
!      dist_sort - array to be carried with distance (all swaps of distance elements are
!          matched in dist_sort .  After the sort ind_vec contains the original
!         position of the value i in the unsorted distance array and ind_vec_sort the
!         initial indexes once sorted
!      nest - number of stations to be sorted
! 
  SUBROUTINE Sort (distance,ind_vec,nest,dist_sort,ind_vec_sort) 

  IMPLICIT NONE


! Calling variables  
  INTEGER, INTENT(IN) :: nest
  DOUBLE PRECISION, DIMENSION(nest), INTENT(IN) :: distance 
  DOUBLE PRECISION, DIMENSION(nest), INTENT(OUT) :: dist_sort
  INTEGER, DIMENSION(nest), INTENT(IN) :: ind_vec    
  INTEGER, DIMENSION(nest), INTENT(OUT) :: ind_vec_sort  
! Local variables  
  INTEGER :: i, iswap, itemp, iswap1
  DOUBLE PRECISION, DIMENSION(nest)  ::  temp
  
!
   dist_sort=distance
   ind_vec_sort=ind_vec
      
!    MINLOC is a FORTRAN 90 function that returns the index value for the
!    minimum element in the array
   DO  i=1,nest-1

       iswap=MINLOC(dist_sort(i:nest),dim=1)
       
       iswap1=iswap+i-1
       
       IF(iswap1 .NE. i) THEN
          temp(i)=dist_sort(i)
          dist_sort(i)=dist_sort(iswap1)
          dist_sort(iswap1)=temp(i)
          
          itemp=ind_vec_sort(i)
          ind_vec_sort(i)=ind_vec_sort(iswap1)
          ind_vec_sort(iswap1)=itemp
         ENDIF
    ENDDO

  END SUBROUTINE SORT
  
 









END MODULE SpatialInterpolation