!!  Copyright (C) Stichting Deltares, 2005-2023.
!!
!!  This file is part of iMOD.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!  Contact: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.
!!
MODULE MOD_KRIGING

!## simple kriging assumes a mean constant over the entire domain
!## ordinary kriging assumes that the mean is constant in the neighborhoud of the estimated point
USE WINTERACTER
USE RESOURCE
USE MOD_DBL
USE IMODVAR
USE MOD_IDF_PAR, ONLY : IDFOBJ
USE MOD_IDF
USE MOD_UTL
USE MOD_OSD
USE MOD_GRAPH
USE MOD_ASC2IDF_PAR, ONLY : IDFFILE,ASSF_INDICATOR,ELLIPS_IDF,ZONE_IDF,ASSF_IDEPTH,ASSF_TOP,ASSF_BOT,ASSF_ZPLUS, &
   IINT_IDF,INT_IDF
USE MOD_INTERSECT_PAR
USE MOD_INTERSECT
USE MOD_KRIGING_PAR
USE MOD_POLYGON_UTL
USE MOD_LUDCMP, ONLY : LUDCMP_CALC,LUDCMP_INVERSE_MATRIX

CONTAINS

 !###======================================================================
 SUBROUTINE KRIGING_MAIN(MD,XD,YD,ZD,PD,WD,IDF,IBATCH,BO_VALUE,ZONE,IDFV,SEMVAR) 
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(INOUT) :: IDF
 TYPE(IDFOBJ),INTENT(IN) :: ZONE
 TYPE(IDFOBJ),INTENT(INOUT),OPTIONAL :: IDFV
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:),OPTIONAL :: SEMVAR
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: INVSYS,SEMSYS
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:) :: KRANGE
 INTEGER,INTENT(IN) :: IBATCH,MD
 REAL(KIND=DP_KIND),INTENT(IN) :: BO_VALUE
 REAL(KIND=DP_KIND),INTENT(INOUT),DIMENSION(:) :: XD,YD,ZD,PD,WD
 INTEGER :: ND,IROW,ICOL,IRAT,IRAT1,NU
 REAL(KIND=DP_KIND) :: X,Y,Z,MP,KEST,KVAR,R,Z1,Z2
 REAL(KIND=DP_KIND),DIMENSION(3) :: AGL,RAN
 REAL(KIND=DP_KIND),DIMENSION(3,3) :: ROT
 INTEGER :: I,J,MNACTP,N
 LOGICAL :: LT
 INTEGER,ALLOCATABLE,DIMENSION(:) :: UNIQUE,IX
 
 !## clean points for duplicates and get mean for simple kriging
 ROT=0.0D0; RAN=0.0D0
 CALL KRIGING_INIT(MD,XD,YD,ZD,PD,WD,ND,IDF,COINCIDENT,COINCIDENTDIST,IBLANKOUT,BO_VALUE,IBLNTYPE,ROT,RAN)  
! DO I=1,ND
!  WRITE(*,*) XD(I),YD(I),ZD(I)
! ENDDO

 IF(MAXPNT.EQ.0)THEN
  MAXPNT=ND
 ENDIF
 
 IF(IBLANKOUT.EQ.1.OR.ASSOCIATED(IXY))THEN
  IF(PNTSEARCH.EQ.0)MAXPNT=ND; PNTSEARCH=1
 ENDIF

 !## compute mean and substract values from mean - simple kriging
 IF(KTYPE.GT.0)THEN
  MP=0.0D0; DO I=1,ND; MP=MP+PD(I); ENDDO; IF(ND.GT.0.0D0)MP=MP/REAL(ND,8)
  DO I=1,ND; PD(I)=PD(I)-MP; ENDDO
 ENDIF

 !## solve system
 IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'Number of kriging points '//TRIM(VTOS(ND)))
 
 CALL WINDOWOUTSTATUSBAR(1,'Press ESC to terminate')
 
 !## thiessen
 IF(ABS(KTYPE).EQ.6)THEN
  ALLOCATE(INVSYS(1,1)); ALLOCATE(UNIQUE(1),KRANGE(1)); UNIQUE=1; KRANGE=0.0D0
 ELSE
  IF(KTYPE.LT.0)THEN
   !## get number of zones - CHECK THEM ON THE PILOT POINTS
   ALLOCATE(IX(ZONE%NROW*ZONE%NCOL))
   N=0; DO I=1,ND
    CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(I),YD(I))
    IF(IROW.NE.0.AND.ICOL.NE.0)THEN
     N=N+1
     IX(N)=INT(ZONE%X(ICOL,IROW))
    ENDIF
   ENDDO
   CALL UTL_GETUNIQUE_INT(IX,N,NU,0)
   ALLOCATE(UNIQUE(NU)); DO I=1,NU; UNIQUE(I)=IX(I); ENDDO; DEALLOCATE(IX)
  ELSE
   ALLOCATE(UNIQUE(1)); UNIQUE=1
  ENDIF
  AGL=0.0D0; RAN=0.0D0; R=RANGE
  N=ND; IF(KTYPE.LT.0)N=N+1; ALLOCATE(INVSYS(N,N)); INVSYS=0.0D0
  IF(PRESENT(SEMVAR))THEN
   CALL KRIGING_MATRIX(XD,YD,ZD,WD,ND,MAXPNT,R,SILL,NUGGET,KTYPE, &
        IQUADRANT,IDF,IBLANKOUT,BO_VALUE,AGL,RAN,IBLNTYPE,INVSYS,KRANGE,ZONE,UNIQUE,SEMSYS=SEMSYS) 
  ELSE
   CALL KRIGING_MATRIX(XD,YD,ZD,WD,ND,MAXPNT,R,SILL,NUGGET,KTYPE, &
        IQUADRANT,IDF,IBLANKOUT,BO_VALUE,AGL,RAN,IBLNTYPE,INVSYS,KRANGE,ZONE,UNIQUE) 
  ENDIF
 ENDIF

 !IF(RANGE.EQ.0.0)WRITE(*,'(A,F15.3,A)') '>>> range determined',KRANGE(1),' <<<'
 IF(RANGE.EQ.0.0)THEN
  IF(N.EQ.1)THEN
   KRANGE=1000.0D0
   WRITE(*,'(A,F15.3,A)') '>>> RANGE determined - single point so irrelevant RANGE set to ',MAXVAL(KRANGE),' <<<'
  ELSE
   WRITE(*,'(A,F15.3,A)') '>>> RANGE determined',MAXVAL(KRANGE),' <<<'
  ENDIF
 ENDIF

    !!## isotropic
    !IF(ASSOCIATED(ELLIPS_IDF(1,1)%X))THEN
    ! !## set up rotation matrix backwards
    ! AGL=0.0D0; DO I=1,SIZE(ELLIPS_IDF,1); DO J=1,SIZE(ELLIPS_IDF,2)
    !  CALL IDFIROWICOL(ELLIPS_IDF(I,J),JROW,JCOL,X,Y)
    !  IF(JROW.GT.0.AND.JROW.LE.ELLIPS_IDF(I,J)%NROW.AND. &
    !     JCOL.GT.0.AND.JCOL.LE.ELLIPS_IDF(I,J)%NCOL)THEN
    !   IF(ELLIPS_IDF(I,J)%X(JCOL,JROW).NE.ELLIPS_IDF(I,J)%NODATA)THEN
    !    IF(J.EQ.1)THEN
    !     AGL(I)=ELLIPS_IDF(I,J)%X(JCOL,JROW)
    !     !## project reverse on major x/y/z-axis
    !     AGL(I)=-1.0D0*AGL(I)
    !     !## convert to radians
    !     AGL(I)=AGL(I)/(360.0D0/(2.0D0*PI))
    !    ELSE
    !     RAN(I)=ELLIPS_IDF(I,J)%X(JCOL,JROW)
    !    ENDIF
    !   ENDIF
    !  ENDIF
    ! ENDDO; ENDDO
    ! R=MAXVAL(RAN)
    !ELSE 
 
 !## compute semivariance matrix - if desired
 IF(PRESENT(SEMVAR))THEN

  ALLOCATE(SEMVAR(ND,ND))
  DO I=1,ND; DO J=1,ND
   SEMVAR(I,J)=SEMSYS(I,J)
  ENDDO; ENDDO

 ENDIF
 
 IRAT=0; IRAT1=IRAT
 DO IROW=1,IDF%NROW
  MNACTP=0
  DO ICOL=1,IDF%NCOL

   !## skip locations allready filled in
   IF(IDF%X(ICOL,IROW).NE.IDF%NODATA)CYCLE 
   CALL IDFGETLOC(IDF,IROW,ICOL,X,Y)

   LT=.TRUE.
   SELECT CASE (ASSF_IDEPTH)
    !## time
    CASE (0)
     Z1=0.0D0; Z2=Z1
    !## uniform thickness
    CASE (1)
     !## depth in centimeters
     Z1=REAL(ASSF_TOP+ASSF_ZPLUS,8)
     Z2=REAL(ASSF_BOT-ASSF_ZPLUS,8)
    !## spatial thickness
    CASE (2)
     !## depth in centimeters
     Z1=IDFGETXYVAL(INT_IDF(IINT_IDF),X,Y);   IF(Z1.EQ.INT_IDF(IINT_IDF)%NODATA  )LT=.FALSE.
     Z2=IDFGETXYVAL(INT_IDF(IINT_IDF+1),X,Y); IF(Z2.EQ.INT_IDF(IINT_IDF+1)%NODATA)LT=.FALSE.
     !## absolute depth of interface
    CASE (3)
     Z1=REAL(IINT_IDF,8)
     Z2=REAL(IINT_IDF,8)
    !## thickness of interface
    CASE (4)
     Z1=-REAL(IINT_IDF,8)
     Z2=-REAL(IINT_IDF,8)
   END SELECT
  
   IF(.NOT.LT)CYCLE
   
   Z=(Z1+Z2)/2.0D0
   
    !!## isotropic
    !IF(ASSOCIATED(ELLIPS_IDF(1,1)%X))THEN
    ! !## set up rotation matrix backwards
    ! AGL=0.0D0; DO I=1,SIZE(ELLIPS_IDF,1); DO J=1,SIZE(ELLIPS_IDF,2)
    !  CALL IDFIROWICOL(ELLIPS_IDF(I,J),JROW,JCOL,X,Y)
    !  IF(JROW.GT.0.AND.JROW.LE.ELLIPS_IDF(I,J)%NROW.AND. &
    !     JCOL.GT.0.AND.JCOL.LE.ELLIPS_IDF(I,J)%NCOL)THEN
    !   IF(ELLIPS_IDF(I,J)%X(JCOL,JROW).NE.ELLIPS_IDF(I,J)%NODATA)THEN
    !    IF(J.EQ.1)THEN
    !     AGL(I)=ELLIPS_IDF(I,J)%X(JCOL,JROW)
    !     !## project reverse on major x/y/z-axis
    !     AGL(I)=-1.0D0*AGL(I)
    !     !## convert to radians
    !     AGL(I)=AGL(I)/(360.0D0/(2.0D0*PI))
    !    ELSE
    !     RAN(I)=ELLIPS_IDF(I,J)%X(JCOL,JROW)
    !    ENDIF
    !   ENDIF
    !  ENDIF
    ! ENDDO; ENDDO
    ! R=MAXVAL(RAN)
    !ELSE
    ! R=RANGE; AGL=0.0D0; RAN=0.0D0
    !ENDIF
     
   CALL KRIGING_VALUES(X,Y,Z,XD,YD,ZD,PD,WD,ND,KRANGE,SILL,NUGGET,KTYPE,MP,KEST,KVAR, &
          IDF,IBLANKOUT,BO_VALUE,AGL,RAN,IBLNTYPE,ZONE,INVSYS,UNIQUE)
   IDF%X(ICOL,IROW) =KEST
   !## standard deviation (L)
   IF(PRESENT(IDFV))THEN
