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

USE WINTERACTER
USE IMODVAR, ONLY : PI
USE MOD_PMANAGER_PAR, ONLY : PBMAN,TOP,BOT,TOPICS,TKHV,PEST,AFWATIDF,THFB,IUEXAMINE
USE MOD_IDF, ONLY : IDFNULLIFY,IDFALLOCATEX,IDFWRITE,IDFGETLOC,IDFIROWICOL,IDFREAD,IDFCOPY
USE MOD_ASC2IDF_PAR, ONLY : ELLIPS_IDF
USE MOD_UTL, ONLY : UTL_READPOINTER,UTL_READINITFILE,VTOS,UTL_REALTOSTRING,UTL_GETGAMMA,UTL_STDEF, &
  UTL_CREATEDIR,UTL_MATMUL,UTL_DIST_3D,UTL_SETUP_ROTATIONMATRIX,YMDHMSTOITIME,UTL_IDATETOGDATE,&
  UTL_COMPLETEDATE
USE MOD_KRIGING, ONLY : KRIGING_DIST  
USE MOD_IDF_PAR
USE MOD_LUDCMP 
USE MOD_IPEST_IES_PAR
USE MOD_QKSORT, ONLY : QKSORT
USE MOD_OSD, ONLY : OSD_TIMER

CHARACTER(LEN=3*256),PRIVATE :: LINE

CONTAINS
 
 !###======================================================================
 SUBROUTINE IMODBATCH_RUNFILE_INITPBMAN()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I
 
 PBMAN%MODFLOW=''
 PBMAN%IMOD_WQ=''
 PBMAN%MODFLOW6=''
 PBMAN%COUPLER=''
 PBMAN%MODFLOW_ADJ=''
 PBMAN%TOPMODEL=0
 PBMAN%IEXPORTMF2005=1
 
 !## initialize all pbman-variables
 PBMAN%BNDFILE=''
 PBMAN%OUTPUT=''
 PBMAN%IPESTPOUTPUT=''
 PBMAN%SEPM_FOLDER=''
 PBMAN%TIMFNAME=''
 PBMAN%GENFNAME=''

 PBMAN%IFORMAT=0; 
 PBMAN%ISOLVE=0
! PBMAN%SSYSTEM=0
 PBMAN%NSTEP=1
 PBMAN%NMULT=1.0D0
 PBMAN%ISTEADY=0
 PBMAN%IWINDOW=1
 PBMAN%ISS=0
 PBMAN%IDOUBLE=0
 PBMAN%ICONSISTENCY=1
 PBMAN%MINTHICKNESS=0.1D0
 PBMAN%THICKSTRT=0.1D0
 PBMAN%DBOTDEPTH=''
 PBMAN%DEFUNCONF=''
 PBMAN%PCKACTIVE=0
 PBMAN%DEFLAYER='0'
 PBMAN%DISTRCOND='0'
 PBMAN%DMMFILE=0
 AFWATIDF%FNAME=''; CALL IDFNULLIFY(AFWATIDF)
 PBMAN%DATEFORMAT=1
 
 !## steady-hfb package
 PBMAN%HFBTRAN=0
 
 PBMAN%NLAMBDASEARCH=3
 ALLOCATE(PBMAN%LAMBDA_TEST(3)); PBMAN%LAMBDA_TEST=[0.1D0,1.0D0,10.0D0]
 PBMAN%NLINESEARCH=1
 ALLOCATE(PBMAN%LINE_SEARCH(1)); PBMAN%LINE_SEARCH=[1.0D0]
 
 PBMAN%IPESTMETHOD=1
 PBMAN%NCPU=1
 PBMAN%NSWAIT=0
 PBMAN%CMDHIDE=1
 PBMAN%USESENS=0
 PBMAN%MU_INI=0.0D0

 !## set all to interpolate
 DO I=1,SIZE(PBMAN%INT); PBMAN%INT(I)=1; ENDDO
 PBMAN%IGENMF6=0
 PBMAN%SAVEEXCH=1
 
 PBMAN%NLAY=0
 PBMAN%NLAYIBND=0
 
 PBMAN%IPKS=0
 PBMAN%PARTOPT=0
 PBMAN%NRPROC=0
 PBMAN%LOADFILE=''
 PBMAN%PKSMERGE=1
 
 PBMAN%NSUBMODEL=1
 PBMAN%ISUBMODEL=1
 PBMAN%IPEST=0
 PBMAN%IPESTP=0
 PBMAN%IADJOINT=0
 PBMAN%IRANDOM=0
 PBMAN%RESTART=0
 PBMAN%SKIPMODEL=0
 PBMAN%IIES=0
 PBMAN%EIGV=99.9D0 
 PBMAN%ICLEAN=0
 PBMAN%AVGOBS=0
 
 PBMAN%IFILLHEAD=1

 PBMAN%IFVDL=0
 PBMAN%INT=1
 PBMAN%SSYSTEM=0
 PBMAN%MINKD=0.0D0
 PBMAN%MINC =0.0D0
 PBMAN%ICHKCHD=0
 PBMAN%ICONCHK=0
 PBMAN%IPUZZLE=0
 PBMAN%DWEL=1
 PBMAN%DISG=1
 PBMAN%DSFR=0
 PBMAN%DDRN=0.0D0
 PBMAN%QWEL=0.0D0
 PBMAN%ISAVEENDDATE=0
 PBMAN%INFFCT=0
 PBMAN%FORCECRD=0
 PBMAN%MERGELAYERS=0
 PBMAN%FLEXD=0
 PBMAN%EXAMINE=0
 PBMAN%SEPMODELS=''
 PBMAN%SYNCPARAM=1
 PBMAN%LAMBDAMODELS=''
 PBMAN%BUFMODELS=0.0D0
 PBMAN%NMODELS=0
