Krige Subroutine

SUBROUTINE Krige &
!
(sitedata, anisotropy, nlags, maxlag, neighbors, varModel, grid, gridvar, &
 points, predictions, predictionsVar)


IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: sitedata
INTEGER (KIND = short), INTENT(IN) :: anisotropy !! if = 1, anisotropy is considered 
INTEGER (KIND = short), INTENT(IN) :: nlags !! the number of discrete distances used for comparing data, if 0 it is set to default value
REAL (KIND = float), INTENT(IN) :: maxlag !! maximum distance to consider for semivariogram assessment (m), if 0 it is automatically computed
INTEGER (KIND = short), INTENT(IN) :: neighbors  !! number of neighbors to consider for interpolation
INTEGER (KIND = short), INTENT(IN) :: varModel !!semivariogram model, autofit for atomatic selecting the best among spherical, exponential, and gaussian
TYPE (Coordinate), OPTIONAL, INTENT(IN) :: points (:) !!list of coordinates for prediction (alternative to regular grid)

!Arguments with intent(inout):
TYPE (grid_real), OPTIONAL, INTENT (INOUT) :: grid  !!interpolated grid
TYPE (grid_real), OPTIONAL, INTENT (INOUT) :: gridvar  !!variance of interpolated grid
REAL ( KIND = float), OPTIONAL, DIMENSION (:), INTENT (INOUT) :: predictions  !interpolated values on not regular spaced points
REAL ( KIND = float), OPTIONAL, DIMENSION (:), INTENT (INOUT) :: predictionsVar !variance of interpolated values on not regular spaced points