!    IDFV%X(ICOL,IROW)=KVAR
    IDFV%X(ICOL,IROW)=IDFV%NODATA; IF(KVAR.GE.0.0D0)IDFV%X(ICOL,IROW)=KVAR !SQRT(KVAR)
   ENDIF
  ENDDO
  
  IF(IBATCH.EQ.0)THEN
   SELECT CASE (ABS(KTYPE))
    CASE (1); CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF%NROW,'Progress Ordinary Kriging Linear Semivariogram ('//TRIM(VTOS(ND))//') points')
    CASE (2); CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF%NROW,'Progress Ordinary Kriging Spherical Semivariogram ('//TRIM(VTOS(ND))//') points')
    CASE (3); CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF%NROW,'Progress Ordinary Kriging Exponential Semivariogram ('//TRIM(VTOS(ND))//') points')
    CASE (4); CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF%NROW,'Progress Ordinary Kriging Gaussian Semivariogram ('//TRIM(VTOS(ND))//') points')
    CASE (5); CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF%NROW,'Progress Ordinary Kriging Power Semivariogram ('//TRIM(VTOS(ND))//') points')
    CASE (6); CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF%NROW,'Progress Ordinary Kriging Thiessen Interpolation ('//TRIM(VTOS(ND))//') points')
   END SELECT
  ELSE
   WRITE(6,'(A,F10.3,A)') '+Progress ',REAL(100*IROW)/REAL(IDF%NROW),' % ('//TRIM(VTOS(ND))//') points           '
  ENDIF
 ENDDO
 
 IF(ALLOCATED(XLAG))DEALLOCATE(XLAG); IF(ALLOCATED(ZLAG))DEALLOCATE(ZLAG)
 IF(ASSOCIATED(XY))DEALLOCATE(XY);    IF(ASSOCIATED(IXY))DEALLOCATE(IXY)
 IF(ASSOCIATED(FFCT))DEALLOCATE(FFCT)
 IF(ASSOCIATED(INVSYS)) DEALLOCATE(INVSYS)
 IF(ASSOCIATED(SEMSYS)) DEALLOCATE(SEMSYS)
 IF(ASSOCIATED(KRANGE)) DEALLOCATE(KRANGE)
 IF(ALLOCATED(UNIQUE))DEALLOCATE(UNIQUE)
 
 END SUBROUTINE KRIGING_MAIN

 !###======================================================================
 SUBROUTINE KRIGING_READGEN(ILAY,NGEN,GENFNAME,FGEN,GENILAY)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ILAY,NGEN
 CHARACTER(LEN=*),INTENT(IN),DIMENSION(NGEN) :: GENFNAME
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(NGEN) :: FGEN
 INTEGER,INTENT(IN),DIMENSION(NGEN) :: GENILAY
 REAL(KIND=DP_KIND) :: XMIN,YMIN,XMAX,YMAX
 INTEGER :: IU,IOS,MAXP,MAXN,I,J,K,N,M,IGEN,IFORMAT,MAXPOL,MAXCOL,STRLEN,ITYPE,NP
 CHARACTER(LEN=12) :: CID
 CHARACTER(LEN=11),ALLOCATABLE,DIMENSION(:) :: COLNAMES
 INTEGER,ALLOCATABLE,DIMENSION(:) :: COLWIDTH
 TYPE LBLOBJ
  CHARACTER(LEN=:),ALLOCATABLE :: STRING
 END TYPE LBLOBJ
 TYPE(LBLOBJ),DIMENSION(:),ALLOCATABLE :: LBL
 
 IF(NGEN.LE.0)RETURN
 
 !## store coordinaten per line i nxy and ixy
 MAXP=0; MAXN=0

 !## read polygon in GENfile is available
 N=0; M=0
 DO IGEN=1,NGEN

  !## include proper layers only
  IF(GENILAY(IGEN).NE.ILAY)CYCLE
  IF(LEN_TRIM(GENFNAME(IGEN)).EQ.0)CYCLE
  IF(.NOT.POLYGON_UTL_OPENGEN(GENFNAME(IGEN),IFORMAT,IU))CYCLE
  
  IF(IFORMAT.EQ.0)THEN
  
   !## read overall extent of polygons
   READ(IU) XMIN,YMIN,XMAX,YMAX
   READ(IU) MAXPOL,MAXCOL
   !## labels available
   IF(MAXCOL.GT.0)THEN
    ALLOCATE(COLWIDTH(MAXCOL),LBL(MAXCOL),COLNAMES(MAXCOL))
    READ(IU) (COLWIDTH(I),I=1,MAXCOL)
    DO I=1,MAXCOL; STRLEN=COLWIDTH(I); ALLOCATE(CHARACTER(LEN=STRLEN) :: LBL(I)%STRING); ENDDO
    READ(IU) (COLNAMES(I),I=1,MAXCOL)
   ENDIF

   N=N+MAXPOL

   !## increase memory
   IF(N.GT.MAXN)THEN
    ALLOCATE(IXYDUMMY(N)); DO I=1,MAXN; IXYDUMMY(I)=IXY(I); ENDDO
    IF(ASSOCIATED(IXY))DEALLOCATE(IXY); IXY=>IXYDUMMY
    ALLOCATE(FFCTDUMMY(N)); DO I=1,MAXN; FFCTDUMMY(I)=FFCT(I); ENDDO
    IF(ASSOCIATED(FFCT))DEALLOCATE(FFCT); FFCT=>FFCTDUMMY; MAXN=N
   ENDIF

   N=N-MAXPOL
   
   DO J=1,MAXPOL
    READ(IU) NP,ITYPE
    !## read labels
    IF(MAXCOL.GT.0)READ(IU) (LBL(I)%STRING,I=1,MAXCOL)

    M=M+NP
    
    !## increase memory
    IF(M.GT.MAXP)THEN
     ALLOCATE(XYDUMMY(M,2)); DO I=1,MAXP; DO K=1,2; XYDUMMY(I,K)=XY(I,K); ENDDO; ENDDO
     IF(ASSOCIATED(XY))DEALLOCATE(XY); XY=>XYDUMMY; MAXP=M
    ENDIF

    READ(IU) XMIN,YMIN,XMAX,YMAX
    READ(IU) (XY(I,1),XY(I,2),I=(M-NP)+1,M)

    N=N+1
    IXY(N)=M
    FFCT(N)=FGEN(IGEN)
    
   ENDDO
  
   IF(ALLOCATED(COLWIDTH))DEALLOCATE(COLWIDTH)
   IF(ALLOCATED(COLNAMES))DEALLOCATE(COLNAMES)
   IF(ALLOCATED(LBL))THEN; DO I=1,MAXCOL; DEALLOCATE(LBL(I)%STRING); ENDDO; DEALLOCATE(LBL); ENDIF
  
  ELSEIF(IFORMAT.EQ.1)THEN

   DO

    !## header
    READ(IU,'(A)',IOSTAT=IOS) CID; IF(IOS.NE.0)EXIT
    IF(TRIM(UTL_CAP(CID,'U')).EQ.'END')EXIT
   
    N=N+1

    !## increase memory
    IF(N.GT.MAXN)THEN
     ALLOCATE(IXYDUMMY(MAXN+10)); DO I=1,MAXN; IXYDUMMY(I)=IXY(I); ENDDO
     IF(ASSOCIATED(IXY))DEALLOCATE(IXY); IXY=>IXYDUMMY
     ALLOCATE(FFCTDUMMY(MAXN+10)); DO I=1,MAXN; FFCTDUMMY(I)=FFCT(I); ENDDO
     IF(ASSOCIATED(FFCT))DEALLOCATE(FFCT); FFCT=>FFCTDUMMY; MAXN=MAXN+10
    ENDIF
   
    DO
     M=M+1

     !## increase memory
     IF(M.GT.MAXP)THEN
      ALLOCATE(XYDUMMY(MAXP+1000,2)); DO I=1,MAXP; DO J=1,2; XYDUMMY(I,J)=XY(I,J); ENDDO; ENDDO
      DEALLOCATE(XY); XY=>XYDUMMY; MAXP=MAXP+1000
     ENDIF

     READ(IU,*,IOSTAT=IOS) XY(M,1),XY(M,2)
     !## read END
     IF(IOS.NE.0)EXIT 

    END DO
    M=M-1; IXY(N)=M; FFCT(N)=FGEN(IGEN)

   END DO
 
   CLOSE(IU)
  ENDIF
  
 ENDDO
 
 IF(M.GT.0)THEN
 
  !## fix vector length to 
  IF(N.LT.MAXN)THEN
   ALLOCATE(IXYDUMMY(N));  DO I=1,N; IXYDUMMY(I)=IXY(I);   ENDDO; DEALLOCATE(IXY);  IXY =>IXYDUMMY
   ALLOCATE(FFCTDUMMY(N)); DO I=1,N; FFCTDUMMY(I)=FFCT(I); ENDDO; DEALLOCATE(FFCT); FFCT=>FFCTDUMMY
  ENDIF
  IF(M.LT.MAXP)THEN
   ALLOCATE(XYDUMMY(M,2)); DO I=1,M; DO J=1,2; XYDUMMY(I,J)=XY(I,J); ENDDO; ENDDO; DEALLOCATE(XY); XY=>XYDUMMY
  ENDIF
 ELSE
  DEALLOCATE(IXY,XY,FFCT)
 ENDIF
 
 END SUBROUTINE KRIGING_READGEN

 !###======================================================================
 SUBROUTINE KRIGING_VARIOGRAM(NDD,XD,YD,ZD,ND,IDF,LAGINTERVAL,LAGDISTANCE,IBATCH,SNAME,IFILE,IFIT)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: IDF
 INTEGER,INTENT(IN) :: NDD,LAGINTERVAL
 INTEGER,INTENT(OUT) :: ND
 INTEGER,INTENT(IN),OPTIONAL :: IBATCH,IFILE,IFIT
 CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: SNAME
 REAL(KIND=DP_KIND),INTENT(INOUT),DIMENSION(NDD) :: XD,YD,ZD
 REAL(KIND=DP_KIND),INTENT(IN) :: LAGDISTANCE
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: DVG,ZVG,NLAG
 REAL(KIND=DP_KIND) :: VAR,MEAN,LAGDIST,DX,DY,MZ,ERROR,FIT
 INTEGER :: I,J,K,IU,NPOP,LAGINT,IRAT,IRAT1,IS,IOS,ITYPE,IEXIT,JFIT
 INTEGER(KIND=8) :: N,NN
 TYPE(WIN_MESSAGE) :: MESSAGE
 CHARACTER(LEN=15),DIMENSION(5) :: STYPE=['LINEAR','SPHERICAL','EXPONENTIAL','GAUSSIAN','POWER']
 
 IS=1
 DO
  !## get row/column data for number of unique data sets
  J=0; DO I=1,NDD,IS
   IF(IDF%NCOL.GT.0.AND.IDF%NROW.GT.0)THEN
    IF(XD(I).GE.IDF%XMIN.AND.XD(I).LE.IDF%XMAX.OR. &
       YD(I).GE.IDF%YMIN.AND.YD(I).LE.IDF%YMAX)J=J+1
   ELSE
    J=J+1
   ENDIF
  ENDDO
  NN=J

  !## n.pairs is n(n-1)/2
  N=(NN*(NN-1))/2
  ALLOCATE(DVG(N),STAT=IOS)
  IF(IOS.EQ.0)ALLOCATE(ZVG(N),STAT=IOS)
  IF(IOS.EQ.0)EXIT

  !## try again with less points
  IF(ALLOCATED(DVG))DEALLOCATE(DVG); IF(ALLOCATED(ZVG))DEALLOCATE(ZVG); IS=IS+1
 ENDDO
 
 ND=INT(NN,4)
 
 !## get row/column data for number of unique data sets
 J=0; DO I=1,NDD,IS
  IF(IDF%NCOL.GT.0.AND.IDF%NROW.GT.0)THEN
   IF(XD(I).GE.IDF%XMIN.AND.XD(I).LE.IDF%XMAX.OR. &
      YD(I).GE.IDF%YMIN.AND.YD(I).LE.IDF%YMAX)THEN
    J=J+1; XD(J)=XD(I); YD(J)=YD(I); ZD(J)=ZD(I)
   ENDIF 
  ELSE
   J=J+1; XD(J)=XD(I); YD(J)=YD(I); ZD(J)=ZD(I)
  ENDIF
 ENDDO
 ND=J

 !## compute mean and substract values from mean
 MZ=0.0D0; DO I=1,ND; MZ=MZ+ZD(I); ENDDO; MZ=MZ/REAL(ND)
 DO I=1,ND; ZD(I)=ZD(I)-MZ; ENDDO
  
 DVG=0.0D0; ZVG=0.0D0; K=0
 DO I=1,ND,IS
  DO J=I+1,ND,IS
   DX=(XD(I)-XD(J))**2.0D0; DY=(YD(I)-YD(J))**2.0D0
   IF(DX+DY.GT.0.0D0)THEN
    K=K+1
    DVG(K)=SQRT(DX+DY)
    ZVG(K)=(ZD(I)-ZD(J))**2.0D0
   ENDIF
  ENDDO
 ENDDO
 
 !## get statistics in data
 CALL UTL_STDEF(DVG,K,-999.99D0,VAR,MEAN,NPOP)
 
 !## laginterval=number of intervals
 IF(LAGINTERVAL.NE.0.0D0)THEN
  LAGINT =LAGINTERVAL
 ELSE
  LAGINT=20
 ENDIF

 !## lagdistance=size of interval
 IF(LAGDISTANCE.NE.0.0D0)THEN 
  LAGDIST=LAGDISTANCE
 ELSE
  LAGDIST=(MEAN+(1.0D0*VAR))/LAGINT
 ENDIF

 IF(ALLOCATED(XLAG))DEALLOCATE(XLAG)
 IF(ALLOCATED(ZLAG))DEALLOCATE(ZLAG)
  
 ALLOCATE(ZLAG(0:LAGINT),XLAG(0:LAGINT),NLAG(LAGINT)) 
 XLAG(0)=0.0D0; ZLAG=0.0D0; NLAG=0.0D0
 
 !## compute semivariogram
 IRAT=0; IRAT1=IRAT
 DO I=1,LAGINT
  XLAG(I)=XLAG(I-1)+LAGDIST
  DO J=1,N 
   IF(DVG(J).GE.XLAG(I-1).AND.DVG(J).LT.XLAG(I))THEN
    NLAG(I)=NLAG(I)+1.0D0
    ZLAG(I)=ZLAG(I)+ZVG(J)
   ENDIF
  ENDDO
  IF(NLAG(I).GT.0.0D0)ZLAG(I)=ZLAG(I)/(2.0D0*NLAG(I))
  CALL UTL_WAITMESSAGE(IRAT,IRAT1,I,LAGINT,'Progress computing Semivariance')
 ENDDO
 
 IF(IBATCH.EQ.0)THEN
   IF(ALLOCATED(GRAPH))CALL GRAPH_DEALLOCATE()
  CALL GRAPH_ALLOCATE(1,1)
  ALLOCATE(GRAPH(1,1)%RX(LAGINT))
  ALLOCATE(GRAPH(1,1)%RY(LAGINT))
  GRAPH(1,1)%RX(1:LAGINT)=XLAG(1:LAGINT)
  GRAPH(1,1)%RY(1:LAGINT)=ZLAG(1:LAGINT)
  GRAPH(1,1)%NP=LAGINT
  GRAPH(1,1)%GTYPE=1 
  GRAPH(1,1)%LEGTXT='Value'
  GRAPH(1,1)%ICLR=WRGB(56,180,176)
  CALL UTL_MESSAGEHANDLE(1) 
  GRAPHDIM(1)%GRAPHNAMES=''; GRAPHDIM(1)%IFIXX=0; GRAPHDIM(1)%IFIXY=0; GRAPHDIM(1)%XTITLE='Distance'; GRAPHDIM(1)%YTITLE='Semivariance (m2)'; GRAPHDIM%LDATE=.FALSE.
  GRAPHDIM(1)%IGROUP=1; GRAPHDIM(1)%TEXTSIZE=5.0D0
  CALL GRAPH_INIT(3)
  DO
   CALL WMESSAGE(ITYPE,MESSAGE)
   CALL GRAPH_MAIN(ITYPE,MESSAGE,IEXIT=IEXIT)
   IF(IEXIT.EQ.1)EXIT
  ENDDO
  IF(ALLOCATED(GRAPH))CALL GRAPH_DEALLOCATE()
 ELSE
  IF(TRIM(SNAME).NE.'')THEN
   IU=UTL_GETUNIT()
   CALL OSD_OPEN(IU,FILE=SNAME(:INDEX(SNAME,'.',.TRUE.)-1)//'variogram.txt',STATUS='REPLACE',ACTION='WRITE')
   WRITE(IU,'(3A15)') 'LAG_DISTANCE','VARIOGRAM','NUMBER'
   DO I=1,LAGINT; WRITE(IU,'(2F15.7,I15)') XLAG(I),ZLAG(I),INT(2*NLAG(I)); ENDDO
   !## fit variogram models
   JFIT=0; IF(PRESENT(IFIT))THEN
    JFIT=IFIT
   ENDIF
   IF(JFIT.EQ.1)THEN
    WRITE(IU,'(/A/)') 'SEMIVARIOGRAM FITTED PARAMETERS'
    WRITE(IU,'(6(A15,1X))') 'TYPE','NUGGET','SILL','RANGE','ERROR','GOODNESS_OF_FIT'
    DO I=1,SIZE(STYPE)
     CALL KRIGING_FIT_VARIOGRAMS(LAGINT,XLAG,ZLAG,I,NUGGET,SILL,RANGE,ERROR,FIT)
     WRITE(IU,'(A15,4(1X,F15.1),1X,F15.7)') ADJUSTR(STYPE(I)),NUGGET,SILL,RANGE,ERROR,FIT
    ENDDO
   ENDIF  
   CLOSE(IU)
  ENDIF
 ENDIF
 
 DEALLOCATE(DVG,ZVG) 
 
 END SUBROUTINE KRIGING_VARIOGRAM

 !###======================================================================
 SUBROUTINE KRIGING_FIT_VARIOGRAMS(NLAG,XLAG,ZLAG,ITYPE,NUGGET,SILL,RANGE,ERROR,FIT)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),DIMENSION(3) :: STOPCRIT
 INTEGER,INTENT(IN) :: NLAG,ITYPE
 REAL(KIND=DP_KIND),INTENT(INOUT),DIMENSION(NLAG) :: XLAG,ZLAG
 REAL(KIND=DP_KIND),INTENT(OUT) :: NUGGET,SILL,RANGE,ERROR,FIT
 REAL(KIND=DP_KIND),DIMENSION(3,3) :: F
 REAL(KIND=DP_KIND) :: ME,Z
 REAL(KIND=DP_KIND),DIMENSION(3) :: F1,DF
 REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: E
 INTEGER :: I,J,K,IE,N
 INTEGER,DIMENSION(7) :: IOPT1=[2,1,3,2,2,2,2]
 INTEGER,DIMENSION(7) :: IOPT2=[2,2,2,1,3,2,2]
 INTEGER,DIMENSION(7) :: IOPT3=[2,2,2,2,2,1,3]
 
 ALLOCATE(E(7)); E=0.0D0
 
 !## remove zlag=0.0
 J=0; DO I=1,NLAG
  IF(ZLAG(I).NE.0.0D0)THEN
   J=J+1
   XLAG(J)=XLAG(I)
   ZLAG(J)=ZLAG(I)
  ENDIF
 ENDDO
 !## set rest to zero
 DO I=J+1,NLAG
  XLAG(I)=0.0D0
  ZLAG(I)=0.0D0
 ENDDO
 N=J
 
 !## initial values
 F1=10.0D0; ME=HUGE(1.0)

 !## initial value
 NUGGET=MAXVAL(ZLAG(1:N))/100.0D0
 F1(1)=NUGGET;            DF(1)=NUGGET/10.0D0
 F1(2)=MAXVAL(ZLAG(1:N)); DF(2)=MAXVAL(ZLAG(1:N))/10.0D0
 F1(3)=MAXVAL(XLAG(1:N)); DF(3)=MAXVAL(XLAG(1:N))/10.0D0
 
 DO I=1,3
  STOPCRIT(I)=F1(I)/1000.0D0
 ENDDO
