!!  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_CREATEIZONE

USE MOD_IDF, ONLY : IDFREAD,IDFWRITE,IDFNULLIFY,IDFOBJ,IDFDEALLOCATE,IDFCOPY
USE MOD_UTL, ONLY : ITOS,UTL_GETUNIT,UTL_CREATEDIR
USE IMODVAR, ONLY : DP_KIND,SP_KIND
USE MOD_QKSORT, ONLY : QKSORT_INT

CONTAINS

 !###======================================================================
 SUBROUTINE CREATEIZONE_MAIN(FIDF,PFOLDER,OFOLDER,TPARAMETER,NLAY,MINF,IZONEOFFSET,IGROUPOFFSET)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: PFOLDER,OFOLDER,TPARAMETER
 INTEGER,INTENT(IN) :: NLAY,IZONEOFFSET,IGROUPOFFSET
 REAL(KIND=DP_KIND),INTENT(IN) :: MINF
 CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: FIDF
 TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: IDF
 CHARACTER(LEN=256) :: LINE
 INTEGER :: I,J,JZ,NZ,TZ,IROW,ICOL,IG,IU,JU
 REAL(KIND=DP_KIND) :: F
 
 CALL UTL_CREATEDIR(OFOLDER)
 IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\PARAM.TXT',STATUS='UNKNOWN',ACTION='WRITE')
 JU=UTL_GETUNIT(); OPEN(JU,FILE=TRIM(OFOLDER)//'\ZONES.TXT',STATUS='UNKNOWN',ACTION='WRITE')
  
 ALLOCATE(IDF(NLAY)); DO I=1,SIZE(IDF); CALL IDFNULLIFY(IDF(I)); ENDDO
  
 IG=IGROUPOFFSET; NZ=IZONEOFFSET 
 DO J=1,SIZE(FIDF)
  IG=IG+1
  TZ=0
  DO I=1,NLAY
   CALL IDFDEALLOCATE(IDF,SIZE(IDF))
   IDF(I)%FNAME=TRIM(PFOLDER)//'\FRACTIONS_L'//TRIM(ITOS(I))//'\'//FIDF(J)
   IF(IDFREAD(IDF(I),IDF(I)%FNAME,1,IQ=1))THEN
    JZ=0
    DO IROW=1,IDF(I)%NROW; DO ICOL=1,IDF(I)%NCOL
     !## take fraction only
     F=MOD(IDF(I)%X(ICOL,IROW),1.0D0)
     IF(F.GE.MINF)THEN
      JZ=JZ+1
      IF(F.GT.0.99)THEN
       IDF(I)%X(ICOL,IROW)=REAL(NZ+1) 
      ELSE
       IDF(I)%X(ICOL,IROW)=REAL(NZ+1)+IDF(I)%X(ICOL,IROW)
      ENDIF
     ELSEIF(IDF(I)%X(ICOL,IROW).EQ.1.0D0)THEN
      JZ=JZ+1; IDF(I)%X(ICOL,IROW)=REAL(NZ+1)
     ELSE
      IDF(I)%X(ICOL,IROW)=0.0D0 
     ENDIF
    ENDDO; ENDDO
    TZ=TZ+JZ
    IF(JZ.GT.0)THEN
     NZ=NZ+1
     IDF(I)%FNAME=TRIM(OFOLDER)//'\FRACTIONS_L'//TRIM(ITOS(I))//'\'//TRIM(FIDF(J))
     IF(.NOT.IDFWRITE(IDF(I),IDF(I)%FNAME,1))RETURN

     WRITE(LINE,'(I3,A4,2(I3,A1),5(F5.2,A1),I3)') 1,','//TRIM(TPARAMETER)//',',I,',',NZ,',',1.0D0,',',1.1,',',0.1,',',10.0D0,',',10.0D0,',',IG
     WRITE(IU,'(A)') TRIM(LINE)
     WRITE(JU,'(A)') TRIM(IDF(I)%FNAME)
    ENDIF
    
   ENDIF
  ENDDO
  IF(TZ.LE.0)IG=IG-1
 ENDDO

 END SUBROUTINE CREATEIZONE_MAIN
 
 !###======================================================================
 SUBROUTINE CREATEIZONE_LITHOS(FIDF,OFOLDER,TPARAMETER,NLAY,IZONEOFFSET,IGROUPOFFSET)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: OFOLDER,TPARAMETER
 INTEGER,INTENT(IN) :: NLAY,IZONEOFFSET,IGROUPOFFSET
 CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: FIDF
 TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: LIDF,ZIDF
 INTEGER :: ILAY,IROW,ICOL,I,J,K,NULIST,IU,IZ,IG
 INTEGER,DIMENSION(:),ALLOCATABLE :: ULIST,PLIST
 INTEGER,DIMENSION(:,:),ALLOCATABLE :: NUZON
 CHARACTER(LEN=15),DIMENSION(:),ALLOCATABLE :: CULIST
 
 ALLOCATE(LIDF(NLAY)); DO ILAY=1,NLAY; CALL IDFNULLIFY(LIDF(ILAY)); IF(.NOT.IDFREAD(LIDF(ILAY),FIDF(ILAY),1))STOP; ENDDO
 ALLOCATE(ZIDF(NLAY)); DO ILAY=1,NLAY; CALL IDFNULLIFY(ZIDF(ILAY)); CALL IDFCOPY(LIDF(ILAY),ZIDF(ILAY)); ENDDO

 !## read lithos and determine unique if not specified
 I=MAXVAL(LIDF%NCOL)*MAXVAL(LIDF%NROW); ALLOCATE(ULIST(I)); ULIST=0; NULIST=0
 DO ILAY=1,NLAY; DO IROW=1,LIDF(ILAY)%NROW; DO ICOL=1,LIDF(ILAY)%NCOL
  IF(LIDF(ILAY)%X(ICOL,IROW).NE.LIDF(ILAY)%NODATA)THEN
   I=INT(LIDF(ILAY)%X(ICOL,IROW))
   LIDF(ILAY)%X(ICOL,IROW)=REAL(I,8)
   DO J=1,NULIST; IF(I.EQ.ULIST(J))EXIT; ENDDO
   IF(J.GT.NULIST)THEN; NULIST=NULIST+1; ULIST(NULIST)=I; ENDIF
  ENDIF
 ENDDO; ENDDO; ENDDO

 !## sort ulist
 CALL QKSORT_INT(NULIST,ULIST)
 
 !## make a pivottable
 ALLOCATE(PLIST(MAXVAL(ULIST))); PLIST=0
 DO J=1,NULIST
  PLIST(ULIST(J))=J
 ENDDO

 !## echo found unique-values
 ALLOCATE(CULIST(NULIST),NUZON(NULIST,NLAY)); CULIST=''; NUZON=0
 DO ILAY=1,NLAY
  DO IROW=1,LIDF(ILAY)%NROW; DO ICOL=1,LIDF(ILAY)%NCOL
   IF(LIDF(ILAY)%X(ICOL,IROW).NE.LIDF(ILAY)%NODATA)THEN
    I=INT(LIDF(ILAY)%X(ICOL,IROW))
    IZ=PLIST(I)
!    DO J=1,NULIST; IF(I.EQ.ULIST(J))EXIT; ENDDO
    !## get new zone number
!    IZ=J; 
    NUZON(IZ,ILAY)=NUZON(IZ,ILAY)+1 !IZ
   ENDIF
  ENDDO; ENDDO
 ENDDO 
 
 
 WRITE(*,'(/5(A15,1X))') 'LAYER','ZONENUMBER','LITHOCODE','GROUPNUMBER','ACRONYM'
 IZ=IZONEOFFSET
 DO ILAY=1,NLAY
  IG=IGROUPOFFSET
  DO I=1,NULIST
   IG=IG+1
   IZ=IZ+1
   IF(NUZON(I,ILAY).GT.0)THEN
    CULIST(I)=TRIM(TPARAMETER)//'_LITHO'//TRIM(ITOS(ULIST(I)))
    WRITE(*,'(4(I15,1X),A15)') ILAY,IZ,ULIST(I),IG,TRIM(CULIST(I))
   ENDIF
  ENDDO
 ENDDO
 
 CALL UTL_CREATEDIR(OFOLDER)
 IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\IPEST.PRJ',STATUS='UNKNOWN',ACTION='WRITE')
 I=0; DO ILAY=1,NLAY
  DO IZ=1,NULIST
   IF(NUZON(IZ,ILAY).GT.0)I=I+1
  ENDDO
 ENDDO
 WRITE(IU,'(A)') TRIM(ITOS(I))//',(PST) PARAMETER ESTIMATION [-],1'
 WRITE(IU,'(A)') '<no_measure>[<0=interpolation or >0]'
 WRITE(IU,'(A)') '<measure_ipf>,<itype[1-2]>,<ixcol>,<iycol>,<ilcol>,<imcol>,<iwcol>'
 WRITE(IU,'(A)') '<niter>,0.0,0.1,0,0,1.0,0.0,3,<accept.residue>,0.0,2'

 DO ILAY=1,NLAY
  IG=IGROUPOFFSET
  IZ=IZONEOFFSET
  DO I=1,NULIST
   IG=IG+1
   IZ=IZ+1
   IF(NUZON(I,ILAY).GT.0)THEN
    WRITE(IU,'(A)') ' 1,'//TRIM(TPARAMETER)//','//TRIM(ITOS(ILAY))//','// &
                                                  TRIM(ITOS(IZ+(ILAY-1)*NULIST))//',1.0,1.1,0.1,10.0,10.0,'//TRIM(ITOS(IG))//',1,'//TRIM(CULIST(I))
    !## fill in idf with zone number
    DO IROW=1,LIDF(ILAY)%NROW; DO ICOL=1,LIDF(ILAY)%NCOL
     if(icol.eq.552.and.irow.eq.207.and.ilay.eq.27)then
      write(*,*)
     endif
     IF(LIDF(ILAY)%X(ICOL,IROW).NE.LIDF(ILAY)%NODATA)THEN
      K=INT(LIDF(ILAY)%X(ICOL,IROW))
      J=PLIST(K)
      IF(J.EQ.IZ)ZIDF(ILAY)%X(ICOL,IROW)=IZ+(ILAY-1)*NULIST
     ENDIF  
    ENDDO; ENDDO
   ENDIF
  ENDDO
 ENDDO

 WRITE(IU,'(A)') TRIM(ITOS(NLAY))
 
 DO ILAY=1,NLAY
  IF(.NOT.IDFWRITE(ZIDF(ILAY),TRIM(OFOLDER)//'\ZONE_L'//TRIM(ITOS(ILAY))//'.IDF',1))STOP
  WRITE(IU,'(A)') TRIM(OFOLDER)//'\ZONE_L'//TRIM(ITOS(ILAY))//'.IDF'
 ENDDO

 CLOSE(IU); DEALLOCATE(CULIST,ULIST,PLIST)
 CALL IDFDEALLOCATE(ZIDF,SIZE(ZIDF)); CALL IDFDEALLOCATE(LIDF,SIZE(LIDF))

 END SUBROUTINE CREATEIZONE_LITHOS
 
 END MODULE MOD_CREATEIZONE