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

USE MOD_IDF, ONLY : IDFREAD,IDFWRITE,IDFNULLIFY,IDFOBJ,IDFDEALLOCATE,IDFCOPY,IDFALLOCATEX,IDFREADSCALE
USE MOD_UTL, ONLY : VTOS,UTL_GETUNIT,UTL_CREATEDIR,UTL_GETUNIQUE_DINT,UTL_SUBST
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,SUBDIVIDE,LAYER,GRPN,OUTPUTTYPE)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: PFOLDER,OFOLDER,TPARAMETER
 INTEGER,INTENT(IN) :: NLAY,IZONEOFFSET,IGROUPOFFSET,OUTPUTTYPE
 INTEGER,INTENT(IN),DIMENSION(NLAY) :: LAYER
 REAL(KIND=DP_KIND),INTENT(IN) :: MINF
 CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: FIDF,GRPN
 CHARACTER(LEN=256),DIMENSION(:),ALLOCATABLE :: ZONEFILES
 TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: IDF
 TYPE(IDFOBJ),INTENT(INOUT) :: SUBDIVIDE
 CHARACTER(LEN=256) :: LINE
 CHARACTER(LEN=15) :: GRPNAME
 CHARACTER(LEN=32) :: GRPNAMEFULL
 INTEGER :: I,II,J,JZ,NZ,TZ,IROW,ICOL,IG,IU,JU,NSUB,N,IS,IIS,NF
 REAL(KIND=DP_KIND) :: F
 INTEGER(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: DXINT
 
 CALL UTL_CREATEDIR(OFOLDER)
 IF(OUTPUTTYPE.EQ.0)THEN
  IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\PARAM_'//TRIM(TPARAMETER)//'.TXT',STATUS='REPLACE',ACTION='WRITE')
  JU=UTL_GETUNIT(); OPEN(JU,FILE=TRIM(OFOLDER)//'\ZONES_'//TRIM(TPARAMETER)//'.TXT',STATUS='REPLACE',ACTION='WRITE')
 ELSE
  IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\IPEST_'//TRIM(TPARAMETER)//'.PRJ',STATUS='REPLACE',ACTION='WRITE')
  JU=0
 ENDIF
 
 ALLOCATE(IDF(NLAY)); DO I=1,SIZE(IDF); CALL IDFNULLIFY(IDF(I)); ENDDO
 ALLOCATE(ZONEFILES(NLAY*SIZE(IDF))); ZONEFILES=''
 
 NF=0

FLOOP: DO J=1,SIZE(FIDF)
  DO II=1,NLAY

   I=ABS(LAYER(II))
  
   IF(LAYER(II).GT.0)THEN
    IDF(II)%FNAME=TRIM(PFOLDER)//'\FRACTIONS_L'//TRIM(VTOS(I))//'\'//FIDF(J)
   ELSE
    IDF(II)%FNAME=TRIM(PFOLDER)//'\'//FIDF(1)
   ENDIF
  
   IF(IDFREAD(IDF(II),IDF(II)%FNAME,1,IQ=1))THEN
    CALL IDFCOPY(IDF(II),SUBDIVIDE)
    IF(TRIM(SUBDIVIDE%FNAME).EQ.'')THEN
     IF(.NOT.IDFALLOCATEX(SUBDIVIDE))STOP 'ERROR ALLOCATING SUBDIVIDE'
     SUBDIVIDE%X=1.0D0; NSUB=1; ALLOCATE(DXINT(1)); DXINT(1)=INT(1,8)
    ELSE
     IF(.NOT.IDFREADSCALE(SUBDIVIDE%FNAME,SUBDIVIDE,2,2,0.0D0,0))STOP 'ERROR READING GIVEN SUBDIVIDE'
     N=SUBDIVIDE%NCOL*SUBDIVIDE%NROW
     ALLOCATE(DXINT(N))
     N=0; DO IROW=1,SUBDIVIDE%NROW; DO ICOL=1,SUBDIVIDE%NCOL
      IF(SUBDIVIDE%X(ICOL,IROW).EQ.SUBDIVIDE%NODATA)SUBDIVIDE%X(ICOL,IROW)=0.0D0
      N=N+1; DXINT(N)=INT(SUBDIVIDE%X(ICOL,IROW),8)
     ENDDO; ENDDO
     !## how many subdivisions
     CALL UTL_GETUNIQUE_DINT(DXINT,N,NSUB,0)
    ENDIF
    EXIT FLOOP
   ENDIF
  ENDDO
 ENDDO FLOOP
 
 WRITE(*,'(/A/)') 'Unique number of subdivisions '//TRIM(VTOS(NSUB))
 WRITE(*,'(A10,A15)') 'No.','Number'
 DO I=1,NSUB;  WRITE(*,'(I10,I15)') I,DXINT(I); ENDDO

 IG=IGROUPOFFSET; NZ=IZONEOFFSET 
 !## formation
 DO J=1,SIZE(FIDF)
  WRITE(*,'(A)') 'Processing '//TRIM(FIDF(J))//' ('//TRIM(VTOS(REAL(J,8)/REAL(SIZE(FIDF),8)*100.0D0,'F',2))//'%)           '
  DO IIS=1,NSUB
   
   IS=DXINT(IIS)
   
   IF(TRIM(SUBDIVIDE%FNAME).NE.'')WRITE(*,'(A)') '> processing subdivide '//TRIM(VTOS(IS))//' ...'
   
   IG=IG+1
   TZ=0
   DO II=1,NLAY

    I=ABS(LAYER(II))

    CALL IDFDEALLOCATE(IDF,SIZE(IDF))
    IF(LAYER(II).GT.0)THEN
     IDF(II)%FNAME=TRIM(PFOLDER)//'\FRACTIONS_L'//TRIM(VTOS(I))//'\'//FIDF(J)
    ELSE
     IDF(II)%FNAME=TRIM(PFOLDER)//'\'//FIDF(J)
    ENDIF
 
    IF(IDFREAD(IDF(II),IDF(II)%FNAME,1,IQ=1))THEN
     WRITE(*,'(A)') ' > '//TRIM(IDF(II)%FNAME) //' ...'
     JZ=0
     DO IROW=1,IDF(II)%NROW; DO ICOL=1,IDF(II)%NCOL
      
      IF(LAYER(II).LT.0)IDF(II)%X(ICOL,IROW)=1.0D0

      !## skip if not in current subdivide
      IF(INT(SUBDIVIDE%X(ICOL,IROW)).NE.IS)THEN
       IDF(II)%X(ICOL,IROW)=0.0D0; CYCLE
      ENDIF
      
      !## take fraction only
      F=MOD(IDF(II)%X(ICOL,IROW),1.0D0)
      IF(F.GE.MINF)THEN
       JZ=JZ+1
       IF(F.GT.0.99D0)THEN
        IDF(II)%X(ICOL,IROW)=REAL(NZ+1,8) 
       ELSE
        IDF(II)%X(ICOL,IROW)=REAL(NZ+1,8)+IDF(II)%X(ICOL,IROW)
       ENDIF
      ELSEIF(IDF(II)%X(ICOL,IROW).EQ.1.0D0)THEN
       JZ=JZ+1; IDF(II)%X(ICOL,IROW)=REAL(NZ+1,8)
      ELSE
       IDF(II)%X(ICOL,IROW)=0.0D0 
      ENDIF
     
     ENDDO; ENDDO
     
     TZ=TZ+JZ
     IF(JZ.GT.0)THEN
      NZ=NZ+1
      IF(NSUB.EQ.1)THEN
       IDF(II)%FNAME=TRIM(OFOLDER)//'\FRACTIONS_L'//TRIM(VTOS(I))//'\'//TRIM(FIDF(J))
      ELSE
       IDF(II)%FNAME=TRIM(OFOLDER)//'\FRACTIONS_L'//TRIM(VTOS(I))//'\SUBDIVIDE'//TRIM(VTOS(IS))//'\'//TRIM(FIDF(J))
      ENDIF
      IF(.NOT.IDFWRITE(IDF(II),IDF(II)%FNAME,1))RETURN
      IF(NSUB.GT.1)THEN
       GRPNAMEFULL=TRIM(GRPN(J))//'_S'//TRIM(VTOS(IS))
      ELSE
       GRPNAMEFULL=TRIM(GRPN(J))
      ENDIF
      !## see whether to decrease length of parameter by removing special characters
      DO
       IF(GRPNAMEFULL(16:16).NE.' ')THEN
        IF(INDEX(GRPNAMEFULL,'_').EQ.0)THEN
         WRITE(*,*) 'Groupname too long, iMOD cannot shorten it, change the groupname yourself'; PAUSE; STOP
        ENDIF
        CALL UTL_SUBST(GRPNAMEFULL,'_','')
       ELSE
        EXIT
       ENDIF
      ENDDO
      GRPNAME=GRPNAMEFULL(1:15)
      WRITE(LINE,'(I1,A4,I3,A1,I7,A1,5(F5.2,A1),I7,A3,A15)') 1,','//TRIM(TPARAMETER)//',',I,',',NZ,',',1.0,',',1.1,',',0.1,',',10.0,',',10.0,',',IG,'1,',TRIM(GRPNAME)
      WRITE(IU,'(A)') TRIM(LINE)
      NF=NF+1; ZONEFILES(NF)=IDF(II)%FNAME
      IF(JU.NE.0)WRITE(JU,'(A)') TRIM(IDF(II)%FNAME)
     ENDIF
    
    ENDIF
   ENDDO
   IF(TZ.LE.0)IG=IG-1
  ENDDO
 ENDDO
 
 IF(OUTPUTTYPE.EQ.1)THEN
  WRITE(IU,*) NF
  DO J=1,NF
   WRITE(IU,'(A)') TRIM(ZONEFILES(J))
  ENDDO
 ENDIF
 
 DEALLOCATE(DXINT); CLOSE(IU); IF(JU.NE.0)CLOSE(JU)
 
 END SUBROUTINE CREATEIZONE_MAIN
 
 !###======================================================================
 SUBROUTINE CREATEIZONE_LITHOS(FIDF,ILS,ACRONYM,OFOLDER,TPARAMETER,NLAY,IZONEOFFSET,IGROUPOFFSET,OUTPUTTYPE)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: OFOLDER,TPARAMETER,ACRONYM
 INTEGER,INTENT(IN) :: NLAY,IZONEOFFSET,IGROUPOFFSET,OUTPUTTYPE
 INTEGER,INTENT(IN),DIMENSION(NLAY) :: ILS
 CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: FIDF
 TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: LIDF,ZIDF
 INTEGER :: ILAY,IROW,ICOL,I,J,K,NULIST,IU,JU,IZ,IG
 INTEGER,DIMENSION(:),ALLOCATABLE :: ULIST,PLIST
 INTEGER,DIMENSION(:,:),ALLOCATABLE :: NUZON
 CHARACTER(LEN=15),DIMENSION(:),ALLOCATABLE :: CULIST
 
 CALL UTL_CREATEDIR(OFOLDER)
 IF(OUTPUTTYPE.EQ.0)THEN
  IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\PARAM_'//TRIM(TPARAMETER)//'.TXT',STATUS='REPLACE',ACTION='WRITE')
  JU=UTL_GETUNIT(); OPEN(JU,FILE=TRIM(OFOLDER)//'\ZONES_'//TRIM(TPARAMETER)//'.TXT',STATUS='REPLACE',ACTION='WRITE')
 ELSE
  IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\IPEST_'//TRIM(TPARAMETER)//'.PRJ',STATUS='REPLACE',ACTION='WRITE')
  JU=0
 ENDIF
 
 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)
    !## get new zone number
    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
    IF(LEN_TRIM(ACRONYM).EQ.0)THEN
     CULIST(I)=TRIM(TPARAMETER)//'_LITHO'//TRIM(VTOS(ULIST(I)))
    ELSE
     CULIST(I)=TRIM(TPARAMETER)//'_'//TRIM(ACRONYM)//'_L'//TRIM(VTOS(ULIST(I)))
    ENDIF
    WRITE(*,'(4(I15,1X),A15)') ILAY,IZ,ULIST(I),IG,TRIM(CULIST(I))
   ENDIF
  ENDDO
 ENDDO
 
 I=0; DO ILAY=1,NLAY
  DO IZ=1,NULIST
   IF(NUZON(IZ,ILAY).GT.0)I=I+1
  ENDDO
 ENDDO
 IF(OUTPUTTYPE.EQ.1)THEN
  WRITE(IU,'(A)') TRIM(VTOS(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'
 ELSE
  WRITE(IU,'(A)') TRIM(VTOS(I))
 ENDIF
 
 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
    J=ILS(ILAY)
    WRITE(IU,'(A)') ' 1,'//TRIM(TPARAMETER)//','//TRIM(VTOS(J))//','// &
                                                   TRIM(VTOS(IZ+(ILAY-1)*NULIST))//',1.0,1.1,0.1,10.0,10.0,'//TRIM(VTOS(IG))//',1,'//TRIM(CULIST(I))
    !## fill in idf with zone number
    DO IROW=1,LIDF(ILAY)%NROW; DO ICOL=1,LIDF(ILAY)%NCOL
     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-IZONEOFFSET)ZIDF(ILAY)%X(ICOL,IROW)=IZ+(ILAY-1)*NULIST
     ENDIF  
    ENDDO; ENDDO
   ENDIF
  ENDDO
 ENDDO

 IF(OUTPUTTYPE.EQ.0)THEN; CLOSE(IU); IU=JU; ENDIF
 WRITE(IU,'(A)') TRIM(VTOS(NLAY))
 
 DO ILAY=1,NLAY
  J=ILS(ILAY)
  IF(.NOT.IDFWRITE(ZIDF(ILAY),TRIM(OFOLDER)//'\ZONE_L'//TRIM(VTOS(J))//'.IDF',1))STOP
  WRITE(IU,'(A)') TRIM(OFOLDER)//'\ZONE_L'//TRIM(VTOS(J))//'.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