! STOPCRIT(2)=F1(2)/1000.0D0
! STOPCRIT(3)=MAXVAL(XLAG(1:N))/1000.0D0
 
 WRITE(*,'(A5,2A15)') 'PNT','ERROR','NUGGET','SILL','RANGE'
 
 !## determine three points 
 DOLOOP: DO
  !## all values need to be positive
  DO I=1,3
   F(I,1)=MAX(F1(I)-DF(I),0.01D0)
   F(I,2)=F1(I)
   F(I,3)=F1(I)+DF(I)
  ENDDO
!  do i=1,3
!   write(*,'(3f15.3)') (f(i,j),j=1,3)
!  enddo
  !## evaluate all functions (seven options)
  E=0.0D0
  DO I=1,7
   NUGGET=F(1,IOPT1(I))
   SILL  =F(2,IOPT2(I))
   RANGE =F(3,IOPT3(I))
   DO K=1,N
    SELECT CASE (ITYPE)
     !## linear
     CASE (1); Z=XLAG(K)/RANGE
     !## spherical
     CASE (2); Z=(3.0D0*XLAG(K))/(2.0D0*RANGE)-(0.5D0*(XLAG(K)/RANGE)**3.0D0)
     !## exponential
     CASE (3); Z=1.0D0-10.0D0**(-3.0D0*(XLAG(K)/RANGE))
     !## gaussian
     CASE (4); Z=1.0D0-10.0D0**(-1.0D0*(XLAG(K)**2.0D0)/(RANGE**2.0D0))
     !## power
     CASE (5); Z=(XLAG(K)/RANGE)**0.5D0
    END SELECT
    Z=NUGGET+SILL*Z
    E(I)=E(I)+(ZLAG(K)-Z)**2.0D0
   ENDDO
  ENDDO
  !## see how e evolves - search for lowest point - exclude f()=0.0
  ME=HUGE(1.0); DO I=1,SIZE(E)
   IF(E(I).LT.ME)THEN; IE=I; ME=E(I); ENDIF
  ENDDO
  SELECT CASE (IE)
   CASE (2)
    F1(1)=F(1,1)
   CASE (3)
    F1(1)=F(1,3)
   CASE (4)
    F1(2)=F(2,1)
   CASE (5)
    F1(2)=F(2,3)
   CASE (6)
    F1(3)=F(3,1)
   CASE (7)
    F1(3)=F(3,3)
   CASE (1)
    DF=DF/2.0D0
    DO I=1,3; IF(DF(I).LT.STOPCRIT(I))EXIT DOLOOP; ENDDO
  END SELECT
  WRITE(*,'(I5,G15.7,3F15.1)') IE,ME,(F1(I),I=1,3)
 ENDDO DOLOOP
 
 NUGGET=F1(1)
 SILL  =F1(2)
 RANGE =F1(3)
 ERROR =SQRT(E(IE))
 
 DEALLOCATE(E); ALLOCATE(E(N)); E=0.0D0
 DO K=1,N
  SELECT CASE (ITYPE)
   !## linear
   CASE (1); Z=XLAG(K)/RANGE
   !## spherical
   CASE (2); Z=(3.0D0*XLAG(K))/(2.0D0*RANGE)-(0.5D0*(XLAG(K)/RANGE)**3.0D0)
   !## exponential
   CASE (3); Z=1.0D0-10.0D0**(-3.0D0*(XLAG(K)/RANGE))
   !## gaussian
   CASE (4); Z=1.0D0-10.0D0**(-1.0D0*(XLAG(K)**2.0D0)/(RANGE**2.0D0))
   !## power
   CASE (5); Z=(XLAG(K)/RANGE)**0.5D0
  END SELECT
  E(K)=NUGGET+SILL*Z
 ENDDO
 FIT=UTL_GOODNESS_OF_FIT(ZLAG,E,N)
 DEALLOCATE(E)
 
 END SUBROUTINE KRIGING_FIT_VARIOGRAMS
 
 !###======================================================================
 SUBROUTINE KRIGING_MATRIX(XD,YD,ZD,WD,ND,MAXPNT,RANGE,SILL,NUGGET,KTYPE,IQUADRANT,IDF,IBLANKOUT, &
             BO_VALUE,AGL,RAN,IBLNTYPE,INVSYS,KRANGE,ZONE,UNIQUE,SEMSYS)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: IDF
 TYPE(IDFOBJ),INTENT(IN),OPTIONAL :: ZONE
 INTEGER,INTENT(IN) :: ND,MAXPNT,KTYPE,IQUADRANT,IBLANKOUT,IBLNTYPE
 INTEGER,INTENT(IN),DIMENSION(:) :: UNIQUE
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:),INTENT(INOUT) :: INVSYS
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:),INTENT(INOUT) :: KRANGE
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:),INTENT(INOUT),OPTIONAL :: SEMSYS
 REAL(KIND=DP_KIND),INTENT(IN) :: SILL,NUGGET,RANGE,BO_VALUE
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(ND) :: XD,YD,ZD,WD
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(3) :: RAN,AGL
 INTEGER :: I,J,NP,N,ID,JD,IQ,IROW,ICOL,IZONE,JZONE,KZONE,NUNIQUE
 REAL(KIND=DP_KIND) :: DXY,GAMMA,USERANGE,C0,C1,NODATA
 REAL(KIND=DP_KIND),DIMENSION(3,3) :: ROT
 INTEGER,ALLOCATABLE,DIMENSION(:,:) :: PNTQID
 INTEGER,ALLOCATABLE,DIMENSION(:) :: PNTID
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: PNTD,KSYS
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: NDXY,RDXY
 
 NODATA=IDF%NODATA
 
 ALLOCATE(PNTID(ND),PNTQID(4,MAXPNT),PNTD(4,MAXPNT)); PNTID=0; PNTQID=0; PNTD=0.0D0
 
