!! 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 .
!!
!! 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)
!###======================================================================
IMPLICIT NONE
CHARACTER(LEN=*),INTENT(IN) :: PFOLDER,OFOLDER,TPARAMETER
INTEGER,INTENT(IN) :: NLAY,IZONEOFFSET,IGROUPOFFSET
INTEGER,INTENT(IN),DIMENSION(NLAY) :: LAYER
REAL(KIND=DP_KIND),INTENT(IN) :: MINF
CHARACTER(LEN=*),INTENT(IN),DIMENSION(:) :: FIDF,GRPN
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
REAL(KIND=DP_KIND) :: F
INTEGER(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: DXINT
CALL UTL_CREATEDIR(OFOLDER)
IU=UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OFOLDER)//'\PARAM.TXT',STATUS='REPLACE',ACTION='WRITE')
JU=UTL_GETUNIT(); OPEN(JU,FILE=TRIM(OFOLDER)//'\ZONES.TXT',STATUS='REPLACE',ACTION='WRITE')
ALLOCATE(IDF(NLAY)); DO I=1,SIZE(IDF); CALL IDFNULLIFY(IDF(I)); ENDDO
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
GRPNAMEFULL=TRIM(GRPN(J))//'_S'//TRIM(VTOS(IS))
!## 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)
WRITE(JU,'(A)') TRIM(IDF(II)%FNAME)
ENDIF
ENDIF
ENDDO
IF(TZ.LE.0)IG=IG-1
ENDDO
ENDDO
DEALLOCATE(DXINT)
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(VTOS(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='REPLACE',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(VTOS(I))//',(PST) PARAMETER ESTIMATION [-],1'
WRITE(IU,'(A)') '[<0=interpolation or >0]'
WRITE(IU,'(A)') ',,,,,,'
WRITE(IU,'(A)') ',0.0,0.1,0,0,1.0,0.0,3,,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(VTOS(ILAY))//','// &
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(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(VTOS(NLAY))
DO ILAY=1,NLAY
IF(.NOT.IDFWRITE(ZIDF(ILAY),TRIM(OFOLDER)//'\ZONE_L'//TRIM(VTOS(ILAY))//'.IDF',1))STOP
WRITE(IU,'(A)') TRIM(OFOLDER)//'\ZONE_L'//TRIM(VTOS(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