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

 USE MOD_MATH_MERGE_PAR
 USE WINTERACTER
 USE IMODVAR, ONLY : DP_KIND,SP_KIND
 USE MOD_IDF_PAR, ONLY : IDFOBJ
 USE MOD_IDF, ONLY : IDFGETLOC,IDFOPEN,IDFREAD,IDFWRITE,IDFGETVAL,IDFWRITEDIM,IDFNULLIFY,IDFPUTVAL,IDFIROWICOL, &
                     IDFDEALLOCATE 
 USE MOD_UTL, ONLY : UTL_GETUNIT,UTL_IDFSNAPTOGRID,UTL_WAITMESSAGE,UTL_CREATEDIR 
 USE MODPLOT, ONLY : MPW
 TYPE(IDFOBJ),DIMENSION(:),ALLOCATABLE :: MATH

CONTAINS

 !###======================================================================
 LOGICAL FUNCTION MATH1MERGE(IBATCH)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IBATCH
 INTEGER :: IROW,ICOL,IR,IC,I,M,JMASK,IRAT,IRAT1,NIDF
 REAL(KIND=DP_KIND) :: IDFVAL,TW,TD,MD,X,Y,MAXD

 MATH1MERGE=.FALSE.

 NIDF=SIZE(IDFNAMES)

 ALLOCATE(MATH(ABS(IMSK-1):NIDF+1))
 DO I=1,SIZE(MATH); CALL IDFNULLIFY(MATH(I)); ENDDO
 
 !## use/read Mask-IDF
 IF(IMSK.EQ.1)THEN
  IF(.NOT.IDFREAD(MATH(0),MSKNAME,0))THEN
   IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot read/find Mask IDF '//CHAR(13)// &
     TRIM(MSKNAME)//'.','Error')
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'Cannot read Mask IDF '//TRIM(MSKNAME)//'.'
   RETURN
  ENDIF
 ENDIF

 M=0; DO I=1,NIDF !## nidf is output filename
  IF(.NOT.IDFREAD(MATH(I),IDFNAMES(I),0,1))THEN
   IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot read/find IDF '//CHAR(13)// &
     TRIM(IDFNAMES(I))//'.','Error')
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'Cannot read IDF '//TRIM(IDFNAMES(I))//'.'
   MATH(I)%IU=0
  ELSE
   M=M+1
  ENDIF
 ENDDO
 IF(M.LE.0)THEN
  IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'No IDF files available.','Error')
  IF(IBATCH.EQ.1)WRITE(*,'(A)') 'No IDF files available.'
  RETURN
 ENDIF
  
 DO I=1,NIDF
  IF(MATH(I)%IU.GT.0)THEN
   MATH(NIDF+1)%DX  =MATH(I)%DX
   MATH(NIDF+1)%DY  =MATH(I)%DY
   MATH(NIDF+1)%IEQ =MATH(I)%IEQ
   MATH(NIDF+1)%ITB =MATH(I)%ITB
   IF(MATH(NIDF+1)%ITB.EQ.1)THEN
    MATH(NIDF+1)%TOP =MATH(I)%TOP
    MATH(NIDF+1)%BOT =MATH(I)%BOT
   ENDIF
   MATH(NIDF+1)%IXV =0
   MATH(NIDF+1)%NODATA=MATH(I)%NODATA
   EXIT
  ENDIF
 ENDDO
 
 !## it is not allowed to merge non-equidistant IDF files
 IF(MATH(NIDF+1)%IEQ.EQ.1)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You cannot merge non-equidistantial IDF files.','Warning')
  RETURN
 ENDIF
 
 !## zoom for all idfs
 IF(IEXT.EQ.2)THEN

  MATH(NIDF+1)%XMIN= 10.0D10; MATH(NIDF+1)%YMIN= 10.0D10
  MATH(NIDF+1)%XMAX=-10.0D10; MATH(NIDF+1)%YMAX=-10.0D10
  DO I=1,NIDF
   IF(MATH(I)%IU.GT.0)THEN
    MATH(NIDF+1)%XMIN=MIN( MATH(NIDF+1)%XMIN,MATH(I)%XMIN)
    MATH(NIDF+1)%XMAX=MAX( MATH(NIDF+1)%XMAX,MATH(I)%XMAX)
    MATH(NIDF+1)%YMIN=MIN( MATH(NIDF+1)%YMIN,MATH(I)%YMIN)
    MATH(NIDF+1)%YMAX=MAX( MATH(NIDF+1)%YMAX,MATH(I)%YMAX)
    MATH(NIDF+1)%NCOL=INT((MATH(NIDF+1)%XMAX-MATH(NIDF+1)%XMIN)/MATH(NIDF+1)%DX)
    MATH(NIDF+1)%NROW=INT((MATH(NIDF+1)%YMAX-MATH(NIDF+1)%YMIN)/MATH(NIDF+1)%DY)
   ENDIF
  ENDDO
  
 !## current zoom window
 ELSEIF(IEXT.EQ.1)THEN

  IF(MATH(NIDF+1)%IEQ.EQ.1)THEN
   IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot combine IEXT=2 and non-equidistantial IDF"s!','Error')
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'Cannot combine IEXT=2 and non-equidistantial IDF"s!'
   RETURN
  ENDIF
  MATH(NIDF+1)%XMIN=MPW%XMIN; MATH(NIDF+1)%XMAX=MPW%XMAX
  MATH(NIDF+1)%YMIN=MPW%YMIN; MATH(NIDF+1)%YMAX=MPW%YMAX
  CALL UTL_IDFSNAPTOGRID(MATH(NIDF+1)%XMIN,MATH(NIDF+1)%XMAX,MATH(NIDF+1)%YMIN,MATH(NIDF+1)%YMAX, &
                         MATH(NIDF+1)%DX,MATH(NIDF+1)%NCOL,MATH(NIDF+1)%NROW)
 ENDIF

 MATH(NIDF+1)%IU=UTL_GETUNIT()
 CALL UTL_CREATEDIR(OUTNAME(:INDEX(OUTNAME,'\',.TRUE.)-1))
 IF(.NOT.IDFOPEN(MATH(NIDF+1)%IU,OUTNAME,'WO',MATH(NIDF+1)%ITYPE,0,1))RETURN
 
 IF(IBATCH.EQ.1)WRITE(*,*)
 
 MATH(NIDF+1)%DMIN= 10.0D10
 MATH(NIDF+1)%DMAX=-10.0D10

 !## get merged values for constructed idf math(nidf)%idf
 IRAT1=0
 DO IROW=1,MATH(NIDF+1)%NROW
  DO ICOL=1,MATH(NIDF+1)%NCOL

   CALL IDFGETLOC(MATH(NIDF+1),IROW,ICOL,X,Y)

   JMASK=1
   IF(IMSK.EQ.1)THEN
    CALL IDFIROWICOL(MATH(0),IR,IC,X,Y)
    IF(IC.GE.1.AND.IC.LE.MATH(0)%NCOL.AND.IR.GE.1.AND.IR.LE.MATH(0)%NROW)THEN
     JMASK=0
     IF(IDFGETVAL(MATH(0),IR,IC).NE.MATH(0)%NODATA)JMASK=1
    ELSE
     JMASK=0
    ENDIF
   ENDIF
   
   TD=0.0D0
   TW=0.0D0

   IF(JMASK.EQ.1)THEN
   
    MAXD=0.0D0
    DO I=1,NIDF
     IF(MATH(I)%IU.GT.0)THEN
      IF(X.GE.MATH(I)%XMIN.AND.X.LE.MATH(I)%XMAX.AND. &
         Y.GE.MATH(I)%YMIN.AND.Y.LE.MATH(I)%YMAX)THEN
       !## get minimal distance to ANY cross-boundary
       MD=MATH(NIDF+1)%XMAX-MATH(NIDF+1)%XMIN
       IF(MATH(I)%XMIN.NE.MATH(NIDF+1)%XMIN)MD=MIN(MD,X-MATH(I)%XMIN)
       IF(MATH(I)%XMAX.NE.MATH(NIDF+1)%XMAX)MD=MIN(MD,MATH(I)%XMAX-X)
       IF(MATH(I)%YMIN.NE.MATH(NIDF+1)%YMIN)MD=MIN(MD,Y-MATH(I)%YMIN)
       IF(MATH(I)%YMAX.NE.MATH(NIDF+1)%YMAX)MD=MIN(MD,MATH(I)%YMAX-Y)
       CALL IDFIROWICOL(MATH(I),IR,IC,X,Y)
       IDFVAL=IDFGETVAL(MATH(I),IR,IC)
       IF(IDFVAL.NE.MATH(I)%NODATA)THEN
        IF(IINT.EQ.0)THEN
         IF(MD.GT.MAXD)THEN
          MAXD=MD; TW=1.0D0; TD=IDFVAL
         ENDIF
        ELSE
         TW=TW+MD
         TD=TD+IDFVAL*MD
        ENDIF
       ENDIF
      ENDIF
     ENDIF
    ENDDO
   
   ENDIF
   
   IF(TW.NE.0.0D0)THEN
    IDFVAL=TD/TW
    MATH(NIDF+1)%DMIN=MIN(MATH(NIDF+1)%DMIN,IDFVAL)
    MATH(NIDF+1)%DMAX=MAX(MATH(NIDF+1)%DMAX,IDFVAL)
   ELSE
    IDFVAL=MATH(NIDF+1)%NODATA
   ENDIF
   CALL IDFPUTVAL(MATH(NIDF+1),IROW,ICOL,IDFVAL)
  
  ENDDO
  IF(IBATCH.EQ.0)CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,MATH(NIDF+1)%NROW,'Progress ')
  IF(IBATCH.EQ.1)WRITE(6,'(A,F10.2)') '+Progress ',REAL(IROW)/REAL(MATH(NIDF+1)%NROW)*100.0D0
 ENDDO

 IF(.NOT.IDFWRITEDIM(0,MATH(NIDF+1)))RETURN
 CLOSE(MATH(NIDF+1)%IU)

 MATH1MERGE=.TRUE.

 END FUNCTION MATH1MERGE

! 
!
!!###====================================================================
! subroutine pks_imod_utl_idfmerge(fname,ipksdelete)
!!###====================================================================
! ! arguments
! character(len=*), intent(in) :: fname
! integer,intent(in) :: ipksdelete
! ! functions
! logical :: pks7mpimasterwrite
! ! locals
! integer :: iu
! logical :: lpks, lex
! type(idfobj) :: gidf, lidf
! integer :: ios, gncol, gnrow, nrproc, iproc, ic1, ic2, ir1, ir2, ip0, i, j, n, m, iirow, iicol, icol, irow
!! real :: xmin, xmax, ymin, ymax, dx, dy
! character(len=1024) :: s, fp0, f
! real :: nodata
! character(len=1) :: slash 
! REAL(KIND=8),POINTER,DIMENSION(:) :: XCOL=>null(),YROW=>null(),XCOL_BU=>null(),YROW_BU=>null()
! REAL(KIND=8) :: DX1,DX2,DY1,DY2,XC,YC
! 
! !......................................................................
!
! call pks7mpibarrier()
! if (pks7mpimasterwrite()) then
!    call imod_utl_printtext('Start merging PKS output IDF files',0)
!    
!    inquire(file=fname,exist=lex)
!    if(.not.lex) call imod_utl_printtext('File '//trim(fname)//' does not exist',2)
!    
!    iu = getunit(); open(unit=iu,file=fname,status='old')
!    read(iu,*,iostat=ios) nrproc
!    if(ios.ne.0) call imod_utl_printtext('Error reading '//trim(fname),2)
!
!    CALL IDFNULLIFY(GIDF); CALL IDFNULLIFY(LIDF)
!        
!    !## determine common network
!    do
!    
!       read(iu,'(a1024)',iostat=ios) s; if (ios.ne.0) exit
!       if(len_trim(s).eq.0) cycle; if(s(1:1).eq.'#') cycle
!       if(index(s,'"').gt.0)then
!        read(s,*) fp0  != s
!       else
!        fp0=s
!       endif
!       
!!       !## skip all others than heads
!!       if(index(fp0,'HEAD_').le.0)cycle
!       
!       !## start point to insert process-number
!       ip0=len_trim(fp0)-9
!       
!       N=0; M=0; IIROW=0; IICOL=0
!       
!       do iproc = 1, nrproc
!          f = fp0(1:ip0)
!          write(s,'(a,i3.3)') '_p',iproc-1
!          f=trim(f)//trim(s)//'.idf'
!          ! check is file exist
!          inquire(file=f,exist=lex)
!          if (.not.lex)cycle
!
!          ! read the IDF header info only
!          if (.not.idfread(lidf,f,0)) call imod_utl_printtext('Error reading '//trim(f),2); close(lidf%iu)
!
!          N=N+LIDF%NCOL+1
!          IF(ASSOCIATED(XCOL))THEN
!           IF(N.GT.SIZE(XCOL))THEN
!            ALLOCATE(XCOL_BU(N)); DO ICOL=1,SIZE(XCOL); XCOL_BU(ICOL)=XCOL(ICOL); ENDDO; DEALLOCATE(XCOL); XCOL=>XCOL_BU
!           ENDIF
!          ELSE
!           IF(.NOT.ASSOCIATED(XCOL))ALLOCATE(XCOL(N))
!          ENDIF
!          M=M+LIDF%NROW+1
!          IF(ASSOCIATED(YROW))THEN
!           IF(M.GT.SIZE(YROW))THEN
!            ALLOCATE(YROW_BU(M)); DO ICOL=1,SIZE(YROW); YROW_BU(ICOL)=YROW(ICOL); ENDDO; DEALLOCATE(YROW); YROW=>YROW_BU
!           ENDIF
!          ELSE
!           IF(.NOT.ASSOCIATED(YROW))ALLOCATE(YROW(M))
!          ENDIF
!          
!          IF(.NOT.IDFFILLSXSY(LIDF))call imod_utl_printtext('Error allocating sx/sy vectors '//trim(f),2)
!          !## store x and y's
!          DO IROW=0,LIDF%NROW
!           IIROW=IIROW+1; YROW(IIROW)=LIDF%SY(IROW)
!          ENDDO
!          DO ICOL=0,LIDF%NCOL
!           IICOL=IICOL+1; XCOL(IICOL)=LIDF%SX(ICOL)
!          ENDDO
!          CALL IDFDEALLOCATEX(LIDF)
!          
!       end do
!    
!       IF(.NOT.ASSOCIATED(XCOL))THEN
!        STOP 'NOTHING TO MERGE, PROBABLY MERGE-FILE IS EMPTY !!!!'
!       ENDIF
!    
!       !## get network - sort x and y-coordinates
!       CALL IMOD_UTL_QKSORT(N,N,XCOL)
!       CALL IMOD_UTL_QKSORT(M,M,YROW)
!       !## remove duplicates
!       IICOL=1; DO ICOL=2,SIZE(XCOL)
!        IF(XCOL(ICOL).GT.XCOL(IICOL))THEN; IICOL=IICOL+1; XCOL(IICOL)=XCOL(ICOL); ENDIF
!       ENDDO
!       IIROW=1; DO IROW=2,SIZE(YROW)
!        IF(YROW(IROW).GT.YROW(IIROW))THEN; IIROW=IIROW+1; YROW(IIROW)=YROW(IROW); ENDIF
!       ENDDO
!       LIDF%NCOL=IICOL-1; LIDF%NROW=IIROW-1
!       IF(.NOT.IDFALLOCATESXY(LIDF))call imod_utl_printtext('Error allocating new network sx/sy vector',2)
!       DO ICOL=1,IICOL; LIDF%SX(ICOL-1)=XCOL(ICOL); ENDDO; LIDF%XMIN=LIDF%SX(0); LIDF%XMAX=LIDF%SX(LIDF%NCOL)
!       IR1=IIROW+1; DO IROW=1,IIROW
!        IR1=IR1-1; LIDF%SY(IROW-1)=YROW(IR1)
!       ENDDO; LIDF%YMAX=LIDF%SY(0); LIDF%YMIN=LIDF%SY(LIDF%NROW)
!       !## determine whether equi- or non-equi
!       DX1=10.0D10; DX2=0.0D0
!       DO ICOL=1,LIDF%NCOL
!        DX1=MIN(DX1,LIDF%SX(ICOL)-LIDF%SX(ICOL-1))
!        DX2=MAX(DX2,LIDF%SX(ICOL)-LIDF%SX(ICOL-1))
!       ENDDO
!       DY1=10.0D10; DY2=0.0D0
!       DO IROW=1,LIDF%NROW
!        DY1=MIN(DY1,LIDF%SY(IROW-1)-LIDF%SY(IROW))
!        DY2=MAX(DY2,LIDF%SY(IROW-1)-LIDF%SY(IROW))
!       ENDDO
!       !## equidistantial network
!       IF(DX1.EQ.DX2.AND.DY1.EQ.DY2)THEN
!        LIDF%IEQ=0; LIDF%DX=DX1; LIDF%DY=DY1
!       !## non-equidistantial network
!       ELSE
!        LIDF%IEQ=1; LIDF%DX=0.0D0; LIDF%DY=0.0D0
!       ENDIF
!       !## copy settings to the global IDF-file
!       CALL IDFDEALLOCATEX(GIDF); CALL IDFCOPY(LIDF,GIDF)
!       IF(.NOT.IDFALLOCATEX(GIDF))call imod_utl_printtext('Error allocating new network x array',2)
!       
!       i=0; do iproc = 1, nrproc
!          f = fp0(1:ip0)
!          write(s,'(a,i3.3)') '_p',iproc-1
!          f=trim(f)//trim(s)//'.idf'
!       
!          ! check is file exist
!          inquire(file=f,exist=lex)
!          if (.not.lex)cycle
!
!          i=i+1
!          
!          !## read the IDF
!          if (.not.idfread(lidf,f,1))call imod_utl_printtext('Error reading '//trim(f),2)   
!          !## get indices of model to fit in
!          CALL UTL_IDFGETLOC(LIDF,1        ,1        ,XC,YC); CALL UTL_IDFIROWICOL(GIDF,IR1,IC1,XC,YC)
!          CALL UTL_IDFGETLOC(LIDF,LIDF%NROW,LIDF%NCOL,XC,YC); CALL UTL_IDFIROWICOL(GIDF,IR2,IC2,XC,YC)
!          GIDF%NODATA =LIDF%NODATA
!
!          !## set all to nodata for the first time
!          IF(i.eq.1)GIDF%X=GIDF%NODATA
!
!          GIDF%COMMENT=LIDF%COMMENT 
!
!          GIDF%X(IC1:IC2,IR1:IR2) = LIDF%X
!          IF(IPKSDELETE.EQ.1)THEN
!           LIDF%IU=IMOD_UTL_GETUNIT(); OPEN(UNIT=LIDF%IU,FILE=F); CLOSE(LIDF%IU,STATUS='DELETE')
!          ENDIF
!          CALL IDFDEALLOCATEX(LIDF)
!       END DO
!       
!       !## write the global IDF
!       f = fp0(1:ip0)//'.idf'; F=IMOD_UTL_CAPF(F,'U')
!       call imod_utl_printtext('Writing '//trim(f)//'...',0)
!       lidfout=2
!       if (.not.idfwrite(gidf,f,0)) call imod_utl_printtext('Could not write '//trim(f),2)
!
!       IF(ASSOCIATED(XCOL))DEALLOCATE(XCOL)
!       IF(ASSOCIATED(YROW))DEALLOCATE(YROW)
!       
!    end do
!    
!    !## close file
!    call idfdeallocatex(gidf)
!    close(iu) !,status='delete')
!        
!    call imod_utl_printtext('Done merging PKS output IDF files',0)
! end if
! call pks7mpibarrier() 
!
! return
! end subroutine pks_imod_utl_idfmerge 
 
 
 !###======================================================================
 SUBROUTINE MATH1MERGECLOSE(IBATCH) 
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IBATCH 

 IF(IBATCH.EQ.0)THEN; CALL WINDOWSELECT(0); CALL WINDOWOUTSTATUSBAR(4,''); ENDIF

 !## module used outside imod
 IF(ALLOCATED(IDFNAMES))DEALLOCATE(IDFNAMES)
 IF(ALLOCATED(MATH))THEN
  CALL IDFDEALLOCATE(MATH,SIZE(MATH))
  DEALLOCATE(MATH)
 ENDIF

 END SUBROUTINE MATH1MERGECLOSE
 
 END MODULE MOD_MATH_MERGE