! !## if 2d interpolation, skip zcd
! IF(RAN(3).EQ.0.0D0)ZCD=Z

 !## set up rotation matrice
 CALL UTL_SETUP_ROTATIONMATRIX(AGL,ROT)
 
 !## collect all if range is zero
 USERANGE=RANGE; IF(USERANGE.EQ.0.0D0)USERANGE=HUGE(1.0)
 
 C0      =NUGGET
 C1      =SILL-NUGGET

 NUNIQUE=SIZE(UNIQUE)
 !## simple kriging (ktype.gt.0) and ordinary kriging (ktype.lt.0)
 N=ND; IF(KTYPE.LT.0)N=ND+NUNIQUE

 IF(ALLOCATED(KSYS))   DEALLOCATE(KSYS);   ALLOCATE(KSYS(N,N))
 IF(ASSOCIATED(KRANGE))DEALLOCATE(KRANGE); ALLOCATE(KRANGE(N)); KRANGE=0.0D0
 
 !## compute range automatically
 IF(RANGE.EQ.0.0D0)THEN !.AND.ABS(KTYPE).NE.6)THEN
  ALLOCATE(NDXY(NUNIQUE),RDXY(NUNIQUE)); NDXY=0.0D0; RDXY=0.0D0
  DO I=1,NUNIQUE
   KZONE=UNIQUE(I)   
   DO JD=1,ND
    IZONE=ZONE%NODATA
    CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(JD),YD(JD))
    IF(IROW.NE.0.AND.ICOL.NE.0)THEN
     IZONE=ZONE%X(ICOL,IROW)
    ENDIF  
    IF(IZONE.NE.KZONE)CYCLE
    !## maximal search is MAXPNT - get them
    DO ID=1,ND
     !## skip current point itself
     IF(ID.EQ.JD)CYCLE
     JZONE=ZONE%NODATA
     CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(ID),YD(ID))
     IF(IROW.NE.0.AND.ICOL.NE.0)THEN
      JZONE=ZONE%X(ICOL,IROW)
     ENDIF
     !## skip this point, not in same zone
     IF(JZONE.NE.KZONE)CYCLE
     !## get distance between points - remove whenever intersected by faults
     DXY=KRIGING_DIST(XD(JD),YD(JD),ZD(JD),XD(ID),YD(ID),ZD(ID),WD(ID),IDF,IBLANKOUT,BO_VALUE,USERANGE,IBLNTYPE,ROT,RAN)
     !## otherside of fault or distance too big, take next, skip this interrelationship
     RDXY(I)=RDXY(I)+DXY
     NDXY(I)=NDXY(I)+1.0D0
    ENDDO
   ENDDO
  ENDDO
  DO I=1,NUNIQUE
   NDXY(I)=NDXY(I)*(2.0D0/3.0D0)
   IF(NDXY(I).GT.0.0D0)RDXY(I)=RDXY(I)/NDXY(I)
  ENDDO
  DEALLOCATE(NDXY)
  DO JD=1,ND
   IZONE=ZONE%NODATA
   CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(JD),YD(JD))
    IF(IROW.NE.0.AND.ICOL.NE.0)THEN
    IZONE=ZONE%X(ICOL,IROW)
   ENDIF  
   !## point outside model
   IF(IZONE.EQ.ZONE%NODATA)CYCLE
   DO I=1,NUNIQUE
    IF(UNIQUE(I).EQ.IZONE)THEN
     KRANGE(JD)=RDXY(I)
     EXIT
    ENDIF
   ENDDO
  ENDDO
 ELSE
  KRANGE=RANGE
 ENDIF
 
 !## fill in A() for all points
 KSYS=0.0D0
 DO JD=1,ND
  
  IZONE=ZONE%NODATA
  CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(JD),YD(JD))
   IF(IROW.NE.0.AND.ICOL.NE.0)THEN
   IZONE=ZONE%X(ICOL,IROW)
  ENDIF  
  
  IF(IZONE.EQ.ZONE%NODATA)THEN
   KSYS(JD,JD)=1.0D0; CYCLE  
  ENDIF
  
  !## maximal search is MAXPNT - get them
  PNTQID=0; PNTD=HUGE(1.0D0); IF(ABS(KTYPE).EQ.6)PNTD=0.0D0
  DO ID=1,ND
  
   !## skip current point itself
   IF(ID.EQ.JD)CYCLE
   
   JZONE=ZONE%NODATA
   CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(ID),YD(ID))
   IF(IROW.NE.0.AND.ICOL.NE.0)THEN
    JZONE=ZONE%X(ICOL,IROW)
   ENDIF
   !## skip this point, not in same zone
   IF(JZONE.NE.IZONE)CYCLE
  
   !## get distance between points - remove whenever intersected by faults
   DXY=KRIGING_DIST(XD(JD),YD(JD),ZD(JD),XD(ID),YD(ID),ZD(ID),WD(ID),IDF,IBLANKOUT,BO_VALUE,USERANGE,IBLNTYPE,ROT,RAN)
