!!  Copyright (C) Stichting Deltares, 2005-2014.
!!
!!  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.

!###====================================================================
REAL FUNCTION MF2005_TSERIE1HMEAN(X,Y,ICOL,IROW,ILAY)
!###====================================================================
USE IMOD_UTL, ONLY: IMOD_UTL_POL1INTMAIN
USE GLOBAL, ONLY: NROW,NCOL,HNEW,iunit
USE GWFBASMODULE, ONLY: HNOFLO,hdry
use gwfmetmodule, only: cdelr, cdelc
USE GWFLPFMODULE, ONLY: laywet
use M_MF2005_IU, only: IULPF

IMPLICIT NONE
INTEGER,INTENT(IN) :: ICOL,IROW,ILAY
INTEGER :: I,J,IR,IC
REAL,INTENT(IN) :: X,Y
REAL,DIMENSION(3) :: XCRD,YCRD
REAL,DIMENSION(3,3) :: ZCRD
REAL,DIMENSION(1,1) :: XINT

XCRD=0.0; YCRD=0.0; ZCRD=HNOFLO

J=0
DO IR=IROW-1,IROW+1
 J=J+1; I=0
 IF(IR.LT.1)THEN
  YCRD(J)=CDELR(IR)+(CDELC(IR)-CDELC(IR+1))/2.0
 ELSEIF(IR.GT.NROW)THEN
  YCRD(J)=CDELC(IR-1)-(CDELC(IR-2)-CDELC(IR-1))/2.0
 ELSE
  YCRD(J)=(CDELC(IR)+CDELC(IR-1))/2.0
 ENDIF
 DO IC=ICOL-1,ICOL+1
  I=I+1
  IF(IC.LT.1)THEN
   XCRD(I)=CDELR(IC)-(CDELR(IC+1)-CDELR(IC))/2.0
  ELSEIF(IC.GT.NCOL)THEN
   XCRD(I)=CDELR(IC-1)+(CDELR(IC-2)-CDELR(IC-1))/2.0
  ELSE
   XCRD(I)=(CDELR(IC)+CDELR(IC-1))/2.0
  ENDIF
  ZCRD(I,J)=HNEW(MIN(NCOL,MAX(1,IC)),MIN(NROW,MAX(1,IR)),ILAY)
 ENDDO
ENDDO

!## skip dry cells as well - make them nodata
IF(IUNIT(IULPF).GT.0)THEN
 IF(LAYWET(ILAY).NE.0)THEN
  DO I=1,SIZE(ZCRD,1); DO J=1,SIZE(ZCRD,2)
   IF(ZCRD(I,J).EQ.HDRY)ZCRD(I,J)=HNOFLO
  ENDDO; ENDDO
 ENDIF
ENDIF

CALL IMOD_UTL_POL1INTMAIN(1,1,3,3,XCRD,YCRD,ZCRD,(/X,X/),(/Y,Y/),XINT,4,HNOFLO)

MF2005_TSERIE1HMEAN=XINT(1,1)

END FUNCTION