!Local declarations
TYPE(ObservationalNetwork) :: active_stations
INTEGER (KIND = short) :: count
INTEGER (KIND = short) :: nclosest !!number of closest points to include in interpolation 
INTEGER (KIND = short) :: i,j  
TYPE (site) :: prediction_point
REAL (KIND = float) :: x, y
TYPE(site), DIMENSION (:), ALLOCATABLE ::  estdist
REAL (KIND = float) :: est !estimated value 
REAL (KIND = float) :: variance
INTEGER (KIND = short) :: outmodel !!semivariogram model selected by varfit
!---------------------------end of declarations--------------------------------

  !first check optional input data
  IF (PRESENT (grid) .AND. PRESENT (points) ) THEN
       CALL Catch ('error', 'Kriging interpolation',  &
          'only grid or points should be set for interpolation, they are mutual exclusive' )
  END IF
  
  
  !create a subset of active stations
  CALL ActualObservations (sitedata, count, active_stations)
  
  !compute pair distance
  CALL PairDistance ( active_stations )
  
  !compute direction angle when anisotropy analysis is required
  IF ( anisotropy == 1 ) THEN
      CALL PairDirection (active_stations )
  END IF
  
  !initialize Kriging 
  CALL KrigingInitialize (anisotropy, nlags, maxlag, count, neighbors, nclosest)
  
  !generate empirical semivariogram
  CALL Semivariogram ( active_stations)  
  
  
  !automatic fitting of empirical semivariogram to the best model
  CALL VarFit (varModel, outModel)
  
  
  !Interpolate
  IF (PRESENT (grid) .OR. PRESENT(points) ) THEN
          !populate distance matrix
          CALL DistMatrix (active_stations)
  
          !set coordinate reference system of prediction point 
          IF (PRESENT (grid)) THEN
              prediction_point % xyz % system = grid % grid_mapping
          ELSE
              prediction_point % xyz % system = points(1) % system
          END IF
  
          !allocate estdist
          ALLOCATE ( estdist ( count ) )
  
  
          !Loop through prediction locations
          IF (PRESENT (grid)) THEN  !interpolate regular spaced points
                  DO i = 1, grid % idim
                      DO j = 1, grid % jdim
                          IF ( grid % mat (i,j) /= grid % nodata) THEN
              
                              !set prediction point
                              CALL GetXY (i, j, grid, x, y)
                              prediction_point % xyz % easting = x
                              prediction_point % xyz % northing = y
              
                              !Build Estimation distance Vector  estdist   
		                       CALL Separation (prediction_point, active_stations, estdist)
               
                              !Sort estdist to find the nearest neighbors of estimation location, and return
		                      !the prediction h vector and h matrix.
                              CALL Sorted (estdist, nclosest) 
              
                              !Calculate est covariance from semivariance model, using the estimate to data
		                      !distance vector
                              CALL Covariance(controlpts%h,controlpts%theta,cvest(1:nclosest), sill = partialsill, range = range, varmodel = outModel)
              
                              !Add extra element for lagrangian minimisation
		                      cvest (nclosest + 1 ) = 1.
              
                              !If control points for this datum are the same as for the last, the next steps
		                      !can be skipped over: Since distances between the observations are unchanged
		                      !the observation covariace matrix and its inverse are also unchanged. 
		                      IF (ANY(last .NE. controlpts%oid)) THEN
			                      !Calculate observation covariance matrix from specified semivariance model, 
			                      !using the observation distance matrix (hobs)
                                  CALL Covariance2d(hobs%h, hobs%theta, cvobs(1:nclosest,1:nclosest), sill = partialsill, range = range, varmodel = outModel)
  
			                      !Add extra row and column for lagrangian minimisation
                                  cvobs ( nclosest + 1, : ) = 1.
			                      cvobs ( :, nclosest + 1 ) = 1.
			                      cvobs ( nclosest + 1, nclosest + 1 ) = 0.
  
			                      !Invert Observation covariance Matrix
                                  CALL Matinv ( cvobs, cvobs )
                  
                              END IF
              
                              !Produce estimate for location
                              CALL Interpolate (cvobs, cvest, controlpts%value, est, variance, nclosest, sill = (partialsill + nugget) )
                              grid % mat (i,j) = est
                      
                              IF (PRESENT(gridvar)) THEN
                                 gridvar % mat (i,j) = variance
                              END IF
              
                  
                              !store last set of control points
		                      last = controlpts % oid
              
                          END IF
                      END DO
                  END DO
  
          ELSE !interpolate not regular spaced points
      
              DO i = 1, SIZE(predictions) 
          
                   prediction_point % xyz % easting = points(i) % easting
                   prediction_point % xyz % northing = points(i) % northing
          
                   !Build Estimation distance Vector  estdist   
		            CALL Separation (prediction_point, active_stations, estdist)
               
                    !Sort estdist to find the nearest neighbors of estimation location, and return
		            !the prediction h vector and h matrix.
                    CALL Sorted (estdist, nclosest) 
              
                    !Calculate est covariance from semivariance model, using the estimate to data
		            !distance vector
                    CALL Covariance(controlpts%h,controlpts%theta,cvest(1:nclosest), sill = partialsill, range = range, varmodel = outModel)
              
                    !Add extra element for lagrangian minimisation
		            cvest (nclosest + 1 ) = 1.
              
                    !If control points for this datum are the same as for the last, the next steps
		            !can be skipped over: Since distances between the observations are unchanged
		            !the observation covariace matrix and its inverse are also unchanged. 
		            IF (ANY(last .NE. controlpts%oid)) THEN
			            !Calculate observation covariance matrix from specified semivariance model, 
			            !using the observation distance matrix (hobs).
                        CALL Covariance2d(hobs%h, hobs%theta, cvobs(1:nclosest,1:nclosest), sill = partialsill, range = range, varmodel = outModel)
  
			            !Add extra row and column for lagrangian minimisation
                        cvobs ( nclosest + 1, : ) = 1.
			            cvobs ( :, nclosest + 1 ) = 1.
			            cvobs ( nclosest + 1, nclosest + 1 ) = 0.
  
			            !Invert Observation covariance Matrix
                        CALL Matinv ( cvobs, cvobs )
                  
                    END IF
              
                    !Produce estimate for location
                    !sill = nugget + partial sill = C0 + C1
		            CALL Interpolate (cvobs, cvest, controlpts%value, est, variance, nclosest, sill = (partialsill + nugget) )
            
                    predictions (i) = est
            
                    IF (PRESENT (predictionsVar)) THEN     
                       predictionsVar (i) = variance
                    END IF
              
                  
                    !store last set of control points
		            last = controlpts%oid
          
              END DO
          END IF
  END IF
  
  !free memory
  DEALLOCATE ( active_stations % obs )
  DEALLOCATE ( dist_pairs )
  IF (ALLOCATED (dir_pairs) ) DEALLOCATE ( dir_pairs )
  DEALLOCATE ( dir )
  DEALLOCATE ( lagNumber ) 
  DEALLOCATE ( dist ) 
  DEALLOCATE ( empVariogram ) 
  DEALLOCATE ( modVariogram)
  DEALLOCATE ( pairNumber ) 
  IF (ALLOCATED (obsdist) )  DEALLOCATE ( obsdist )
  DEALLOCATE ( weights )
  DEALLOCATE ( hest )
  DEALLOCATE ( cvest )
  IF (ALLOCATED ( estdist) )   DEALLOCATE ( estdist )
  DEALLOCATE ( controlpts )
  DEALLOCATE ( hobs )
  DEALLOCATE ( cvobs )
  DEALLOCATE ( last )

  RETURN
END SUBROUTINE Krige