!   !## otherside of fault or distance too big, take next, skip this interrelationship
!   IF(DXY.GT.USERANGE)CYCLE

   !## determine quadrant
   IQ=1; IF(IQUADRANT.EQ.1)IQ=KRIGING_QUADRANT(XD(ID),YD(ID),XD(JD),YD(JD))

   !## does not fit in current list
   IF(DXY.GE.PNTD(IQ,MAXPNT))CYCLE
   !## get position in sort
   DO I=1,MAXPNT; IF(DXY.LT.PNTD(IQ,I))EXIT; ENDDO
   !## make place for new least-distance point
   DO J=MAXPNT,I+1,-1; PNTQID(IQ,J)=PNTQID(IQ,J-1); PNTD(IQ,J)=PNTD(IQ,J-1); ENDDO
   !## put new location in nearest list
   PNTQID(IQ,I)=ID; PNTD(IQ,I)=DXY
  
  ENDDO
  
  PNTID=0; NP=0 
  DO IQ=1,1+(IQUADRANT*3)
   DO J=1,MAXPNT
    IF(PNTQID(IQ,J).GT.0)THEN
     NP=NP+1; PNTID(NP)=PNTQID(IQ,J) 
    ENDIF
   ENDDO 
  ENDDO 
 
  !## semivariance
  KSYS(JD,JD)=SILL 
  
  !## fill in selected set of points
  DO J=1,NP
   ID=PNTID(J)
   DXY=KRIGING_DIST(XD(JD),YD(JD),ZD(JD),XD(ID),YD(ID),ZD(ID),WD(ID),IDF,IBLANKOUT,BO_VALUE,KRANGE(JD),IBLNTYPE,ROT,RAN)
   GAMMA=KRIGING_GETGAMMA(DXY,KRANGE(JD),C1,C0,KTYPE)
   !## truncate on zero
   KSYS(JD,ID)=SILL-GAMMA !MAX(0.0D0,SILL-GAMMA)
  ENDDO 
  
  !## ordinary kriging
  IF(KTYPE.LT.0)THEN
   !## find zone
   DO I=1,NUNIQUE
    IF(UNIQUE(I).EQ.IZONE)THEN
     KSYS(JD,I+ND)=1.0D0
     KSYS(I+ND,JD)=1.0D0
    ENDIF
   ENDDO
  ENDIF

 ENDDO
 
! open(99,file='d:\tmp_imod.txt',status='replace',action='write')
! WRITE(99,*) 'KSYS'
! do i=1,n
!  write(99,'(999f10.2)') krange(i) !,(ksys(j,i),j=1,n)
! enddo
! flush(99)
 
 !## weigth from semivariogram
 IF(PRESENT(SEMSYS))THEN
  ALLOCATE(SEMSYS(N,N)); SEMSYS=KSYS
 ENDIF
 
 ALLOCATE(INVSYS(N,N)) 
 
 !## if gaussian, use eigenvalue decomposition as matrix is not positive definite
 IF(ABS(KTYPE).EQ.4)THEN
  CALL LUDCMP_INVERSE_MATRIX(N,KSYS,INVSYS)
 ELSE 
  CALL LUDCMP_CALC(N,KSYS,AI=INVSYS)
 ENDIF

! WRITE(99,*) 'INVSYS'
! do i=1,n
!  write(99,'(999f10.2)') (INVSYS(j,i),j=1,n)
! enddo
! close(99)

 DEALLOCATE(PNTID,PNTQID,PNTD)
 IF(ALLOCATED(RDXY))DEALLOCATE(RDXY)

 END SUBROUTINE KRIGING_MATRIX

 !###======================================================================
 SUBROUTINE KRIGING_VALUES(X,Y,Z,XD,YD,ZD,PD,WD,ND,KRANGE,SILL,NUGGET, &
          KTYPE,MP,KEST,KVAR,IDF,IBLANKOUT,BO_VALUE,AGL,RAN,IBLNTYPE,ZONE,INVSYS,UNIQUE)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: IDF
 TYPE(IDFOBJ),INTENT(IN) :: ZONE
 INTEGER,INTENT(IN) :: ND,KTYPE,IBLANKOUT,IBLNTYPE
 INTEGER,INTENT(IN),DIMENSION(:) :: UNIQUE
 REAL(KIND=DP_KIND),INTENT(IN),POINTER,DIMENSION(:,:) :: INVSYS
 REAL(KIND=DP_KIND),INTENT(OUT) :: KEST,KVAR
 REAL(KIND=DP_KIND),INTENT(IN) :: X,Y,Z,SILL,NUGGET,MP,BO_VALUE
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(ND) :: XD,YD,ZD,PD,WD,KRANGE
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(3) :: RAN,AGL
 INTEGER :: I,N,ID,IROW,ICOL,IZONE,JZONE,NUNIQUE,MINI
 REAL(KIND=DP_KIND) :: DXY,GAMMA,C0,C1,MIND
 REAL(KIND=DP_KIND),DIMENSION(3,3) :: ROT
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: B,C
 
 NUNIQUE=SIZE(UNIQUE)
 N=ND
 IF(ABS(KTYPE).NE.6)THEN
  IF(KTYPE.LT.0)N=N+NUNIQUE; ALLOCATE(B(N),C(N)); B=0.0D0
 ENDIF
 
 C0      =NUGGET
 C1      =SILL-NUGGET

 !## set up rotation matrice
 CALL UTL_SETUP_ROTATIONMATRIX(AGL,ROT)

 IZONE=ZONE%NODATA
 CALL IDFIROWICOL(ZONE,IROW,ICOL,X,Y)
 !## include them from outside the model
 IF(IROW.NE.0.AND.ICOL.NE.0)THEN
  IZONE=ZONE%X(ICOL,IROW)
 ENDIF

 DO ID=1,ND

  JZONE=ZONE%NODATA
  CALL IDFIROWICOL(ZONE,IROW,ICOL,XD(ID),YD(ID))
  !## include them from outside the model
  IF(IROW.NE.0.AND.ICOL.NE.0)THEN
   JZONE=ZONE%X(ICOL,IROW)
  ENDIF
  !## skip this point, not in same zone
  IF(JZONE.NE.IZONE)CYCLE
  
  !## get distance between points - remove whenever intersected by faults
  DXY=KRIGING_DIST(X,Y,Z,XD(ID),YD(ID),ZD(ID),WD(ID),IDF,IBLANKOUT,BO_VALUE,KRANGE(ID),IBLNTYPE,ROT,RAN)
  IF(ABS(KTYPE).EQ.6)THEN
   IF(DXY.LT.MIND)THEN
    MIND=DXY
    MINI=ID
   ENDIF
  ELSE
   GAMMA=KRIGING_GETGAMMA(DXY,KRANGE(ID),C1,C0,KTYPE)
   B(ID)=SILL-GAMMA  
  ENDIF
  
 ENDDO
 
 IF(ABS(KTYPE).EQ.6)THEN
  KEST=PD(MINI)
  KVAR=0.0D0
  RETURN
 ENDIF
 
 !## ordinary kriging, fill for lambda
 IF(KTYPE.LT.0)THEN
  DO I=1,SIZE(UNIQUE)
   B(ND+I)=1.0D0
  ENDDO
 ENDIF

 C=MATMUL(INVSYS,B) 
 
 KEST=0.0D0
 DO ID=1,ND 
  IF(B(ID).NE.0.0D0)THEN
   KEST=KEST+C(ID)*PD(ID)
  ENDIF
 ENDDO
 !## global mean (simple kriging)
 IF(KTYPE.GT.0)KEST=KEST+MP
 
 !## standard deviation
 KVAR=0.0D0
 DO ID=1,ND
  IF(B(ID).NE.0.0D0)THEN