! PBMAN%INTPARAM=0
 PBMAN%EXCLUDE=1
 PBMAN%BOSS=1
 PBMAN%NWORKERS=0
 PBMAN%RESETWORKERS=0
 PBMAN%BNDSUBMODEL=3
 PBMAN%APPLYCHD=1
 PBMAN%DEPTHUZF=''
 PBMAN%OUTERUPDATE=0
 PBMAN%NCYCLE=0
 PBMAN%KVAISKVV=0
 PBMAN%NWAVESUZF=40
 
 !## default packages
 PBMAN%ITT=9
 !## default number of selected interval
 PBMAN%IDT=1
 
 PBMAN%SDATE=INT(0,8)
 PBMAN%EDATE=INT(0,8)
 
 PBMAN%SMTYPE=0
 
 !## iMOD-WQ    
 PBMAN%BTN%NPROBS=1     !## btn
 PBMAN%BTN%NPRMAS=1     !## btn 
 PBMAN%ADV%MIXELM=-1    !## adv
 PBMAN%ADV%NADVFD=0     !## adv
 PBMAN%ADV%PERCEL=1     !## adv
 PBMAN%SSM%MXSS= 100000 !## ssm

 !## scr method (1)=iso 2=bjerrum
 PBMAN%SCR_IMETHOD=1
 !## pre-consolidation (0,1,2)
 PBMAN%SCR_ISTPCS=0
 !## output control (0,1,2)
 PBMAN%SCR_ISCROC=0
 !## est checking
 PBMAN%SCR_ESTCHK=1
 !## ontroduce feedback creep to hydraulic head
 PBMAN%SCR_FBFLAG=1
 
 !## use specific storage
 PBMAN%SPECIFICSTORAGE=0

 !## use default modflow6 solver congifuration
 PBMAN%NEWTON=0
 PBMAN%COMPLEXITY='SIMPLE'
 
 PBMAN%APPLYTC=0
 
 END SUBROUTINE IMODBATCH_RUNFILE_INITPBMAN
 
 !###======================================================================
 LOGICAL FUNCTION IMODBATCH_RUNFILE_READ(IU,IBATCH)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IBATCH,IU
 INTEGER :: I,J,N,IYR,IMH,IDY,IHR,IMT,ISC,IOS
 REAL(KIND=DP_KIND) :: X
 INTEGER,POINTER,DIMENSION(:) :: FLXILAY
 
 IMODBATCH_RUNFILE_READ=.FALSE.

 !## initialisation variables
 CALL IMODBATCH_RUNFILE_INITPBMAN()

 !## Look for Keywords
 PBMAN%IGENMF6=0; PBMAN%IFORMAT=2
 IF(UTL_READINITFILE('SIM_TYPE',LINE,IU,1))READ(LINE,*) PBMAN%IFORMAT
 SELECT CASE (PBMAN%IFORMAT)
  !## runfile
  CASE (1)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to a iMOD RUNFILE'
   IF(.NOT.UTL_READINITFILE('RUNFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'RUNFILE_OUT='//TRIM(PBMAN%RUNFILE)
  !## mf2005
  CASE (2)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to MODFLOW2005 files' 
   IF(.NOT.UTL_READINITFILE('NAMFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'NAMFILE_OUT='//TRIM(PBMAN%RUNFILE)
  !## mf6
  CASE (3)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to MODFLOW6 files' 
   IF(.NOT.UTL_READINITFILE('NAMFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'NAMFILE_OUT='//TRIM(PBMAN%RUNFILE)
   !## usage of genfiles for submodelling icm with mf6
   IF(UTL_READINITFILE('IGENMF6',LINE,IU,1))READ(LINE,*) PBMAN%IGENMF6
   !## read model network
   IF(PBMAN%IGENMF6.EQ.1)THEN
    PBMAN%IWINDOW=4
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IGENMF6=',PBMAN%IGENMF6
    IF(.NOT.UTL_READINITFILE('GENFNAME',LINE,IU,0))RETURN
    READ(LINE,*) PBMAN%GENFNAME; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'GENFNAME='//TRIM(PBMAN%GENFNAME)
    PBMAN%NSUBMODEL=0; IF(UTL_READINITFILE('NSUBMODEL',LINE,IU,1))READ(LINE,*) PBMAN%NSUBMODEL
    !## entre number of layers per submodel
    IF(PBMAN%NSUBMODEL.GT.0)THEN
     PBMAN%SMTYPE=1
     ALLOCATE(PBMAN%SM(PBMAN%NSUBMODEL))
     DO I=1,PBMAN%NSUBMODEL
      IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%SM(I)%ILAY,'ILAY_SM'//TRIM(VTOS(I)),0,EXCLVALUE=0))RETURN
      !## allocate memory for storage of submodel-idf files
      ALLOCATE(PBMAN%SM(I)%IDF(N)); DO J=1,N; CALL IDFNULLIFY(PBMAN%SM(I)%IDF(J)); ENDDO
     ENDDO
    ENDIF
    IF(UTL_READINITFILE('FORCECRD',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%FORCECRD
     IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'FORCECRD=',PBMAN%FORCECRD
    ENDIF    
    IF(UTL_READINITFILE('SAVEEXCH',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%SAVEEXCH
     IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SAVEEXCH=',PBMAN%SAVEEXCH
    ENDIF    
   ENDIF
  !## seawat - runfile
  CASE (4)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to an iMOD SEAWAT RUNFILE'
   IF(.NOT.UTL_READINITFILE('RUNFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'RUNFILE_OUT='//TRIM(PBMAN%RUNFILE)
  !## mt3d
  CASE (5)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to an iMOD MT3D RUNFILE'
   IF(.NOT.UTL_READINITFILE('RUNFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'RUNFILE_OUT='//TRIM(PBMAN%RUNFILE)
   IF(.NOT.UTL_READINITFILE('FLOW_RESULT_DIR',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%FLOW_RESULT_DIR; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'FLOW_RESULT_DIR='//TRIM(PBMAN%FLOW_RESULT_DIR)

   IF(UTL_READINITFILE('NPROBS',LINE,IU,1))READ(LINE,*) PBMAN%BTN%NPROBS
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'NPROBS='//TRIM(VTOS(PBMAN%BTN%NPROBS))
   IF(UTL_READINITFILE('NPRMAS',LINE,IU,1))READ(LINE,*) PBMAN%BTN%NPRMAS
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'NPRMAS='//TRIM(VTOS(PBMAN%BTN%NPRMAS))
   IF(UTL_READINITFILE('MIXELM',LINE,IU,1))READ(LINE,*) PBMAN%ADV%MIXELM
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'MIXELM='//TRIM(VTOS(PBMAN%ADV%MIXELM))
   IF(UTL_READINITFILE('NADVFD',LINE,IU,1))READ(LINE,*) PBMAN%ADV%NADVFD
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'NADVFD='//TRIM(VTOS(PBMAN%ADV%NADVFD))
   IF(UTL_READINITFILE('PERCEL',LINE,IU,1))READ(LINE,*) PBMAN%ADV%PERCEL
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'PERCEL='//TRIM(UTL_REALTOSTRING(PBMAN%ADV%PERCEL))
   IF(UTL_READINITFILE('MXSS',LINE,IU,1))READ(LINE,*) PBMAN%SSM%MXSS
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'MXSS='//TRIM(VTOS(PBMAN%SSM%MXSS))
  !## seawat - namfile
  CASE (6)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to an iMOD SEAWAT NAMFILE'
   IF(.NOT.UTL_READINITFILE('NAMFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'NAMFILE_OUT='//TRIM(PBMAN%RUNFILE)
  CASE (7)
   IF(IBATCH.EQ.1)WRITE(*,'(/A/)') 'Export to MODFLOW1988 files'
   IF(.NOT.UTL_READINITFILE('BASFILE_OUT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%RUNFILE; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'BASFILE_OUT='//TRIM(PBMAN%RUNFILE)
  CASE DEFAULT
   IF(IBATCH.EQ.1)WRITE(*,'(/A,I5/)') 'Incorrect Output option, iformat=', PBMAN%IFORMAT   
   IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Incorrect Output option, iformat='//TRIM(VTOS(PBMAN%IFORMAT)),'Error'); RETURN
 END SELECT

 IF(UTL_READINITFILE('OUTPUT_FOLDER',LINE,IU,1))THEN  
  READ(LINE,*) PBMAN%OUTPUT; IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'OUTPUT_FOLDER='//TRIM(PBMAN%OUTPUT)
 ELSE !## create foldername
  PBMAN%OUTPUT=TRIM(PBMAN%RUNFILE(:INDEX(PBMAN%RUNFILE,'\',.TRUE.)-1))
 ENDIF
 IF(UTL_READINITFILE('IPESTP_FOLDER',LINE,IU,1))THEN  
  READ(LINE,*) PBMAN%IPESTPOUTPUT; IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'IPESTP_FOLDER='//TRIM(PBMAN%IPESTPOUTPUT)
 ELSE !## create foldername
  PBMAN%IPESTPOUTPUT=TRIM(PBMAN%RUNFILE(:INDEX(PBMAN%RUNFILE,'\',.TRUE.)-1))//'\IPEST'
 ENDIF
 
 IF(PBMAN%IGENMF6.EQ.0)THEN
  IF(UTL_READINITFILE('NETWORKIDF',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%BNDFILE
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'NETWORKIDF='//TRIM(PBMAN%BNDFILE)
   PBMAN%IWINDOW=3
  ELSE
   IF(UTL_READINITFILE('WINDOW',LINE,IU,1))THEN
    PBMAN%IWINDOW=2
    READ(LINE,*) PBMAN%XMIN,PBMAN%YMIN,PBMAN%XMAX,PBMAN%YMAX
    IF(IBATCH.EQ.1)WRITE(*,'(A,4F15.3)') 'WINDOW=',PBMAN%XMIN,PBMAN%YMIN,PBMAN%XMAX,PBMAN%YMAX
    IF(.NOT.UTL_READINITFILE('CELLSIZE',LINE,IU,0))RETURN
    READ(LINE,*) PBMAN%CELLSIZE; IF(IBATCH.EQ.1)WRITE(*,'(A,F15.3)') 'CELLSIZE=',PBMAN%CELLSIZE
   ELSE
    !## cellsize optional
    IF(UTL_READINITFILE('CELLSIZE',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%CELLSIZE
     IF(IBATCH.EQ.1)THEN
      WRITE(*,'(A,F15.3)') 'CELLSIZE=',PBMAN%CELLSIZE
      IF(PBMAN%CELLSIZE.LE.0.0)STOP 'CELLSIZE LE 0.0'
     ENDIF
    ENDIF
   ENDIF
   PBMAN%BUFFER=0.0D0
   IF(UTL_READINITFILE('BUFFER',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%BUFFER
    IF(IBATCH.EQ.1)WRITE(*,'(A,F15.3)') 'BUFFER=',PBMAN%BUFFER
   ENDIF
   IF(PBMAN%BUFFER.GT.0.0D0)THEN
    PBMAN%BUFFERCS=0.0D0; IF(UTL_READINITFILE('BUFFERCS',LINE,IU,1))READ(LINE,*) PBMAN%BUFFERCS
    IF(IBATCH.EQ.1)THEN
     WRITE(*,'(A,F15.3)') 'BUFFERCS=',PBMAN%BUFFERCS
     IF(PBMAN%BUFFERCS.LT.PBMAN%CELLSIZE)THEN
      WRITE(*,'(/A,2(F10.3,A)/)') 'BUFFERCS (',PBMAN%BUFFERCS,') NEED TO BE AT LEAST EQUAL OR LARGER THAN CELLSIZE (',PBMAN%CELLSIZE,')'; STOP
     ENDIF
    ENDIF
   ENDIF
  ENDIF
  IF(UTL_READINITFILE('APPLYCHD',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%APPLYCHD
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'APPLYCHD='//TRIM(VTOS(PBMAN%APPLYCHD))
  ENDIF
 ENDIF
    
 SELECT CASE (PBMAN%IFORMAT)
  CASE (1,2,3,4,6,7)
   !## more specific options to create a runfile
   IF(UTL_READINITFILE('ISS',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%ISS
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ISS=',PBMAN%ISS
   ENDIF
  !## always transient for mt3d
  CASE (5)
   PBMAN%ISS=1
 END SELECT
 
 !## more specific options to create a runfile
 IF(UTL_READINITFILE('NLAY',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%NLAY
  IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NLAY=',PBMAN%NLAY
  IF(UTL_READINITFILE('NLAYIBND',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%NLAYIBND
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NLAYIBND=',PBMAN%NLAYIBND
  ENDIF
 ENDIF

 !## transient hfb in prj-file
 IF(UTL_READINITFILE('HFBTRAN',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%HFBTRAN; IF(IBATCH.EQ.1)WRITE(*,'(A,A)') 'HFBTRAN=',TRIM(VTOS(PBMAN%HFBTRAN))
  IF(PBMAN%HFBTRAN.EQ.1)TOPICS(THFB)%TIMDEP=.TRUE.
 ENDIF

 !## transient                              
 IF(PBMAN%ISS.EQ.1)THEN
  !## tim-file?
  IF(UTL_READINITFILE('TIMFNAME',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%TIMFNAME; IF(IBATCH.EQ.1)WRITE(*,'(A,A)') 'TIMFNAME=',TRIM(PBMAN%TIMFNAME)
  ELSE
   !## sdate
   IF(.NOT.UTL_READINITFILE('SDATE',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%SDATE; PBMAN%SDATE=UTL_COMPLETEDATE(PBMAN%SDATE); IF(IBATCH.EQ.1)WRITE(*,'(A,I14)') 'SDATE=',PBMAN%SDATE

   !## change 31st of December at 0hr,0mt and 0s into 1st of January next year
   CALL UTL_IDATETOGDATE(PBMAN%SDATE,IYR,IMH,IDY,IHR,IMT,ISC)
   IF(IMH.EQ.12.AND.IDY.EQ.31.AND.IHR.EQ.0.AND.IMT.EQ.0.AND.ISC.EQ.0)THEN
    IYR=IYR+1
    IDY=1
    IMH=1
    PBMAN%SDATE=YMDHMSTOITIME(IYR,IMH,IDY,IHR,IMT,ISC)
    IF(IBATCH.EQ.1)WRITE(*,'(A,I14)') 'SDATE (corrected)=',PBMAN%SDATE
   ENDIF

   !## edate
   IF(.NOT.UTL_READINITFILE('EDATE',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%EDATE; PBMAN%EDATE=UTL_COMPLETEDATE(PBMAN%EDATE); IF(IBATCH.EQ.1)WRITE(*,'(A,I14)') 'EDATE=',PBMAN%EDATE

   !## change 31st of December at 0hr,0mt and 0s into 1st of January next year
   CALL UTL_IDATETOGDATE(PBMAN%EDATE,IYR,IMH,IDY,IHR,IMT,ISC)
   IF(IMH.EQ.12.AND.IDY.EQ.31.AND.IHR.EQ.0.AND.IMT.EQ.0.AND.ISC.EQ.0)THEN
    IYR=IYR+1
    IDY=1
    IMH=1
    PBMAN%EDATE=YMDHMSTOITIME(IYR,IMH,IDY,IHR,IMT,ISC)
    IF(IBATCH.EQ.1)WRITE(*,'(A,I14)') 'EDATE (corrected)=',PBMAN%EDATE
   ENDIF
   
   !## type of interval
   IF(.NOT.UTL_READINITFILE('ITT',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%ITT; IF(IBATCH.EQ.1)WRITE(*,'(A,I14)') 'ITT=',PBMAN%ITT
   SELECT CASE (PBMAN%ITT)
    CASE (1:8)
     !## timesteps
     IF(.NOT.UTL_READINITFILE('IDT',LINE,IU,0))RETURN
     READ(LINE,*) PBMAN%IDT; IF(IBATCH.EQ.1)WRITE(*,'(A,I14)') 'IDT=',PBMAN%IDT
   END SELECT
   SELECT CASE (PBMAN%ITT)
    CASE (1)
     IF(IBATCH.EQ.1)WRITE(*,'(A)')      'Time step: ',PBMAN%IDT,' minutes'
    CASE (2)
     IF(IBATCH.EQ.1)WRITE(*,'(A,I5,A)') 'Time step: ',PBMAN%IDT,' hours'
    CASE (3)
     IF(IBATCH.EQ.1)WRITE(*,'(A,I5,A)') 'Time step: ',PBMAN%IDT,' days'
    CASE (4)
     IF(IBATCH.EQ.1)WRITE(*,'(A,I5,A)') 'Time step: ',PBMAN%IDT,' weeks'
    CASE (5)
     IF(IBATCH.EQ.1)WRITE(*,'(A)')      'Time step: decades (10 days)' 
    CASE (6)
     IF(IBATCH.EQ.1)WRITE(*,'(A)')      'Time step: every 14/28'
    CASE (7)
     IF(IBATCH.EQ.1)WRITE(*,'(A,I5,A)') 'Time step: ',PBMAN%IDT,' months'
    CASE (8)
     IF(IBATCH.EQ.1)WRITE(*,'(A,I5,A)') 'Time step: ',PBMAN%IDT,' years'
    CASE (9)
     IF(IBATCH.EQ.1)WRITE(*,'(A)')      'Time step: depending on packages'; PBMAN%IDT=1
   END SELECT

   IF(UTL_READINITFILE('NSTEP',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%NSTEP
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NSTEP=',PBMAN%NSTEP
   ENDIF
   IF(UTL_READINITFILE('NMULT',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%NMULT
    IF(IBATCH.EQ.1)WRITE(*,'(A,F10.2)') 'NMULT=',PBMAN%NMULT
   ENDIF
   SELECT CASE (PBMAN%IFORMAT)
    CASE (1,2,3,6)
    IF(UTL_READINITFILE('ISTEADY',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%ISTEADY
     IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ISTEADY=',PBMAN%ISTEADY
    ENDIF
   END SELECT 
   
  ENDIF    
 
 ENDIF
  
 IF(UTL_READINITFILE('SSYSTEM',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%SSYSTEM
  IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SSYSTEM=',PBMAN%SSYSTEM
 ENDIF
 !## get isave options
 DO I=1,SIZE(TOPICS)
  IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%ISAVE(I)%ILAY,'SAVE'//TRIM(TOPICS(I)%CMOD),1,EXCLVALUE=0))RETURN
 ENDDO
 !## try to read old fashioned save settings
 IF(.NOT.UTL_READPOINTER(IU,N,FLXILAY,'SAVEFLX',1,EXCLVALUE=0))RETURN
 IF(ASSOCIATED(FLXILAY))THEN
  IF(ASSOCIATED(PBMAN%ISAVE(TKHV)%ILAY))DEALLOCATE(PBMAN%ISAVE(TKHV)%ILAY)
  ALLOCATE(PBMAN%ISAVE(TKHV)%ILAY(SIZE(FLXILAY)))
  PBMAN%ISAVE(TKHV)%ILAY=FLXILAY; DEALLOCATE(FLXILAY)
 ENDIF
 
 !## get interpolation options
 SELECT CASE (PBMAN%IFORMAT)
  CASE (2,3,6,7)
   IF(UTL_READINITFILE('PCKACTIVE',LINE,IU,1))READ(LINE,*) PBMAN%PCKACTIVE
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'PCKACTIVE='//TRIM(VTOS(PBMAN%PCKACTIVE))
   IF(UTL_READINITFILE('DISTRCOND',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%DISTRCOND
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'DISTRCOND='//TRIM(PBMAN%DISTRCOND)
   ENDIF
   IF(UTL_READINITFILE('IDEFLAYER',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%DEFLAYER
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'IDEFLAYER='//TRIM(PBMAN%DEFLAYER)
   ENDIF
   DO I=1,SIZE(TOPICS)
    IF(UTL_READINITFILE('INT'//TRIM(TOPICS(I)%CMOD),LINE,IU,1))READ(LINE,*) PBMAN%INT(I)
    IF(PBMAN%INT(I).EQ.0.AND.IBATCH.EQ.1)WRITE(*,'(A,I10)') 'INT'//TRIM(TOPICS(I)%CMOD)//'=',PBMAN%INT(I)
   ENDDO
 END SELECT
  
 IF(PBMAN%IFORMAT.EQ.1.OR.PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.6)THEN
  IF(UTL_READINITFILE('ICONCHK',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%ICONCHK
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ICONCHK=',PBMAN%ICONCHK
  ENDIF
  IF(UTL_READINITFILE('IDOUBLE',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%IDOUBLE
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IDOUBLE=',PBMAN%IDOUBLE
  ENDIF
  IF(UTL_READINITFILE('IPKS',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%IPKS
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IPKS=',PBMAN%IPKS
  ENDIF
  IF(PBMAN%IPKS.EQ.1)THEN
   IF(UTL_READINITFILE('NRPROC',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%NRPROC
    PBMAN%NRPROC=MAX(1,PBMAN%NRPROC)
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NRPROC=',PBMAN%NRPROC
   ENDIF
   IF(UTL_READINITFILE('PKSMERGE',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%PKSMERGE
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'PKSMERGE='//TRIM(VTOS(PBMAN%PKSMERGE))
   ENDIF
  ENDIF
 ENDIF
 
 SELECT CASE (PBMAN%IFORMAT)
  CASE (1,2,3)
   IF(UTL_READINITFILE('ISAVEENDDATE',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%ISAVEENDDATE
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ISAVEENDDATE=',PBMAN%ISAVEENDDATE
   ENDIF
  !## for iMOD-WQ allways apply timestamp at the end of the simulation
  CASE (4,5,6,7)
   PBMAN%ISAVEENDDATE=1
 END SELECT

 !## mf2005/modflow6/seawat
 SELECT CASE(PBMAN%IFORMAT)
  CASE(2,3,6)
   IF(UTL_READINITFILE('KVAISKVV',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%KVAISKVV; WRITE(*,'(A)') 'KVAISKVV='//TRIM(VTOS(PBMAN%KVAISKVV))
   ENDIF
 END SELECT
 
 !## mf2005/modflow6 
 PBMAN%EXAMINE=0.0D0; IUEXAMINE=0
 IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.3)THEN
  IF(UTL_READINITFILE('EXAMINE',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%EXAMINE(1),PBMAN%EXAMINE(2)
   IF(PBMAN%EXAMINE(1).LT.0.AND.PBMAN%EXAMINE(2).LT.0)THEN
    IF(IBATCH.EQ.1)WRITE(*,'(A,2I10)') 'EXAMINE=',ABS(INT(PBMAN%EXAMINE(1))),ABS(INT(PBMAN%EXAMINE(2)))
   ELSEIF(PBMAN%EXAMINE(1).GT.0.0D0.AND.PBMAN%EXAMINE(2).GT.0.0D0)THEN
    IF(IBATCH.EQ.1)WRITE(*,'(A,2F15.2)') 'EXAMINE=',PBMAN%EXAMINE(1),PBMAN%EXAMINE(2)
   ELSE
    WRITE(*,'(/A/)') '>>> ERROR READING EXAMINE, BOTH NEGATIVE OR BOTH POSITIVE VALUES <<<'; STOP
   ENDIF
  ENDIF
 ENDIF
 SELECT CASE (PBMAN%IFORMAT)
  CASE (2,3,6)
   IF(UTL_READINITFILE('NCYCLE',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%NCYCLE; WRITE(*,'(A)') 'NCYCLE='//TRIM(VTOS(PBMAN%NCYCLE))
    
    IF(UTL_READINITFILE('SEPM_FOLDER',LINE,IU,1))THEN  
     READ(LINE,*) PBMAN%SEPM_FOLDER; IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'SEPM_FOLDER='//TRIM(PBMAN%SEPM_FOLDER)
    ELSE !## create foldername
     PBMAN%SEPM_FOLDER= PBMAN%OUTPUT
    ENDIF
    
    IF(UTL_READINITFILE('OUTERUPDATE',LINE,IU,1))READ(LINE,*) PBMAN%OUTERUPDATE
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'OUTERUPDATE='//TRIM(VTOS(PBMAN%OUTERUPDATE))
    IF(PBMAN%OUTERUPDATE.EQ.1)THEN
     IF(UTL_READINITFILE('SYNCPARAM',LINE,IU,1))READ(LINE,*) PBMAN%SYNCPARAM
     IF(PBMAN%SYNCPARAM.LT.1.OR.PBMAN%SYNCPARAM.GT.3)THEN
      WRITE(*,'(1X,A/)') '>>> SYNPARAM need to be 1,2 or 3 <<<'; STOP
     ENDIF
     IF(IBATCH.EQ.1)WRITE(*,'(A)') 'SYNCPARAM='//TRIM(VTOS(PBMAN%SYNCPARAM))
    ENDIF
    IF(UTL_READINITFILE('NMODELS',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%NMODELS; WRITE(*,'(A)') 'NMODELS='//TRIM(VTOS(PBMAN%NMODELS))
    ELSEIF(UTL_READINITFILE('SEPMODELS',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%SEPMODELS
     IF(IBATCH.EQ.1)WRITE(*,'(A)') 'SEPMODELS='//TRIM(PBMAN%SEPMODELS)    
    ELSEIF(UTL_READINITFILE('LAMBDAMODELS',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%LAMBDAMODELS
     IF(IBATCH.EQ.1)WRITE(*,'(A)') 'LAMBDAMODELS='//TRIM(PBMAN%LAMBDAMODELS)    
    ELSE
     WRITE(*,'(/A/)') '>>> ERROR, NEED TO SPECIFY NMODELS OR SEPMODELS OR LAMBDAMODELS <<<'; STOP   
    ENDIF
    IF(UTL_READINITFILE('BUFMODELS',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%BUFMODELS
     IF(IBATCH.EQ.1)WRITE(*,'(A,F15.2)') 'BUFMODELS=',PBMAN%BUFMODELS
    ENDIF
!    IF(UTL_READINITFILE('INTPARAM',LINE,IU,1))READ(LINE,*) PBMAN%INTPARAM
!    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'INTPARAM='//TRIM(VTOS(PBMAN%INTPARAM))
    IF(UTL_READINITFILE('SKIPMODEL',LINE,IU,1))READ(LINE,*) PBMAN%SKIPMODEL
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SKIPMODEL=',PBMAN%SKIPMODEL
    IF(UTL_READINITFILE('EXCLUDE',LINE,IU,1))READ(LINE,*) PBMAN%EXCLUDE
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'EXCLUDE='//TRIM(VTOS(PBMAN%EXCLUDE))
    IF(UTL_READINITFILE('NWORKERS',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%NWORKERS; PBMAN%NWORKERS=MAX(0,PBMAN%NWORKERS)
     IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NWORKERS=',PBMAN%NWORKERS
    ENDIF
    IF(PBMAN%NWORKERS.LE.0)THEN
     IF(UTL_READINITFILE('WORKER',LINE,IU,1))THEN
      READ(LINE,*) PBMAN%WORKER; IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'WORKER=',PBMAN%WORKER
      !## it is a worker
      PBMAN%BOSS=0
      IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%IMODEL,'IMODEL',0))RETURN
     ELSE
      IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%IMODEL,'IMODEL',1))RETURN
     ENDIF
    ENDIF
    IF(UTL_READINITFILE('RESETWORKERS',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%RESETWORKERS; IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'RESETWORKERS=',PBMAN%RESETWORKERS
    ENDIF
    IF(UTL_READINITFILE('BNDSUBMODEL',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%BNDSUBMODEL; IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'BNDSUBMODEL=',PBMAN%BNDSUBMODEL
    ENDIF
   ENDIF
 END SELECT

 !## usage of PKS
 IF(PBMAN%IPKS.EQ.1)THEN 
  IF(UTL_READINITFILE('PARTOPT',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%PARTOPT
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'PARTOPT=',PBMAN%PARTOPT
   PBMAN%PARTOPT=PBMAN%PARTOPT-1
  ENDIF
  SELECT CASE (PBMAN%PARTOPT)
   CASE (0); WRITE(*,'(/A)') 'UNIFORM SUBMODELLING SELECTED'
   CASE (1); WRITE(*,'(/A)') 'RCB SUBMODELLING SELECTED'
   CASE (2); WRITE(*,'(/A)') 'RCB WEIGHTED SUBMODELLING SELECTED'
   CASE DEFAULT; STOP 'WRONG PARTOPT SELECTED, PICK 1,2 OR 3'
  END SELECT
   
  IF(PBMAN%PARTOPT.NE.0)THEN
   IF(UTL_READINITFILE('LOADFILE',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%LOADFILE
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'LOADFILE='//TRIM(PBMAN%LOADFILE)
   ENDIF
  ENDIF
 
 ENDIF
 
 !## specify unconfined/parameter estimation
 IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.3.OR.PBMAN%IFORMAT.EQ.6)THEN
  IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%UNCONFINED,'UNCONFINED',1))RETURN
  !## if any of the entries of unconfined shows a 2 (unconfined simulation), enforce ichkchd=1
  PBMAN%ICHKCHD=-1
  IF(ASSOCIATED(PBMAN%UNCONFINED))THEN
   DO I=1,SIZE(PBMAN%UNCONFINED)
    IF(PBMAN%UNCONFINED(I).EQ.1)THEN; PBMAN%ICHKCHD=1; EXIT; ENDIF
   ENDDO
   IF(IBATCH.EQ.1.AND.PBMAN%ICHKCHD.EQ.1)WRITE(*,'(A,I10)') '(enforced for unconfined conditions) ICHKCHD=',PBMAN%ICHKCHD
  ENDIF
  IF(PBMAN%ICHKCHD.EQ.-1)THEN
   IF(UTL_READINITFILE('ICHKCHD',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%ICHKCHD
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ICHKCHD=',PBMAN%ICHKCHD
   ENDIF
  ENDIF
  IF(PBMAN%ICHKCHD.EQ.-1)PBMAN%ICHKCHD=0
  IF(UTL_READINITFILE('ICONSISTENCY',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%ICONSISTENCY
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ICONSISTENCY=',PBMAN%ICONSISTENCY
  ENDIF
  IF(UTL_READINITFILE('THICKSTRT',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%THICKSTRT
   IF(IBATCH.EQ.1)WRITE(*,'(A,F10.3)') 'THICKSTRT=',PBMAN%THICKSTRT
   IF(PBMAN%THICKSTRT.LE.0.0D0)THEN
    WRITE(*,'(/A/)') 'THICKSTRT needs to be > 0.0'; STOP
   ENDIF
  ENDIF
  IF(UTL_READINITFILE('DBOTDEPTH',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%DBOTDEPTH
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'DBOTDEPTH='//TRIM(PBMAN%DBOTDEPTH)
   READ(PBMAN%DBOTDEPTH,*,IOSTAT=IOS) X
   IF(IOS.EQ.0.AND.X.LE.0.0D0)THEN
    WRITE(*,'(/A/)') 'DBOTDEPTH needs to be > 0.0'; STOP
   ENDIF
  ENDIF
  IF(PBMAN%ICONSISTENCY.EQ.2)THEN
   IF(UTL_READINITFILE('MINTHICKNESS',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%MINTHICKNESS
    IF(IBATCH.EQ.1)WRITE(*,'(A,G15.7)') 'MINTHICKNESS=',PBMAN%MINTHICKNESS
   ENDIF
  ENDIF
  IF(UTL_READINITFILE('SPECIFICSTORAGE',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%SPECIFICSTORAGE
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SPECIFICSTORAGE=',PBMAN%SPECIFICSTORAGE
  ENDIF
  IF(UTL_READINITFILE('APPLYTC',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%APPLYTC
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'APPLYTC=',PBMAN%APPLYTC
  ENDIF
 ENDIF
  
 !## if not defined, than all confined
 IF(.NOT.ASSOCIATED(PBMAN%UNCONFINED))THEN; ALLOCATE(PBMAN%UNCONFINED(1)); PBMAN%UNCONFINED(1)=0; ENDIF
 
 IF(PBMAN%IFORMAT.EQ.3)THEN
  IF(UTL_READINITFILE('DEFUNCONF',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%DEFUNCONF
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'DEFUNCONF='//TRIM(PBMAN%DEFUNCONF)
   READ(PBMAN%DEFUNCONF,*,IOSTAT=IOS) X
   IF(IOS.EQ.0.AND.X.LT.0.0D0)THEN
    WRITE(*,'(/A/)') 'DEFUNCONF needs to be >= 0'; STOP
   ENDIF
  ENDIF
 ENDIF
 
 !## usage of dflow-fm (modflow2005/modflow6)
 IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.3)THEN
  IF(UTL_READINITFILE('DMMFILE',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%DMMFILE
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'DMMFILE=',PBMAN%DMMFILE
  ENDIF
  IF(PBMAN%DMMFILE.EQ.1)THEN
   IF(UTL_READINITFILE('AFWATIDF',LINE,IU,1))THEN
    READ(LINE,*) AFWATIDF%FNAME
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'AFWATIDF='//TRIM(AFWATIDF%FNAME)
   ENDIF
  ENDIF
 ENDIF
 !## usage of subscreep (imod-wq)
 IF(PBMAN%IFORMAT.EQ.6)THEN
  IF(UTL_READINITFILE('SCR_IMETHOD',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%SCR_IMETHOD
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SCR_IMETHOD=',PBMAN%SCR_IMETHOD
   IF(PBMAN%SCR_IMETHOD.LT.1.OR.PBMAN%SCR_IMETHOD.GT.2)THEN
    STOP 'SCR_METHOD needs to be 1 or 2'
   ENDIF
  ENDIF
  IF(UTL_READINITFILE('SCR_ISTPCS',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%SCR_ISTPCS
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SCR_ICTPCS=',PBMAN%SCR_ISTPCS
   IF(PBMAN%SCR_ISTPCS.LT.0.OR.PBMAN%SCR_ISTPCS.GT.3)THEN
    STOP 'SCR_ISTPCS needs to be 0,1,2 or 3'
   ENDIF
  ENDIF
  IF(UTL_READINITFILE('SCR_ISCROC',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%SCR_ISCROC
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SCR_ISCROC=',PBMAN%SCR_ISCROC
   IF(PBMAN%SCR_ISCROC.LT.0.OR.PBMAN%SCR_ISCROC.GT.2)THEN
    STOP 'SCR_ISCROC NEEDS TO BE 0,1 OR 2'
   ENDIF
  ENDIF 
  IF(UTL_READINITFILE('SCR_ESTCHK',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%SCR_ESTCHK
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SCR_ESTCHK=',PBMAN%SCR_ESTCHK
   IF(PBMAN%SCR_ESTCHK.LT.0.OR.PBMAN%SCR_ESTCHK.GT.1)THEN
    WRITE(*,'(/A/)') '>>> SCR_ESTCHK NEED TO BE 0 OR 1 <<<'; PAUSE; STOP
   ENDIF
  ENDIF
  IF(UTL_READINITFILE('SCR_FBFLAG',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%SCR_FBFLAG; 
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'SCR_FBFLAG=',PBMAN%SCR_FBFLAG
   IF(PBMAN%SCR_FBFLAG.LT.0.OR.PBMAN%SCR_FBFLAG.GT.1)THEN
    WRITE(*,'(/A/)') '>>> SCR_FBFLAG NEED TO BE 0 OR 1 <<<'; PAUSE; STOP
   ENDIF
  ENDIF
 ENDIF

 !## usage of ipest
 SELECT CASE (PBMAN%IFORMAT)
  CASE (1,2)
   IF(UTL_READINITFILE('IPEST',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%IPEST
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IPEST=',PBMAN%IPEST
   ENDIF
  CASE (7)
   IF(UTL_READINITFILE('IADJOINT',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%IADJOINT
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IADJOINT=',PBMAN%IADJOINT
   ENDIF
 END SELECT

 IF(UTL_READINITFILE('IEXPORT',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%IEXPORTMF2005
  IF(IBATCH.EQ.1)WRITE(*,'(2A)') 'IEXPORT='//TRIM(VTOS(PBMAN%IEXPORTMF2005))
 ENDIF
 
 !## usage of parrallel ipest or ies for modflow2005/modflow6
 IF(PBMAN%IPEST.EQ.0.AND.PBMAN%IADJOINT.EQ.0)THEN
  IF(UTL_READINITFILE('IIES',LINE,IU,1))READ(LINE,*) PBMAN%IIES
  IF(PBMAN%IIES.AND.IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IIES=',PBMAN%IIES
  IF(PBMAN%IIES.EQ.0)THEN
   IF(UTL_READINITFILE('IPESTP',LINE,IU,1))READ(LINE,*) PBMAN%IPESTP
   IF(PBMAN%IPESTP.EQ.1)THEN
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IPESTP=',PBMAN%IPESTP
   ENDIF
  ENDIF
  IF(PBMAN%IPESTP.EQ.1.OR.PBMAN%IIES.EQ.1)THEN

   IF(UTL_READINITFILE('IRANDOM',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%IRANDOM
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IRANDOM=',PBMAN%IRANDOM
   ENDIF

   IF(UTL_READINITFILE('MU_INI',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%MU_INI
    IF(IBATCH.EQ.1)WRITE(*,'(A,F10.3)') 'MU_INI=',PBMAN%MU_INI
    IF(PBMAN%MU_INI.LE.0.0D0)THEN
     WRITE(*,'(/A/)') 'MU_INI needs to be > 0.0'; STOP
    ENDIF
   ENDIF
 
   IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%LAMBDA_TEST,'LAMBDA_TEST',1))RETURN
   !## sort from small to big
   CALL QKSORT(SIZE(PBMAN%LAMBDA_TEST),PBMAN%LAMBDA_TEST)
   PBMAN%NLAMBDASEARCH=SIZE(PBMAN%LAMBDA_TEST)
   
   IF(.NOT.UTL_READPOINTER(IU,N,PBMAN%LINE_SEARCH,'LINESEARCH',1))RETURN
   !## sort from small to big
   CALL QKSORT(SIZE(PBMAN%LINE_SEARCH),PBMAN%LINE_SEARCH)
   PBMAN%NLINESEARCH=SIZE(PBMAN%LINE_SEARCH)

   IF(UTL_READINITFILE('CMDHIDE',LINE,IU,1))READ(LINE,*) PBMAN%CMDHIDE
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'CMDHIDE=',PBMAN%CMDHIDE
   IF(UTL_READINITFILE('NCPU',LINE,IU,1))READ(LINE,*) PBMAN%NCPU
   PBMAN%NCPU=MAX(1,PBMAN%NCPU)
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NCPU=',PBMAN%NCPU
   IF(UTL_READINITFILE('NSWAIT',LINE,IU,1))READ(LINE,*) PBMAN%NSWAIT
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NSWAIT=',PBMAN%NSWAIT
   IF(UTL_READINITFILE('SVD_EIGV',LINE,IU,1))READ(LINE,*) PBMAN%EIGV
   IF(IBATCH.EQ.1)WRITE(*,'(A,F10.4)') 'SVD_EIGV=',PBMAN%EIGV
   IF(UTL_READINITFILE('USESENS',LINE,IU,1))READ(LINE,*) PBMAN%USESENS
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'USESENS=',PBMAN%USESENS
   IF(UTL_READINITFILE('AVGOBS',LINE,IU,1))READ(LINE,*) PBMAN%AVGOBS
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'AVGOBS=',PBMAN%AVGOBS
   IF(UTL_READINITFILE('ICLEAN',LINE,IU,1))READ(LINE,*) PBMAN%ICLEAN
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ICLEAN=',PBMAN%ICLEAN
   IF(UTL_READINITFILE('RESTART',LINE,IU,1))READ(LINE,*) PBMAN%RESTART
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'RESTART=',PBMAN%RESTART
   IF(UTL_READINITFILE('IPESTMETHOD',LINE,IU,1))READ(LINE,*) PBMAN%IPESTMETHOD
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IPESTMETHOD=',PBMAN%IPESTMETHOD
   SELECT CASE (PBMAN%IPESTMETHOD)
    CASE(1); WRITE(*,'(A/)') '>>> Standard GLM selected <<<'
    CASE(2); WRITE(*,'(A/)') '>>> Pattern search selected <<<'
    CASE(3); WRITE(*,'(A/)') '>>> Neural-Network GLM-selected <<<'
    CASE DEFAULT; WRITE(*,'(A/)') '>>> Pick 1,2,3 as option <<<'; PAUSE; STOP
   END SELECT
  ENDIF
 ENDIF

 !## write warning for system in case of ipest - ssystem=sumsystem
 IF(PBMAN%IPEST+PBMAN%IPESTP+PBMAN%IIES.GT.0.AND.PBMAN%SSYSTEM.EQ.1)THEN
  IF(IBATCH.EQ.1)WRITE(*,'(/1X,A/)') '>>> WARNING: summing systems might influence iPEST(P)/IIES if you optimize packages <<<'
 ENDIF
 IF(PBMAN%IPKS.EQ.1.AND.PBMAN%SSYSTEM.EQ.0)THEN
  IF(IBATCH.EQ.1)WRITE(*,'(/1X,A/)') '>>> WARNING: not summing systems causes PKS to merge the results improperly. It is advised to set ssystem = 1 <<<'
 ENDIF
  
 IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.6)THEN
  IF(UTL_READINITFILE('OVERLANDSTORAGE',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%OVERLANDSTORAGE
   IF(IBATCH.EQ.1)WRITE(*,'(A)') 'OVERLANDSTORAGE='//TRIM(VTOS(PBMAN%OVERLANDSTORAGE))
  ENDIF
 ENDIF
 
 IF(PBMAN%IFORMAT.EQ.1.OR.PBMAN%IFORMAT.EQ.2)THEN
  IF(UTL_READINITFILE('MINKD',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%MINKD
   IF(IBATCH.EQ.1)WRITE(*,'(A,G15.7)') 'MINKD=',PBMAN%MINKD
  ENDIF
  IF(UTL_READINITFILE('MINC',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%MINC
   IF(IBATCH.EQ.1)WRITE(*,'(A,G15.7)') 'MINC=',PBMAN%MINC
  ENDIF
  IF(UTL_READINITFILE('IFVDL',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%IFVDL
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'IFVDL=',PBMAN%IFVDL
  ENDIF
 ENDIF
 IF(UTL_READINITFILE('DWEL',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%DWEL
  IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'DWEL=',PBMAN%DWEL
 ENDIF
 IF(UTL_READINITFILE('DISG',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%DISG
  IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'DISG=',PBMAN%DISG
 ENDIF
 IF(UTL_READINITFILE('DSFR',LINE,IU,1))THEN
  READ(LINE,*) PBMAN%DSFR
  IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'DSFR=',PBMAN%DSFR
 ENDIF
 IF(PBMAN%IFORMAT.EQ.2.OR.PBMAN%IFORMAT.EQ.3)THEN
  IF(UTL_READINITFILE('NWAVESUZF',LINE,IU,1))THEN
   READ(LINE,*) PBMAN%NWAVESUZF
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NWAVESUZF=',PBMAN%NWAVESUZF
  ENDIF
 ENDIF 
 SELECT CASE (PBMAN%IFORMAT)
  !## runfile does not support this
  CASE (1)
   PBMAN%INFFCT=0
  !## modflow6 always used inffct=1
  CASE (3)
   IF(UTL_READINITFILE('DATEFORMAT',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%DATEFORMAT; WRITE(*,'(A)') 'DATEFORMAT='//TRIM(VTOS(PBMAN%DATEFORMAT))
   ENDIF
   PBMAN%INFFCT=1
   IF(UTL_READINITFILE('IFILLHEAD',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%IFILLHEAD
    IF(IBATCH.EQ.1)WRITE(*,'(A,F10.2)') 'IFILLHEAD=',PBMAN%IFILLHEAD
   ENDIF
   IF(UTL_READINITFILE('MERGELAYERS',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%MERGELAYERS
    IF(IBATCH.EQ.1)WRITE(*,'(A,F10.2)') 'MERGELAYERS=',PBMAN%MERGELAYERS
   ENDIF
   IF(UTL_READINITFILE('DDRN',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%DDRN
    IF(IBATCH.EQ.1)WRITE(*,'(A,F10.2)') 'DDRN=',PBMAN%DDRN
   ENDIF
   IF(UTL_READINITFILE('QWEL',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%QWEL
    IF(IBATCH.EQ.1)WRITE(*,'(A,F10.2)') 'QWEL=',PBMAN%QWEL
   ENDIF
   IF(UTL_READINITFILE('DEPTHUZF',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%DEPTHUZF
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'DEPTHUZF='//TRIM(PBMAN%DEPTHUZF)
   ENDIF
   IF(UTL_READINITFILE('NWAVESUZF',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%NWAVESUZF
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NWAVESUZF=',PBMAN%NWAVESUZF
   ENDIF
   IF(UTL_READINITFILE('NEWTON',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%NEWTON
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'NEWTON=',PBMAN%NEWTON
    IF(PBMAN%NEWTON.EQ.1)PBMAN%COMPLEXITY='COMPLEX'
   ENDIF
   IF(UTL_READINITFILE('COUPLER',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%COUPLER
    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'COUPLER='//TRIM(PBMAN%COUPLER)
    IF(UTL_READINITFILE('TOPMODEL',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%TOPMODEL
     IF(IBATCH.EQ.1)WRITE(*,'(A)') 'TOPMODEL='//TRIM(VTOS(PBMAN%TOPMODEL))
    ENDIF
   ENDIF
   IF(PBMAN%NEWTON.EQ.0)THEN
    IF(UTL_READINITFILE('COMPLEXITY',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%COMPLEXITY
     SELECT CASE (PBMAN%COMPLEXITY)
      CASE ('SIMPLE','MODERATE','COMPLEX')
       IF(IBATCH.EQ.1)WRITE(*,'(A)') 'COMPLEXITY='//TRIM(PBMAN%COMPLEXITY)
      CASE DEFAULT
       IF(IBATCH.EQ.1)THEN; WRITE(*,'(A,F10.2)') 'COMPLEXITY NEED TO BE ONE FROM SIMPLE,MODERATE,COMPLEX'; RETURN; ENDIF     
     END SELECT
    ENDIF
   ELSE
    WRITE(*,'(A)') 'COMPLEXITY='//TRIM(PBMAN%COMPLEXITY)
   ENDIF
  !## imod-wq export to standard files
  CASE (6)
   PBMAN%INFFCT=1
  CASE DEFAULT
   IF(UTL_READINITFILE('INFFCT',LINE,IU,1))THEN
    READ(LINE,*) PBMAN%INFFCT
    IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'INFFCT=',PBMAN%INFFCT
   ENDIF
 END SELECT
 
 IF(UTL_READINITFILE('ISOLVE',LINE,IU,1))READ(LINE,*) PBMAN%ISOLVE
 IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'ISOLVE=',PBMAN%ISOLVE
 SELECT CASE (PBMAN%IFORMAT)
  !## mf2005-run,mf2005-nam
  CASE (1,2)
   I=1; IF(PBMAN%ISOLVE.EQ.1)I=0
   IF(.NOT.UTL_READINITFILE('MODFLOW',LINE,IU,I))RETURN
   READ(LINE,*) PBMAN%MODFLOW; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'MODFLOW='//TRIM(PBMAN%MODFLOW)
  !## modflow6
  CASE (3)
   IF(PBMAN%COUPLER.EQ.'')THEN
    I=1; IF(PBMAN%ISOLVE.EQ.1)I=0
    IF(.NOT.UTL_READINITFILE('MODFLOW6',LINE,IU,I))RETURN
    READ(LINE,*) PBMAN%MODFLOW6; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'MODFLOW6='//TRIM(PBMAN%MODFLOW6)
   ELSE
    IF(UTL_READINITFILE('MODFLOW6',LINE,IU,1))THEN
     READ(LINE,*) PBMAN%MODFLOW6
     IF(IBATCH.EQ.1)WRITE(*,'(A)') 'MODFLOW6='//TRIM(PBMAN%MODFLOW6)
    ENDIF
   ENDIF
  !## seawat,mt3d
  CASE (4,5,6)
   IF(.NOT.UTL_READINITFILE('IMOD-WQ',LINE,IU,0))RETURN
   READ(LINE,*) PBMAN%IMOD_WQ; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'IMOD-WQ='//TRIM(PBMAN%IMOD_WQ)
  !## modflow1988 (adjoint)
  CASE (7)
   I=1; IF(PBMAN%ISOLVE.EQ.1)I=0
   IF(.NOT.UTL_READINITFILE('MODFLOW_ADJ',LINE,IU,I))RETURN
   READ(LINE,*) PBMAN%MODFLOW_ADJ; IF(IBATCH.EQ.1)WRITE(*,'(A)') 'MODFLOW_ADJ='//TRIM(PBMAN%MODFLOW_ADJ)
 END SELECT

 IMODBATCH_RUNFILE_READ=.TRUE.

 END FUNCTION IMODBATCH_RUNFILE_READ
 
 !###======================================================================
 SUBROUTINE IMODBATCH_CREATEENSEMBLES_CHOLESKY(STDEV,MEAN,DIR,RANGE,NSIM,FNAME,ILOG)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(INOUT) :: STDEV,MEAN
 CHARACTER(LEN=*),INTENT(IN) :: DIR,FNAME
 REAL(KIND=DP_KIND),INTENT(IN) :: RANGE
 INTEGER,INTENT(IN) :: NSIM,ILOG
 INTEGER :: N,NN,MM,NE,I,II,J,JJ,IROW,ICOL,IROW2,ICOL2,JROW,JCOL,ISIM,KTYPE,NPOP,SEED,IERR,ITOC,ITIC,SVDOPT,IU,IL
 REAL(KIND=DP_KIND),DIMENSION(:,:),ALLOCATABLE :: COV,COR,EV,EVEC
 REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: V,M,X,E,EW,EVAL,WORK
 INTEGER,DIMENSION(:),ALLOCATABLE :: IWORK,IFAIL
 REAL(KIND=DP_KIND) :: GAMMA,XCOR,X1,Y1,X2,Y2,STDV,XT,SMALL,NUGGET,SILL,R,ECRIT,TOTEW,Z2,DXY
 TYPE(IDFOBJ) :: IDF
 LOGICAL :: LCHOL,LANI
 REAL(KIND=DP_KIND),DIMENSION(3) :: AGL,RAN
 REAL(KIND=DP_KIND),DIMENSION(3,3) :: ROT
 
 LCHOL=.FALSE.; IF(PEST%PE_CHOLESKY.EQ.1)LCHOL=.TRUE.
 
 CALL UTL_CREATEDIR(DIR)
 
 !## fill in variance (kriging)
 SILL=1.0D0; NUGGET=0.0D0; KTYPE=3 !exponential STOPS @ 1
 SMALL=0.01D0
 
 !## in log space sill cannot be 1.0
 
 N=STDEV%NROW*STDEV%NCOL; ALLOCATE(COV(N,N),COR(N,N),V(N),M(N)); COR=0.0D0; COV=0.0D0; V=0.0D0; M=0.0D0
 I=0; DO IROW=1,STDEV%NROW; DO ICOL=1,STDEV%NCOL
  IF(STDEV%X(ICOL,IROW).EQ.STDEV%NODATA.OR. &
     MEAN%X(ICOL ,IROW).EQ.MEAN%NODATA)THEN
   STDEV%X(ICOL,IROW)=STDEV%NODATA
   MEAN%X(ICOL ,IROW)=MEAN%NODATA
  ENDIF
  I=I+1; V(I)=STDEV%X(ICOL,IROW); M(I)=MEAN%X(ICOL,IROW)
 ENDDO; ENDDO
  
 COV=0.0D0
 DO I=1,N
  IF(V(I).NE.STDEV%NODATA)THEN
   COV(I,I)=V(I)*V(I)
   COR(I,I)=1.0D0
   COV(I,I)=COV(I,I)*COR(I,I)
  ENDIF
 ENDDO

 LANI=.FALSE.; DO II=1,SIZE(ELLIPS_IDF,1); DO JJ=1,SIZE(ELLIPS_IDF,2)
  IF(ELLIPS_IDF(II,JJ)%FNAME.NE.'')THEN
   IF(.NOT.ASSOCIATED(ELLIPS_IDF(II,JJ)%X))THEN
    CALL IDFCOPY(STDEV,ELLIPS_IDF(II,JJ))
    IF(.NOT.IDFALLOCATEX(ELLIPS_IDF(II,JJ)))RETURN
   ENDIF
   !## read ellipse information
   IF(.NOT.IDFREAD(ELLIPS_IDF(II,JJ),ELLIPS_IDF(II,JJ)%FNAME,1))RETURN; LANI=.TRUE.
  ENDIF
 ENDDO; ENDDO

 !## ani-isotropic (non-symmetric matrix gaat denk ik niet werken - of toch)
 IF(LANI)THEN 
 
  IROW=1; ICOL=0; DO I=1,N
   ICOL=ICOL+1; IF(ICOL.GT.STDEV%NCOL)THEN; ICOL=1; IROW=IROW+1; ENDIF

   !## skip nodata
   IF(COV(I,I).EQ.0.0D0)CYCLE !STDEV%NODATA)CYCLE
 
   CALL IDFGETLOC(STDEV,IROW,ICOL,X1,Y1)
   
   !## set up rotation matrix backwards
   AGL=0.0D0; DO II=1,SIZE(ELLIPS_IDF,1); DO JJ=1,SIZE(ELLIPS_IDF,2)
    CALL IDFIROWICOL(ELLIPS_IDF(II,JJ),JROW,JCOL,X1,Y1)
    IF(JROW.GT.0.AND.JROW.LE.ELLIPS_IDF(II,JJ)%NROW.AND. &
       JCOL.GT.0.AND.JCOL.LE.ELLIPS_IDF(II,JJ)%NCOL)THEN
     IF(ELLIPS_IDF(II,JJ)%X(JCOL,JROW).NE.ELLIPS_IDF(II,JJ)%NODATA)THEN
      IF(JJ.EQ.1)THEN
       AGL(II)=ELLIPS_IDF(II,JJ)%X(JCOL,JROW)
       !## project reverse on major x/y/z-axis
       AGL(II)=-1.0D0*AGL(II)
       !## convert to radians
       AGL(II)=AGL(II)/(360.0D0/(2.0D0*PI))
      ELSE
       RAN(II)=ELLIPS_IDF(II,JJ)%X(JCOL,JROW)
      ENDIF
     ENDIF
    ENDIF
   ENDDO; ENDDO
   R=MAXVAL(RAN)
 
   !## set up rotation matrice
   CALL UTL_SETUP_ROTATIONMATRIX(AGL,ROT)
 
   IROW2=1; ICOL2=0; DO J=1,N !I+1,N
    ICOL2=ICOL2+1; IF(ICOL2.GT.STDEV%NCOL)THEN; ICOL2=1; IROW2=IROW2+1; ENDIF
   
    !## skip nodata
    IF(COV(J,J).EQ.0.0D0)CYCLE !STDEV%NODATA)CYCLE
   
    CALL IDFGETLOC(STDEV,IROW2,ICOL2,X2,Y2)

    !## get distance
    Z2=0.0; DXY=KRIGING_DIST(X2,Y2,Z2,X1,Y1,0.0D0,0.0D0,IDF,0,0.0D0,R,-999,ROT,RAN)
    
    !## inside current range
    IF(DXY.LT.R)THEN 
 
     !## get the variogram values
     GAMMA=UTL_GETGAMMA(X1,Y1,X2,Y2,R,SILL,NUGGET,KTYPE)
   
     !## if log dna log10() ipv log() 
     !## if log also nugget might be 1 en sill 2?
     XCOR=1.0D0-GAMMA
     COR(I,J)=XCOR; COR(J,I)=XCOR
     COV(I,J)=SQRT(COV(I,I))*SQRT(COV(J,J))*XCOR
!     COV(J,I)=COV(I,J)
    
    ENDIF
    
   ENDDO 
  ENDDO 
 
 ELSE
 
  IROW=1; ICOL=0; DO I=1,N
   ICOL=ICOL+1; IF(ICOL.GT.STDEV%NCOL)THEN; ICOL=1; IROW=IROW+1; ENDIF

   !## skip nodata
   IF(COV(I,I).EQ.STDEV%NODATA)CYCLE
  
   CALL IDFGETLOC(STDEV,IROW,ICOL,X1,Y1)
   IROW2=IROW; ICOL2=ICOL; DO J=I+1,N
    ICOL2=ICOL2+1; IF(ICOL2.GT.STDEV%NCOL)THEN; ICOL2=1; IROW2=IROW2+1; ENDIF
   
    !## skip nodata
    IF(COV(J,J).EQ.STDEV%NODATA)CYCLE
   
    CALL IDFGETLOC(STDEV,IROW2,ICOL2,X2,Y2)
   
    !## get the variogram values
    GAMMA=UTL_GETGAMMA(X1,Y1,X2,Y2,RANGE,SILL,NUGGET,KTYPE)
   
    !## if log dna log10() ipv log() 
    !## if log also nugget might be 1 en sill 2?
    XCOR=1.0D0-GAMMA
    COR(I,J)=XCOR; COR(J,I)=XCOR
    COV(I,J)=SQRT(COV(I,I))*SQRT(COV(J,J))*XCOR
    COV(J,I)=COV(I,J)
  
   ENDDO
  ENDDO
 ENDIF
 
 !## save as IDF (for fun)
 CALL IDFNULLIFY(IDF)
 IDF%DX=1.0D0; IDF%DY=1.0D0; IDF%NCOL=N; IDF%NROW=N
 IDF%XMIN=0.0D0; IDF%XMAX=IDF%XMIN+IDF%NCOL*IDF%DX
 IDF%YMIN=0.0D0; IDF%YMAX=IDF%YMIN+IDF%NROW*IDF%DY
 IF(.NOT.IDFALLOCATEX(IDF))RETURN
 DO I=1,N; DO J=1,N; IDF%X(I,J)=COV(I,J); ENDDO; ENDDO
! IDF%FNAME=TRIM(DIR)//'\COV.IDF'; IF(.NOT.IDFWRITE(IDF,IDF%FNAME,1))STOP

 !## test symmetry
 DO I=1,N
  DO J=1,N
   IF(COV(I,J).NE.COV(J,I))THEN
    WRITE(*,*) COV(I,J),COV(J,I)
   ENDIF
  ENDDO
 ENDDO
 
 IF(LCHOL)THEN
 
  !## perform cholesky-decomposition A=LTL - delivers upper triangle is LT
  CALL CHOLESKYDECOMPOSITION(COV,N,1) 
 
  DO I=1,N; DO J=1,N; IDF%X(I,J)=COV(I,J); ENDDO; ENDDO
  IDF%FNAME=TRIM(DIR)//'\CHOL.IDF'; IF(.NOT.IDFWRITE(IDF,IDF%FNAME,1))STOP
  NN=N
  
 ELSE

  ALLOCATE(EV(N,N),E(N),EW(N))
  !## copy covariance matrix
  EV=COV
  
  ECRIT=99.0D0 !9.9D0
  
!  CALL LIS_EIGENVALUE(COV,LIS_EV,LIS_EW,ECRIT)
  CALL OSD_TIMER(ITOC)
  IF(.TRUE.)THEN
   SVDOPT=3
   SELECT CASE (SVDOPT)
    !## dsyev - compute them all
    CASE (1)
     ALLOCATE(WORK(1))
     CALL DSYEV('V','U',N,COV,N,EW,WORK,-1,IERR)
     NN=INT(WORK(1)); DEALLOCATE(WORK); ALLOCATE(WORK(NN))
     CALL DSYEV('V','U',N,COV,N,EW,WORK,NN,IERR)
     DEALLOCATE(WORK)
     EV=COV
     NE=N
    !## dsyevd went wrong (big array)
    CASE (2)
     ALLOCATE(WORK(1),IWORK(1))
     CALL DSYEVD('V','U',N,COV,N,EW,WORK,-1,IWORK,-1,IERR)
     NN=INT(WORK(1));  DEALLOCATE(WORK);  ALLOCATE(WORK(NN)); WRITE(*,*) NN
     MM=INT(IWORK(1)); DEALLOCATE(IWORK); ALLOCATE(IWORK(MM)); WRITE(*,*) MM
     CALL DSYEVD('V','U',N,COV,N,EW,WORK,NN,IWORK,MM,IERR)
     DEALLOCATE(WORK,IWORK)
     EV=COV
     NE=N
    !## dsyevx - # eigenvalues equal to number of ensembles
    CASE (3)

     !## get eigenvalues first
     ALLOCATE(WORK(1))
     EV=COV
     CALL DSYEV('N','U',N,COV,N,EW,WORK,-1,IERR)
     NN=INT(WORK(1)); DEALLOCATE(WORK); ALLOCATE(WORK(NN))
     CALL DSYEV('N','U',N,COV,N,EW,WORK,NN,IERR)
     DEALLOCATE(WORK)
     COV=EV
     
     TOTEW=SUM(EW)
     DO I=N,1,-1; EW(I)=100.0D0*EW(I)/TOTEW; IF(I.LT.N)EW(I)=EW(I+1)+EW(I); ENDDO
     DO I=N,1,-1; IF(EW(I).GT.ECRIT)EXIT; ENDDO
     IU=N
     IL=MAX(1,I)
     NE=IU-IL+1
     WRITE(*,'(/1X,A,F5.2,A1,I10/)') 'Number of remaining eigenvectors (',ECRIT,')',NE ! IU,IL,NE
     
     ALLOCATE(WORK(1),IWORK(1))
     DEALLOCATE(EW); ALLOCATE(EW(NE))
     DEALLOCATE(EV); ALLOCATE(EV(N,NE))
     ALLOCATE(IFAIL(N))
     CALL DSYEVX('V','I','U',N,COV,N,0.0D0,0.0D0,IL,IU,0.0D0,NE,EW,EV,N,WORK,-1,IWORK,IFAIL,IERR)
     NN=INT(WORK(1));  DEALLOCATE(WORK);  ALLOCATE(WORK(NN)) !;  WRITE(*,*) NN
     MM=5*N;           DEALLOCATE(IWORK); ALLOCATE(IWORK(MM))!; WRITE(*,*) MM
     !# compute eigenvectors
     CALL DSYEVX('V','I','U',N,COV,N,0.0D0,0.0D0,IL,IU,0.0D0,NE,EW,EV,N,WORK,NN,IWORK,IFAIL,IERR)
     DEALLOCATE(WORK,IWORK,IFAIL)
   END SELECT

   IF(IERR.NE.0)THEN
    WRITE(*,'(A)') 'LAPACK EIGENVALUE WENT WRONG'; STOP
   ENDIF
   
  else 
   !## perform eigenvalue decomposition
   CALL LUDCMP_TRED2(EV,N,N,EW,E)
   !## ev are the eigenvectors
   IF(.NOT.LUDCMP_TQLI(EW,E,N,N,EV))THEN; PAUSE; STOP; ENDIF
   NE=N
  endif
  CALL OSD_TIMER(ITIC)
  
  WRITE(*,*) ITIC-ITOC
  
  !## sort
  CALL LUDCMP_EIGSRT(EW,EV,N,NE)
  ALLOCATE(EVEC(NE,N))
  EVEC=TRANSPOSE(EV)
  DEALLOCATE(EV)
  
!  open(99,file='d:\tmp.txt',status='REPLACE')
!  do j=1,n; write(99,'(i10,99f15.7)') j,(ev(i,j),i=1,n); enddo
!  do j=1,n; write(99,'(i10,99f15.7)') j,ew(j); enddo
!  close(99)
  
  DO I=1,NE
   IF(EW(I).LE.1.0D-10)EW(I)=0.0D0
   EW(I)=SQRT(EW(I))
  ENDDO

  if(.TRUE.)then
  
  !## apply only percentage (99.9%) is toch nog wel veel ... :-(
  ALLOCATE(EVAL(N)); EVAL=EW
  TOTEW=SUM(EVAL)
  DO I=1,NE
   EVAL(I)=100.0D0*EVAL(I)/TOTEW
   IF(I.GT.1)EVAL(I)=EVAL(I-1)+EVAL(I)
  ENDDO
  DO I=1,NE
   IF(EVAL(I).GT.ECRIT)EXIT
  ENDDO
  
  NN=MIN(NE,I)
  DEALLOCATE(COV); ALLOCATE(COV(NN,N))
  DO I=1,NN
   DO J=1,N
    COV(I,J)=EVEC(I,J)*EW(I)
   ENDDO
  ENDDO
  DEALLOCATE(EVAL,EVEC,EW)
  
  else

  DEALLOCATE(COR); ALLOCATE(COR(N,N)); COR=0.0D0; DO I=1,N
   COR(I,I)=EW(I)
  ENDDO
  CALL UTL_MATMUL(EV,COR,COV)
  NN=N
  DEALLOCATE(EV,EW)
  
  endif

  
 ENDIF

 IF(ILOG.EQ.0)THEN
  WRITE(*,'(A10,2A15)') 'NPOP','MEAN','STDV'
 ELSEIF(ILOG.EQ.1)THEN
  WRITE(*,'(A10,4A15)') 'NPOP','LN(MEAN)','LN(STDV)','MEAN','STDV'
 ELSEIF(ILOG.EQ.2)THEN
  WRITE(*,'(A10,4A15)') 'NPOP','LOG10(MEAN)','LOG10(STDV)','MEAN','STDV'
 ENDIF

 ALLOCATE(X(NN)); X=0.0D0
 
 !## generate ensembles
 SEED=12345
 DO ISIM=1,NSIM

  DO I=1,NN
   !## generates number -3 to +3 - in log-space permeability is normal-distributed
   CALL IPEST_NORMAL_MS_SAMPLE(0.0D0,1.0D0,SEED,X(I))
  ENDDO

  V=0.0D0

  !## z=m*Lx
  DO I=1,N
   V(I)=0.0D0
   IF(LCHOL)THEN
    DO J=1,NN
     !## cholesky uses lower
     V(I)=V(I)+COV(I,J)*X(J)
    ENDDO
   ELSE
    !## eigenvalue decomposition multiplication eigenvalues with random number
    DO J=1,NN
     V(I)=V(I)+COV(J,I)*X(J)
    ENDDO
   ENDIF
   IF(M(I).NE.MEAN%NODATA)THEN
    V(I)=M(I)+V(I)
   ELSE
    V(I)=MEAN%NODATA
   ENDIF
  ENDDO

  !## all of the ensembles go to stdv - statistics is for log10 distribution
  CALL UTL_STDEF(V,N,-999.99D0,STDV,XT,NPOP)
  IF(ILOG.EQ.0)THEN
   WRITE(*,'(I10,2F15.7)') NPOP,XT,STDV
  ELSEIF(ILOG.EQ.1)THEN
   WRITE(*,'(I10,4F15.7)') NPOP,XT,STDV,EXP(XT),EXP(STDV)
  ELSEIF(ILOG.EQ.2)THEN
   WRITE(*,'(I10,4F15.7)') NPOP,XT,STDV,10.0**XT,10.0D0**STDV
  ENDIF
  
  DO I=1,N
   IF(V(I).NE.MEAN%NODATA)THEN
    IF(ILOG.EQ.1)V(I)=EXP(V(I))
    IF(ILOG.EQ.2)V(I)=10.0D0**V(I)
   ENDIF
  ENDDO

  I=0; DO IROW=1,STDEV%NROW; DO ICOL=1,STDEV%NCOL
   I=I+1
   MEAN%X(ICOL,IROW)=V(I)
  ENDDO; ENDDO
  MEAN%FNAME=TRIM(DIR)//'\'//TRIM(FNAME)//'_R'//TRIM(VTOS(ISIM))//'.IDF'
  IF(.NOT.IDFWRITE(MEAN,MEAN%FNAME,1))STOP
 ENDDO

 DEALLOCATE(COV,COR,X,V)

 END SUBROUTINE IMODBATCH_CREATEENSEMBLES_CHOLESKY
 
 END MODULE MOD_BATCH_UTL