!   DXY=KRIGING_DIST(X,Y,Z,XD(ID),YD(ID),ZD(ID),WD(ID),IDF,IBLANKOUT,BO_VALUE,KRANGE(ID),IBLNTYPE,ROT,RAN)
!   GAMMA=KRIGING_GETGAMMA(DXY,KRANGE(ID),C1,C0,KTYPE)
   KVAR=KVAR+C(ID)*B(ID) !(SILL-GAMMA) 
  ENDIF
 ENDDO
 KVAR=SILL-KVAR
 
 DEALLOCATE(B,C)
 
 END SUBROUTINE KRIGING_VALUES

! !###====================================================================
! REAL FUNCTION UTL_KRIGING_RANGE(RANGE,ND,NP,XD,YD,ZD,SELID) 
! !###====================================================================
! IMPLICIT NONE
! INTEGER,INTENT(IN) :: NP,ND
! REAL(KIND=8),INTENT(IN) :: RANGE
! INTEGER,INTENT(IN),DIMENSION(NP) :: SELID
! REAL(KIND=8),INTENT(IN),DIMENSION(ND) :: XD,YD,ZD
! REAL(KIND=8) :: DXY,MDXY
! INTEGER :: ID,JD
! 
! IF(RANGE.GT.0.0D0)THEN
! 
!  UTL_KRIGING_RANGE=RANGE
!
! ELSE
! 
!  !## get largest inter-pilot point distance
!  MDXY=0.0D0; DO ID=1,NP
!   DO JD=1,NP
!    IF(ID.EQ.JD)CYCLE
!    DXY=UTL_DIST_3D(XD(SELID(ID)),YD(SELID(ID)),ZD(SELID(ID)),XD(SELID(JD)),YD(SELID(JD)),ZD(SELID(JD)))
!    MDXY=MAX(DXY,MDXY)
!   ENDDO
!  ENDDO
!  
!  !## increase range with 10%
!  UTL_KRIGING_RANGE=MDXY+0.10D0*MDXY
! 
! ENDIF
! 
!! CALL IMOD_UTL_PRINTTEXT('Kriging Determined Range:'//TRIM(IMOD_UTL_DTOS(MDXY,'F',2))//' meter',1)
!
! END FUNCTION UTL_KRIGING_RANGE
 
 !###======================================================================
 INTEGER FUNCTION KRIGING_QUADRANT(XD,YD,X,Y)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(IN) :: XD,YD,X,Y
 
 IF(XD.GT.X)THEN
  KRIGING_QUADRANT=1; IF(YD.LT.Y)KRIGING_QUADRANT=2
 ELSE
  KRIGING_QUADRANT=4; IF(YD.LT.Y)KRIGING_QUADRANT=3
 ENDIF
 
 END FUNCTION KRIGING_QUADRANT

 !###======================================================================
 REAL(KIND=DP_KIND) FUNCTION KRIGING_DIST(X1,Y1,Z1,X0,Y0,Z0,W,IDF,IBLANKOUT,BO_VALUE, &
             MAXDIST,IBLNTYPE,ROT,RAN)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: IDF
 INTEGER,INTENT(IN) :: IBLANKOUT
 INTEGER,INTENT(IN) :: IBLNTYPE
 REAL(KIND=DP_KIND),DIMENSION(3),INTENT(IN) :: RAN
 REAL(KIND=DP_KIND),DIMENSION(3,3),INTENT(IN) :: ROT
 REAL(KIND=DP_KIND),INTENT(IN) :: X1,Y1,Z1,X0,Y0,Z0,BO_VALUE,MAXDIST,W
 REAL(KIND=DP_KIND) :: X3,Y3,X4,Y4,XINTER,YINTER,A,SSX,DX,DY,D,R,X,Y,Z
 INTEGER :: I,J,I1,I2,ISTATUS,N,IROW,ICOL,IPOL1,IPOL2
 REAL(KIND=DP_KIND),DIMENSION(3) :: XYZE
 
 KRIGING_DIST=UTL_DIST_3D(X1,Y1,Z1,X0,Y0,Z0)
 
 !## to far away (include 1.0D0*maxdist on edge - for smooth edges)
 IF(KRIGING_DIST.GE.MAXDIST)RETURN
 IF(KRIGING_DIST.LE.0.0D0)RETURN

 !## see whether nodata-areas are crossed/similar area
 IF(IBLANKOUT.EQ.1)THEN

  !## include in case point can be "seen" - fault method
  IF(IBLNTYPE.EQ.0)THEN

   !## aspect
   DY=(Y0-Y1); DX=(X0-X1); A=ATAN2(DY,DX)
   SSX=IDF%DX
   D=KRIGING_DIST
   CALL IDFIROWICOL(IDF,IROW,ICOL,X1,Y1)
   DX=COS(A)*SSX
   DY=SIN(A)*SSX
   X3=X1; Y3=Y1
   DO
    X3=X3+DX; Y3=Y3+DY
    !## get new grid location
    CALL IDFIROWICOL(IDF,IROW,ICOL,X3,Y3)
    !## outside window
    IF(IROW.NE.0.AND.ICOL.NE.0)THEN
     !## if crossed-idf value is equal nodata increase distance
     IF(IDF%X(ICOL,IROW).EQ.BO_VALUE)THEN 
      KRIGING_DIST=MAXDIST+1.0D0
      EXIT
     ENDIF
    ENDIF
    SSX=SSX+IDF%DX; IF(SSX.GT.D)EXIT
   ENDDO

  !## include in case of a similar area
  ELSEIF(IBLNTYPE.EQ.1)THEN

   CALL IDFIROWICOL(IDF,IROW,ICOL,X0,Y0); X3=IDF%X(ICOL,IROW)
   CALL IDFIROWICOL(IDF,IROW,ICOL,X1,Y1); X4=IDF%X(ICOL,IROW)
   !## not similar area
   IF(X3.EQ.BO_VALUE.OR.X4.EQ.BO_VALUE)KRIGING_DIST=MAXDIST+1.0D0

  ENDIF
 ENDIF
 
 !## including a genfile
 IF(ASSOCIATED(IXY))THEN

  !## process as a fault
  IF(IBLNTYPE.EQ.0)THEN

   DO J=1,SIZE(IXY) 
    !## see whether point intersect a fault
    I1=1; IF(J.GT.1)I1=IXY(J-1)+1; I2=IXY(J)
    DO I=I1+1,I2
     X3=XY(I-1,1); Y3=XY(I-1,2)
     X4=XY(I  ,1); Y4=XY(I  ,2)
     CALL DBL_IGRINTERSECTLINE(X1,Y1,X0,Y0,X3,Y3,X4,Y4,XINTER,YINTER,ISTATUS)
     !## if line intersect, increase distance as a penalty
     IF(ISTATUS.EQ.5)THEN
      IF(FFCT(J).EQ.0.0D0)THEN
       KRIGING_DIST=MAXDIST+1.0D0
      ELSE
       KRIGING_DIST=KRIGING_DIST*FFCT(J)
      ENDIF
      RETURN
     ENDIF
    ENDDO
   ENDDO

  !## process as a polygon
  ELSEIF(IBLNTYPE.EQ.1)THEN
   
   !## what polygon is x0,y0 ?
   DO J=1,SIZE(IXY)
    I1=1; IF(J.GT.1)I1=IXY(J-1); I2=IXY(J); N=I2-I1+1
    IF(DBL_IGRINSIDEPOLYGON(X0,Y0,XY(I1:I2,1),XY(I1:I2,2),N).EQ.1)EXIT
   ENDDO
   IPOL1=J; IF(IPOL1.GT.SIZE(IXY))IPOL1=0
   
   !## what polygon is x1,y1
   DO J=1,SIZE(IXY)
    I1=1; IF(J.GT.1)I1=IXY(J-1); I2=IXY(J); N=I2-I1+1
    IF(DBL_IGRINSIDEPOLYGON(X1,Y1,XY(I1:I2,1),XY(I1:I2,2),N).EQ.1)EXIT
   ENDDO
   IPOL2=J; IF(IPOL2.GT.SIZE(IXY))IPOL2=0

   IF(IPOL1.NE.IPOL2)KRIGING_DIST=MAXDIST+1.0D0
   
  ENDIF
  
 ENDIF
 
 !## in case of an ellips(oid) see if it is inside the current ellipsoid
 IF(.NOT.UTL_CHECK_UNITY(ROT))THEN
  IF(.NOT.UTL_INSIDEELLIPSOID(X0,Y0,Z0,X1,Y1,Z1,ROT,RAN))THEN

  !## outside ellipsoid, ignore point
   KRIGING_DIST=MAXDIST+1.0D0

  ELSE

   !## http://spatial-analyst.net/ILWIS/htm/ilwisapp/anisotropic_kriging_algorithm.htm  
   !## compute coordinates on sphere with max. range
   R=MAXVAL(RAN)
   
   !## backwards rotate points first
   XYZE(1)=X1-X0
   XYZE(2)=Y1-Y0
   XYZE(3)=Z1-Z0
   XYZE=MATMUL(XYZE,ROT)

   X=X0+XYZE(1)
   R=1.0D0; IF(RAN(2).NE.0.0D0)R=RAN(1)/RAN(2); Y=Y0+XYZE(2)*R
   R=1.0D0; IF(RAN(3).NE.0.0D0)R=RAN(1)/RAN(3); Z=Z0+XYZE(3)*R

   KRIGING_DIST=UTL_DIST_3D(X,Y,Z,X0,Y0,Z0)

  ENDIF
 ENDIF
  
 !## change position of point temporarily
 KRIGING_DIST=W*KRIGING_DIST
 
 END FUNCTION KRIGING_DIST

 !###======================================================================
 FUNCTION UTL_CHECK_UNITY(A) RESULT (LUNITY)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(:,:) :: A
 LOGICAL :: LUNITY
 INTEGER :: I,J
 
 !## false unless proven otherwise
 LUNITY=.FALSE.
 
 IF(SIZE(A,1).EQ.SIZE(A,2))THEN
  DO I=1,SIZE(A,1); DO J=1,SIZE(A,2)
   IF(I.EQ.J)THEN
    IF(A(I,J).NE.1.0D0)RETURN
   ELSE
    IF(A(I,J).NE.0.0D0)RETURN
   ENDIF
  ENDDO; ENDDO
 ENDIF
 
 LUNITY=.TRUE.

 END FUNCTION UTL_CHECK_UNITY

 !###======================================================================
 REAL(KIND=DP_KIND) FUNCTION KRIGING_GETGAMMA(DXY,RANGE,C1,C0,KTYPE)  !c1=sill c0=nugget
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: KTYPE
 REAL(KIND=DP_KIND),INTENT(IN) :: DXY,RANGE,C1,C0
 REAL(KIND=DP_KIND) :: H
 
 !## no part of kriging, beyond given range, equal to sill
 SELECT CASE (ABS(KTYPE))
  CASE (1) !## linear
   H=DXY/RANGE
  CASE (2) !## spherical
   H=(3.0D0*DXY)/(2.0D0*RANGE)-(0.5D0*(DXY/RANGE)**3.0D0)
  CASE (3) !## exponential
   H=1.0D0-10.0D0**(-3.0D0*(DXY/RANGE))
  CASE (4) !## gaussian
   H=1.0D0-10.0D0**(-1.0D0*((DXY**2.0D0)/(RANGE**2.0D0)))
  CASE (5) !## power
   H=(DXY/RANGE)**0.5D0
  CASE DEFAULT
   WRITE(*,*) 'UNKNOWN KTYPE',ABS(KTYPE); PAUSE; STOP
 END SELECT

 H=MIN(H,1.0D0)
 
 KRIGING_GETGAMMA=C0+C1*H
 
 END FUNCTION KRIGING_GETGAMMA
 
 !###======================================================================
 SUBROUTINE KRIGING_INIT(MD,XD,YD,ZD,PD,WD,ND,IDF,COINCIDENT,COINCIDENTDIST, &
    IBLANKOUT,BO_VALUE,IBLNTYPE,AGL,RAN) 
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: IDF
 INTEGER,INTENT(IN) :: MD,COINCIDENT,IBLANKOUT,IBLNTYPE 
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(3) :: RAN,AGL
 REAL(KIND=DP_KIND),DIMENSION(3,3) :: ROT 
 REAL(KIND=DP_KIND),INTENT(INOUT),DIMENSION(:) :: XD,YD,ZD,PD,WD
 REAL(KIND=DP_KIND),INTENT(IN) :: COINCIDENTDIST,BO_VALUE
 INTEGER,INTENT(OUT) :: ND
 INTEGER :: I,J,N,NDBL
 REAL(KIND=DP_KIND) :: XYCRIT,XP,YP,ZP,PP,WP
 
 !## no need to snap
 XYCRIT=0.0D0 !*IDF%DX

 !## set up rotation matrice
 CALL UTL_SETUP_ROTATIONMATRIX(AGL,ROT)

 !## apply a user specified critical distance
 IF(COINCIDENT.EQ.1)XYCRIT=MAX(COINCIDENTDIST,XYCRIT)
 
 !## mark doubles with nodata and compute mean for those double points
 DO I=1,MD
  !## done allready
  IF(PD(I).EQ.IDF%NODATA)CYCLE
  N=1; XP=XD(I); YP=YD(I); ZP=ZD(I); PP=PD(I); WP=WD(I)
  DO J=I+1,MD
   !## done allready
   IF(PD(J).EQ.IDF%NODATA)CYCLE
   !## get distance between points - increase distance whenever points are intersected by fault - is only meant to aggregate if too close
   IF(KRIGING_DIST(XD(I),YD(I),ZD(I),XD(J),YD(J),ZD(J),1.0D0,IDF,IBLANKOUT,BO_VALUE,XYCRIT,IBLNTYPE,ROT,RAN).LE.XYCRIT)THEN !T.XYCRIT)THEN
    N=N+1

    !## take maximal indicator value
    IF(ASSF_INDICATOR.EQ.1)THEN
     PP=MAX(PP,PD(J))
    !## take average value
    ELSE
     PP=PP+PD(J)
    ENDIF

    !## add point to existing point and turn location off
    XP=XP+XD(J); YP=YP+YD(J); ZP=ZP+ZD(J); PD(J)=IDF%NODATA; WP=WP+WD(J)
   ENDIF
  ENDDO
  !## determine average spot
  XD(I)=XP/REAL(N,8); YD(I)=YP/REAL(N,8); ZD(I)=ZP/REAL(N,8); WD(I)=WP/REAL(N,8)

  IF(ASSF_INDICATOR.EQ.1)THEN
   PD(I)=PP
  ELSE
   PD(I)=PP/REAL(N,8)
  ENDIF

 END DO
  
 !## eliminate doubles
 J=0; NDBL=0
 DO I=1,MD
  IF(PD(I).NE.IDF%NODATA)THEN
   J=J+1
   IF(J.NE.I)THEN
    XD(J)=XD(I); YD(J)=YD(I); ZD(J)=ZD(I); PD(J)=PD(I); WD(J)=WD(I)
   ENDIF
  ELSE
   NDBL=NDBL+1
  ENDIF
 END DO
 ND=J
 IF(NDBL.GT.0)THEN
  WRITE(*,'(/1X,A,I10,A/)') 'Found ',NDBL,' double points'
 ENDIF
 
 END SUBROUTINE KRIGING_INIT

 !###======================================================================
 SUBROUTINE KRIGING_UNITTEST() 
 !###======================================================================
 !http://www.unigis.ac.at/fernstudien/UNIGIS_professional/traun/spatial_interpolation/kriging4.htm
 IMPLICIT NONE
 INTEGER,PARAMETER :: MD=3
 REAL(KIND=DP_KIND),DIMENSION(MD) :: XD,YD,ZD,PD,WD !,PCD
 REAL(KIND=DP_KIND),DIMENSION(3) :: RAN,AGL
! REAL(KIND=DP_KIND),DIMENSION(3,3) :: ROT
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: INVSVAR
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:) :: KRANGE
 DATA XD/100.0D0, 50.0D0,250.0D0/
 DATA YD/100.0D0,250.0D0,200.0D0/
 DATA ZD/  0.0D0,  0.0D0,  0.0D0/
 DATA PD/5.0D0  , 20.0D0, 32.0D0/
 DATA WD/1.0D0  ,  1.0D0,  1.0D0/
! DATA XD/2.0,3.0,9.0,6.0,5.0/
! DATA YD/2.0,7.0,9.0,5.0,3.0/
! DATA ZD/3.0,4.0,2.0,4.0,6.0/
! DATA XD/825.0 ,2700.0D0,3700.0D0,2300.0D0,825.0 ,500.0D0/
! DATA YD/3700.0D0,4300.0D0,5200.0D0,5700.0D0,5150.0D0,4950.0D0/
! DATA ZD/13.84 ,12.15 ,12.87 ,12.68 ,14.41 ,14.59/
 REAL(KIND=DP_KIND) :: RANGE,NUGGET,SILL,MP,X,Y,Z,KEST,KVAR,LAGDISTANCE,COINCIDENTDIST
 INTEGER :: MAXPNT,KTYPE,IROW,ICOL,I,PNTSEARCH,ND,NR,NC,MX,LAGINTERVAL,COINCIDENT,IQUADRANT !,IZONE
 TYPE(IDFOBJ) :: IDF,IDFV,ZONE
 INTEGER,DIMENSION(1) :: UNIQUE
 
 AGL=0.0D0; RAN=1.0D0; UNIQUE=1

 IDF%DX=1.0D0; IDF%DY=1.0D0
 IDF%XMIN=MINVAL(XD)-5.0D0; IDF%XMAX=MAXVAL(XD)+5.0D0
 IDF%YMIN=MINVAL(YD)-5.0D0; IDF%YMAX=MAXVAL(YD)+5.0D0
 IDF%NCOL=(IDF%XMAX-IDF%XMIN)/IDF%DX; IDF%NROW=(IDF%YMAX-IDF%YMIN)/IDF%DY
 IF(.NOT.IDFALLOCATEX(IDF))RETURN
 CALL IDFCOPY(IDF,IDFV)
 CALL IDFCOPY(IDF,ZONE); ZONE%X=1.0D0; ZONE%NODATA=-999.0D0

! MAXPNT=MD; RANGE=4141.0D0; NUGGET=0.01D0; SILL=0.78; KTYPE=2; PNTSEARCH=0; NR=0; NC=0
! LAGINTERVAL=20; LAGDISTANCE=RANGE/REAL(LAGINTERVAL); COINCIDENT=0; IQUADRANT=0

 MAXPNT=MD; RANGE=400.0D0; NUGGET=4.0; SILL=9.0; KTYPE=-2; PNTSEARCH=0; NR=0; NC=0
 LAGINTERVAL=20; LAGDISTANCE=RANGE/REAL(LAGINTERVAL); COINCIDENT=0; IQUADRANT=0; COINCIDENT=0.0D0
 MX=RANGE/IDF%DX
 
! SILL=SILL-NUGGET
  
! sill=0.78
! sill=10.0D0
! nugget=5.0
  
 CALL KRIGING_INIT(MD,XD,YD,ZD,PD,WD,ND,IDF,COINCIDENT,COINCIDENTDIST,0,0.0D0,0,AGL,RAN) !,1)
 
 !## compute mean and substract values from mean - simple kriging
 IF(KTYPE.GT.0)THEN
  MP=0.0D0; DO I=1,ND; MP=MP+PD(I); ENDDO; MP=MP/REAL(ND)
  DO I=1,ND; PD(I)=PD(I)-MP; ENDDO
 ENDIF

 ND=MD

 X=150.0D0; Y=150.0D0; Z=0.0D0 
 CALL KRIGING_MATRIX(XD,YD,ZD,WD,ND,MAXPNT,RANGE,SILL,NUGGET,KTYPE,IQUADRANT,IDF,0,0.0D0,AGL,RAN, &
                     IBLNTYPE,INVSVAR,KRANGE,ZONE,UNIQUE) 
    
 DO IROW=1,IDF%NROW
  DO ICOL=1,IDF%NCOL
   X=5.0D0; Y=5.0D0; Z=0.0D0
   CALL IDFGETLOC(IDF,IROW,ICOL,X,Y)
   CALL KRIGING_VALUES(X,Y,Z,XD,YD,ZD,PD,WD,ND,KRANGE,SILL,NUGGET,KTYPE,MP,KEST,KVAR, &
           IDF,0,0.0D0,AGL,RAN,IBLNTYPE,ZONE,INVSVAR,UNIQUE)
   IDF%X(ICOL,IROW)=KEST
   IDFV%X(ICOL,IROW)=SQRT(KVAR)
  ENDDO
 ENDDO
 WRITE(*,*) KEST,KVAR,SQRT(KVAR)
 IF(.NOT.IDFWRITE(IDF ,'D:\unittest_kriging.idf',1))RETURN
 IF(.NOT.IDFWRITE(IDFV,'D:\unittest_krigingstdev.idf',1))RETURN
 OPEN(10,FILE='D:\unittest_kriging_point.ipf',STATUS='REPLACE')
 WRITE(10,*) MD; WRITE(10,*) 3; WRITE(10,*) 'X'; WRITE(10,*) 'Y'; WRITE(10,*) 'Z'; WRITE(10,*) '0,TXT'
 DO I=1,MD; WRITE(10,*) XD(I),YD(I),ZD(I); ENDDO
 CLOSE(10)

 STOP
 END SUBROUTINE KRIGING_UNITTEST

 !###====================================================================
 SUBROUTINE KRIGINGSETTINGS(MAXPNT,KTYPE,RANGE,SILL,NUGGET,PNTSEARCH,COINCIDENT,COINCIDENTDIST,IQUADRANT,NGEN,IBREAK)
 !###====================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(INOUT) :: RANGE,SILL,NUGGET,COINCIDENTDIST
 INTEGER,INTENT(OUT),OPTIONAL :: IBREAK
 INTEGER,INTENT(INOUT) :: MAXPNT,KTYPE,PNTSEARCH,COINCIDENT,IQUADRANT
 INTEGER,INTENT(IN) :: NGEN
 TYPE(WIN_MESSAGE) :: MESSAGE
 INTEGER :: ITYPE,I,IERROR
 
 CALL WDIALOGLOAD(ID_DKRIGING,ID_DKRIGING)
 CALL UTL_DIALOGSHOW(-1,-1,0,3)
 
 IF(NGEN.GT.0)THEN; PNTSEARCH=1; CALL WDIALOGFIELDSTATE(IDF_CHECK1,0); ENDIF
 CALL WDIALOGPUTCHECKBOX(IDF_CHECK1,PNTSEARCH)
 CALL WDIALOGPUTCHECKBOX(IDF_CHECK2,COINCIDENT)
 CALL WDIALOGPUTCHECKBOX(IDF_CHECK3,IQUADRANT)

 CALL WDIALOGPUTINTEGER(IDF_INTEGER1,MAXPNT)
 CALL WDIALOGPUTDOUBLE(IDF_REAL4,COINCIDENTDIST,'(F15.3)')  
 CALL WDIALOGPUTDOUBLE(IDF_REAL1,SILL,'(F15.3)')  
 CALL WDIALOGPUTDOUBLE(IDF_REAL2,RANGE,'(F15.3)') 
 CALL WDIALOGPUTDOUBLE(IDF_REAL3,NUGGET,'(F15.3)')  
 !## put semivariogram type
 IF(KTYPE.GT.5)KTYPE=2; IF(KTYPE.LT.-5)KTYPE=-2; IF(KTYPE.EQ.0)KTYPE=2
 CALL WDIALOGPUTOPTION(IDF_MENU1,ABS(KTYPE))
 !## type of kriging ktype>0 simple kriging, ktype<0 ordinary kriging
 IF(KTYPE.GT.0)CALL WDIALOGPUTOPTION(IDF_MENU2,1)
 IF(KTYPE.LT.0)CALL WDIALOGPUTOPTION(IDF_MENU2,2)
 
 IF(PRESENT(IBREAK))IBREAK=0
 
 CALL WDIALOGFIELDOPTIONS(IDF_INTEGER1,EDITFIELDCHANGED,ENABLED)

 !## choice has been made in calling program
 CALL KRIGINGSETTINGSFIELDS()
 DO
  CALL WMESSAGE(ITYPE,MESSAGE)
  SELECT CASE (ITYPE)
   CASE(FIELDCHANGED)
    SELECT CASE (MESSAGE%VALUE2)
     CASE (IDF_INTEGER1,IDF_CHECK1,IDF_CHECK2,IDF_CHECK3,IDF_MENU1,IDF_MENU2)
      CALL KRIGINGSETTINGSFIELDS()
    END SELECT

   CASE(PUSHBUTTON)
    SELECT CASE (MESSAGE%VALUE1)
     CASE (IDOK)
      CALL WDIALOGGETCHECKBOX(IDF_CHECK1,PNTSEARCH)
      CALL WDIALOGGETCHECKBOX(IDF_CHECK2,COINCIDENT)
      CALL WDIALOGGETCHECKBOX(IDF_CHECK3,IQUADRANT) 
      CALL WDIALOGGETINTEGER(IDF_INTEGER1,MAXPNT)
      CALL WDIALOGGETDOUBLE(IDF_REAL1,SILL) 
      CALL WDIALOGGETDOUBLE(IDF_REAL2,RANGE)  
      CALL WDIALOGGETDOUBLE(IDF_REAL3,NUGGET)  
      CALL WDIALOGGETDOUBLE(IDF_REAL4,COINCIDENTDIST)  
      !## get semivariogram type
      CALL WDIALOGGETMENU(IDF_MENU1,KTYPE)
      CALL WDIALOGGETMENU(IDF_MENU2,I)
      IF(I.EQ.2)KTYPE=-1*KTYPE
      IERROR=0
      IF(PNTSEARCH.EQ.1.AND.MAXPNT.LE.0)THEN
       CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Maximum number of points is '//TRIM(VTOS(MAXPNT))//CHAR(13)// &
         'and need to larger than zero','Error'); IERROR=1
      ENDIF
      IF(RANGE.LE.0.0D0)THEN
       CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Range is '//TRIM(VTOS(RANGE,'F',2))//CHAR(13)// &
         'and need to larger than zero','Error'); IERROR=1
      ENDIF
      IF(SILL.LE.0.0D0)THEN
       CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Sill is '//TRIM(VTOS(SILL,'F',2))//CHAR(13)// &
         'and need to larger than zero','Error'); IERROR=1
      ENDIF
      IF(NUGGET.LT.0.0D0)THEN
       CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Nugget is '//TRIM(VTOS(NUGGET,'F',2))//CHAR(13)// &
         'and need to be larger or equal to zero','Error'); IERROR=1
      ENDIF
      IF(SILL-NUGGET.LE.0.0D0)THEN
       CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Sill minus Nugget is '//TRIM(VTOS(SILL-NUGGET,'F',2))//CHAR(13)// &
         'and need to larger than zero','Error'); IERROR=1
      ENDIF
      IF(IERROR.EQ.0)EXIT
     CASE (IDCANCEL)
      IF(PRESENT(IBREAK))IBREAK=1; EXIT 
    END SELECT
  END SELECT

 ENDDO

 CALL WDIALOGUNLOAD()
 
 END SUBROUTINE KRIGINGSETTINGS

 !###====================================================================
 SUBROUTINE KRIGINGSETTINGSFIELDS()
 !###====================================================================
 IMPLICIT NONE
 INTEGER :: I,J
 
 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,I)
 CALL WDIALOGFIELDSTATE(IDF_INTEGER1,I)
 CALL WDIALOGFIELDSTATE(IDF_CHECK3,I)
 
 CALL WDIALOGGETCHECKBOX(IDF_CHECK2,J)
 CALL WDIALOGFIELDSTATE(IDF_REAL4,J)
 CALL WDIALOGFIELDSTATE(IDF_LABEL13,J)
 
 CALL WDIALOGGETCHECKBOX(IDF_CHECK3,I)
 IF(I.EQ.1)THEN
  CALL WDIALOGGETINTEGER(IDF_INTEGER1,J)
  CALL WDIALOGPUTSTRING(IDF_CHECK3,'Apply Quadrants (total points maximal '//TRIM(VTOS(J*4))//')')
 ELSE
  CALL WDIALOGPUTSTRING(IDF_CHECK3,'Apply Quadrants')
 ENDIF
 
 CALL WDIALOGGETMENU(IDF_MENU1,I)
 SELECT CASE (I)
  CASE (1)
    CALL WDIALOGPUTSTRING(IDF_LABEL11,'Linear act linearly. Appropriate for representing properties with higher level of short-range variability.')
  CASE (2)
    CALL WDIALOGPUTSTRING(IDF_LABEL11,'Spherical acts linear in the beginning. '// &
       'Appropriate for representing properties with higher level of short-range variability')
  CASE (3)
    CALL WDIALOGPUTSTRING(IDF_LABEL11,'Exponential acts linear in the beginning but curves rapidly for points further away and approaces the sill asymptotically.')
  CASE (4)
    CALL WDIALOGPUTSTRING(IDF_LABEL11,'Gaussian weights the points nearest least heavy as a "S"-curve and approaces the sill asymptotically. '// &
       'It gives very smooth varying properties.')
  CASE (5)
    CALL WDIALOGPUTSTRING(IDF_LABEL11,'Power models are appropriate for properties exhibiting fractal behaviour.')
 END SELECT

 CALL WDIALOGGETMENU(IDF_MENU2,I)
 SELECT CASE (I)
  CASE (1)
    CALL WDIALOGPUTSTRING(IDF_LABEL12,'Simple Kriging assumes a mean constant over the entire domain. The result tend to fullfill the total mean.')
  CASE (2)
    CALL WDIALOGPUTSTRING(IDF_LABEL12,'Ordinary Kriging assumes that the mean is constant in the neighborhoud of the estimated point. '// &
      'This option gives more smooth result.')
 END SELECT

 END SUBROUTINE KRIGINGSETTINGSFIELDS

END MODULE MOD_KRIGING