MODULE MOD_TOOLS

 USE IMOD_UTL, ONLY : IMOD_UTL_CREATEDIR
 USE IMOD_IDF
 CHARACTER(LEN=256) :: FNAME,OUTPUTDIR
 CHARACTER(LEN=256),ALLOCATABLE,DIMENSION(:) :: FLIST
 REAL, PARAMETER :: EPS=1.0
 REAL, PARAMETER :: KERROR=1.0,CERROR=1.0
 REAL :: CMIN
 
 CHARACTER(LEN=256),PARAMETER :: REGISDIR='d:\rakotoni\work\lhm_multi_layers\REGIS_V2.2\IDF'
 CHARACTER(LEN=256),PARAMETER :: MODELDIR='d:\rakotoni\work\lhm_multi_layers\REGIS_V2.2\CONSISTENT_BLOCK'
 CHARACTER(LEN=256),PARAMETER :: WORKDIR ='d:\rakotoni\work\lhm_multi_layers'
 CHARACTER(LEN=256),PARAMETER :: OUTPTDIR='d:\rakotoni\work\lhm_multi_layers\LHM\LHM_AGGREGATION'

 CONTAINS

 !###===========================
 SUBROUTINE LHM_CONVERTREGIS_INIT()
 !###===========================
 IMPLICIT NONE
 INTEGER :: I,N,IU,IOS

 !## open Unit for Same Line Printing of echo (is equal to screen or '*')
 OPEN(UNIT=6,CARRIAGECONTROL='FORTRAN')    

 IU=GETUNIT(); OPEN(IU,FILE=TRIM(REGISDIR)//'\CHRONO_ORDER.TXT',STATUS='OLD',ACTION='READ')
 DO I=1,2; N=0; DO
   READ(IU,*,IOSTAT=IOS) FNAME; IF(IOS.NE.0)EXIT
   IF(TRIM(FNAME).EQ.'-')EXIT
   N=N+1; IF(I.EQ.2)FLIST(N)=TRIM(REGISDIR)//'\'//TRIM(FNAME)
  ENDDO; IF(I.EQ.1)ALLOCATE(FLIST(N))
  REWIND(IU)
 ENDDO; CLOSE(IU)
 
 END SUBROUTINE LHM_CONVERTREGIS_INIT

 
 !###===========================
 SUBROUTINE LHM_CONVERTREGIS()
 !###===========================
 IMPLICIT NONE
 INTEGER :: IU,IOS,I,J,N,IREC,ICOL,IROW
 REAL :: T,B,K1,K2
 TYPE(IDFOBJ),DIMENSION(:),ALLOCATABLE :: TP,BT,KH,KV
 TYPE(IDFOBJ) :: TPIDF,BTIDF,KHIDF,VAIDF,IBIDF
  
 CALL IMOD_UTL_CREATEDIR(MODELDIR)
 
 N=SIZE(FLIST); ALLOCATE(TP(N),BT(N),KH(N),KV(N))
 DO I=1,N; CALL IDFNULLIFY(TP(I)); CALL IDFNULLIFY(BT(I)); CALL IDFNULLIFY(KH(I)); CALL IDFNULLIFY(KV(I)); ENDDO
 CALL IDFNULLIFY(TPIDF); CALL IDFNULLIFY(BTIDF); CALL IDFNULLIFY(KHIDF); CALL IDFNULLIFY(VAIDF); CALL IDFNULLIFY(IBIDF)
 
 DO I=1,N
  WRITE(6,*) 'Reading '//TRIM(FLIST(I))//' ...'
  FNAME=TRIM(FLIST(I))//'-T-C.IDF'; IF(.NOT.IDFREAD(TP(I),FNAME,0))THEN; WRITE(6,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
  FNAME=TRIM(FLIST(I))//'-B-C.IDF'; IF(.NOT.IDFREAD(BT(I),FNAME,0))THEN; WRITE(6,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
  FNAME=TRIM(FLIST(I))//'-KH-S.IDF'; IF(.NOT.IDFREAD(KH(I),FNAME,0))THEN; ENDIF
  FNAME=TRIM(FLIST(I))//'-KV-S.IDF'; IF(.NOT.IDFREAD(KV(I),FNAME,0))THEN; ENDIF
 ENDDO

 DO I=2,N
  IF(.NOT.IDFEQUAL(TP(I),TP(1)))THEN; WRITE(*,'(/A/)') 'IDF '//TRIM(TP(I)%FNAME)//' NE TO '//TRIM(TP(1)%FNAME); STOP; ENDIF
  IF(.NOT.IDFEQUAL(BT(I),TP(1)))THEN; WRITE(*,'(/A/)') 'IDF '//TRIM(BT(I)%FNAME)//' NE TO '//TRIM(BT(1)%FNAME); STOP; ENDIF
  IF(KH(I)%NCOL.GT.0)THEN
   IF(.NOT.IDFEQUAL(KH(I),TP(1)))THEN; WRITE(*,'(/A/)') 'IDF '//TRIM(KH(I)%FNAME)//' NE TO '//TRIM(KH(1)%FNAME); STOP; ENDIF
  ENDIF
  IF(KV(I)%NCOL.GT.0)THEN
   IF(.NOT.IDFEQUAL(KV(I),TP(1)))THEN; WRITE(*,'(/A/)') 'IDF '//TRIM(KV(I)%FNAME)//' NE TO '//TRIM(KV(1)%FNAME); STOP; ENDIF
  ENDIF
 ENDDO

 CALL IDFCOPY(TP(1),TPIDF); CALL IDFCOPY(TP(1),BTIDF); CALL IDFCOPY(TP(1),KHIDF); CALL IDFCOPY(TP(1),VAIDF); CALL IDFCOPY(TP(1),IBIDF)
 IF(.NOT.IDFALLOCATEX(TPIDF))STOP; IF(.NOT.IDFALLOCATEX(BTIDF))STOP; IF(.NOT.IDFALLOCATEX(KHIDF))STOP; IF(.NOT.IDFALLOCATEX(VAIDF))STOP
 IF(.NOT.IDFALLOCATEX(IBIDF))STOP

 !## fill in values
 TPIDF%X=TP(1)%NODATA; BTIDF%X=BT(1)%NODATA; KHIDF%X=TP(1)%NODATA; VAIDF%X=TP(1)%NODATA; IBIDF%X=TP(1)%NODATA
 DO I=1,N
  WRITE(*,'(/A)') 'Processing '//TRIM(FLIST(I))
  IREC=ICF +10  +ABS(TPIDF%IEQ-1) *2    +TPIDF%IEQ*(TPIDF%NROW+TPIDF%NCOL) +TPIDF%ITB*2
  DO IROW=1,TP(1)%NROW; DO ICOL=1,TP(1)%NCOL
   IREC=IREC+1; IBIDF%X(ICOL,IROW)=0
   !## nodata remains nodata 
   IF(I.GT.1.AND.TPIDF%X(ICOL,IROW).EQ.TPIDF%NODATA)CYCLE
   READ(TP(I)%IU,REC=IREC) T; READ(BT(I)%IU,REC=IREC) B
   !## look for first no data
   IF(T.EQ.TP(I)%NODATA)THEN
    IF(I.EQ.1)THEN
     !## check below
     DO J=I+1,N
      READ(TP(J)%IU,REC=IREC) T
      IF(T.NE.TP(J)%NODATA)THEN
       B=T; EXIT
      ENDIF
     ENDDO
    ELSE
     T=BTIDF%X(ICOL,IROW); B=BTIDF%X(ICOL,IROW)
    ENDIF
    K1=TP(1)%NODATA; K2=TP(1)%NODATA
   ELSE
    K1=EPS; IF(KH(I)%IU.GT.0)READ(KH(I)%IU,REC=IREC) K1; IF(K1.EQ.KH(I)%NODATA)K1=EPS
    K2=EPS; IF(KV(I)%IU.GT.0)READ(KV(I)%IU,REC=IREC) K2; IF(K2.EQ.KV(I)%NODATA)K2=EPS
    IF(K1.EQ.0.0)K1=K2; IF(K2.EQ.0.0)K2=K1
    IF(K1.EQ.0.0.AND.K2.EQ.0.0)THEN
     K1=TP(1)%NODATA; K2=K1; T=BTIDF%X(ICOL,IROW); B=BTIDF%X(ICOL,IROW)
    ELSE
     !## make sure it connects
     IF(I.GT.1)T=BTIDF%X(ICOL,IROW)
    ENDIF
   ENDIF
   TPIDF%X(ICOL,IROW)=T; BTIDF%X(ICOL,IROW)=B; IBIDF%X(ICOL,IROW)=1
   IF(K1.NE.TP(1)%NODATA.AND.K2.NE.TP(1)%NODATA)THEN
    KHIDF%X(ICOL,IROW)=K1; VAIDF%X(ICOL,IROW)=K2/K1
   ELSE
    KHIDF%X(ICOL,IROW)=TP(1)%NODATA; VAIDF%X(ICOL,IROW)=1.0D0
   ENDIF
  ENDDO; WRITE(6,'(A,F10.3,A)') '+Progress ',REAL(100*IROW)/REAL(TP(1)%NROW),'%     '; ENDDO

  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-IB_L'//TRIM(ITOS(I))//'.IDF'; WRITE(*,*) TRIM(FNAME)
  IF(.NOT.IDFWRITE(IBIDF,FNAME,0,1))THEN; WRITE(*,'(/A/)') 'ERROR WRITING '//TRIM(FNAME); STOP; ENDIF  
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-T_L'//TRIM(ITOS(I))//'.IDF'; WRITE(*,*) TRIM(FNAME)
  IF(.NOT.IDFWRITE(TPIDF,FNAME,0,1))THEN; WRITE(*,'(/A/)') 'ERROR WRITING '//TRIM(FNAME); STOP; ENDIF  
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-KH_L'//TRIM(ITOS(I))//'.IDF'; WRITE(*,*) TRIM(FNAME)
  IF(.NOT.IDFWRITE(KHIDF,FNAME,0,1))THEN; WRITE(*,'(/A/)') 'ERROR WRITING '//TRIM(FNAME); STOP; ENDIF  
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-VA_L'//TRIM(ITOS(I))//'.IDF'; WRITE(*,*) TRIM(FNAME)
  IF(.NOT.IDFWRITE(VAIDF,FNAME,0,1))THEN; WRITE(*,'(/A/)') 'ERROR WRITING '//TRIM(FNAME); STOP; ENDIF  
  IF(I.EQ.N)THEN
   FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-B_L'//TRIM(ITOS(I))//'.IDF'; WRITE(*,*) TRIM(FNAME)
   IF(.NOT.IDFWRITE(BTIDF,FNAME,0,1))THEN; WRITE(*,'(/A/)') 'ERROR WRITING '//TRIM(FNAME); STOP; ENDIF  
  ENDIF

 ENDDO

 END SUBROUTINE LHM_CONVERTREGIS
 
 !###===========================
 SUBROUTINE LHM_CONVERTREGIS_AGGREGATE()
 !###===========================
 IMPLICIT NONE
 INTEGER :: IU,IOS,III,II,I,J,JJ,N,ICOL,IROW,BSIZE,BLOCKSIZE,NODES,INODE
 REAL :: KE,CE
 REAL,DIMENSION(2) :: TKDW,TVCW
 REAL,DIMENSION(:,:),ALLOCATABLE :: TOP,KHV,KVV,KVA,IBD
 REAL,DIMENSION(:),ALLOCATABLE :: VCV,KDW
 TYPE(IDFOBJ),DIMENSION(:,:),ALLOCATABLE :: TP,KH,VA,IB
  
 BLOCKSIZE=500000

 WRITE(OUTPUTDIR,'(A,I5.5)') TRIM(OUTPTDIR)//'\CMIN_',INT(CMIN)
 CALL IMOD_UTL_CREATEDIR(OUTPUTDIR)
 
 N=SIZE(FLIST); ALLOCATE(TP(N+1,2),KH(N,2),VA(N,2),IB(N,2))
 DO I=1,N+1; DO J=1,2; CALL IDFNULLIFY(TP(I,J)); ENDDO; ENDDO
 DO I=1,N;   DO J=1,2; CALL IDFNULLIFY(KH(I,J)); CALL IDFNULLIFY(VA(I,J)); ; CALL IDFNULLIFY(IB(I,J)); ENDDO; ENDDO

 DO I=1,N
  
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-T_L'//TRIM(ITOS(I))//'.IDF'
  IF(.NOT.IDFREAD(TP(I,1),FNAME,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF 
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-T_L'//TRIM(ITOS(I))//'.IDF'
  CALL IDFCOPY(TP(I,1),TP(I,2))
  IF(.NOT.IDFWRITE(TP(I,2),FNAME,0,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
  
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-KH_L'//TRIM(ITOS(I))//'.IDF'
  IF(.NOT.IDFREAD(KH(I,1),FNAME,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF  
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-KH_L'//TRIM(ITOS(I))//'.IDF'
  CALL IDFCOPY(KH(I,1),KH(I,2))
  IF(.NOT.IDFWRITE(KH(I,2),FNAME,0,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF

  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-VA_L'//TRIM(ITOS(I))//'.IDF'
  IF(.NOT.IDFREAD(VA(I,1),FNAME,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-VA_L'//TRIM(ITOS(I))//'.IDF'
  CALL IDFCOPY(VA(I,1),VA(I,2))
  IF(.NOT.IDFWRITE(VA(I,2),FNAME,0,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF

  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-IB_L'//TRIM(ITOS(I))//'.IDF'
  IF(.NOT.IDFREAD(IB(I,1),FNAME,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
  FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-IB_L'//TRIM(ITOS(I))//'.IDF'
  CALL IDFCOPY(IB(I,1),IB(I,2))
  IF(.NOT.IDFWRITE(IB(I,2),FNAME,0,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
 
  IF(I.EQ.N)THEN
   FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(MODELDIR)//'\'//TRIM(FNAME)//'-B_L'//TRIM(ITOS(I))//'.IDF'
   IF(.NOT.IDFREAD(TP(I+1,1),FNAME,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
   FNAME=FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:); FNAME=TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-B_L'//TRIM(ITOS(I))//'.IDF'
   CALL IDFCOPY(TP(I+1,1),TP(I+1,2))
   IF(.NOT.IDFWRITE(TP(I+1,2),FNAME,0,-1))THEN; WRITE(*,'(/A/)') 'ERROR READING '//TRIM(FNAME); STOP; ENDIF
  ENDIF
 
 ENDDO
 
 ALLOCATE(TOP(N+1,BLOCKSIZE),KHV(N,BLOCKSIZE),KVV(N,BLOCKSIZE),KVA(N,BLOCKSIZE),VCV(N),KDW(N), IBD(N,BLOCKSIZE))
 
 NODES=TP(1,1)%NROW*TP(1,1)%NCOL
 INODE=0
 DO 
  
  BSIZE=MIN(BLOCKSIZE,NODES-INODE)

  !## finished
  IF(BSIZE.LE.0)EXIT
  
  !## read variables at currrent location
  DO I=1,N
   READ(TP(I,1)%IU) (TOP(I,J),J=1,BSIZE)
   READ(KH(I,1)%IU) (KHV(I,J),J=1,BSIZE)
   READ(VA(I,1)%IU) (KVA(I,J),J=1,BSIZE)
   READ(IB(I,1)%IU) (IBD(I,J),J=1,BSIZE)
   DO J=1,BSIZE; KVV(I,J)=KHV(I,J)*KVA(I,J); ENDDO
  ENDDO
  READ(TP(N+1,1)%IU) (TOP(N+1,J),J=1,BSIZE)

  DO J=1,BSIZE

   !## compute c- and t-values
   DO I=1,N
    VCV(I)=(TOP(I,J)-TOP(I+1,J))/KVV(I,J)
    KDW(I)=(TOP(I,J)-TOP(I+1,J))*KHV(I,J)
   ENDDO
   TKDW(1)=SUM(KDW); TVCW(1)=SUM(VCV)
  
   !## skip this as it is all nodata
   IF(TKDW(1).LE.0.0)CYCLE
  
   !## aggregate from layer 1 onwards
   I=0; DO
    I=I+1
    DO II=I+1,N
     IF(VCV(II).GT.CMIN.OR.II.EQ.N)THEN
      !## shift bottoms down
      JJ=1; IF(II.EQ.N)JJ=0
      DO III=I+1,II-JJ !1
       TOP(III,J)=TOP(II,J)
       VCV(I)    =VCV(I)+VCV(III); VCV(III)=EPS!0.0
       KDW(I)    =KDW(I)+KDW(III); KDW(III)=EPS!0.0
       ! Make vertical pass through cell here
       IBD(III,J)=-1
      ENDDO
      EXIT
     ENDIF
    ENDDO
    I=II
    IF(I.GE.N)EXIT
   ENDDO
   
   TKDW(2)=SUM(KDW); TVCW(2)=SUM(VCV)
  
   KE=100.0*(TKDW(1)-TKDW(2))/TKDW(1)
   CE=100.0*(TVCW(1)-TVCW(2))/TVCW(1)
   IF(KE.GT.KERROR.OR.CE.GT.CERROR)THEN
    WRITE(*,'(/A)') 'Something went wrong'
    WRITE(*,*) TKDW(1),TKDW(2),TKDW(1)-TKDW(2),KE
    WRITE(*,*) TVCW(1),TVCW(2),TVCW(1)-TVCW(2),CE
   ENDIF
  
   !## compute khv and kvv-values
   DO I=1,N
    IF((TOP(I,J)-TOP(I+1,J)).GT.0.0)THEN
     KHV(I,J)=KDW(I)/(TOP(I,J)-TOP(I+1,J))
     KVV(I,J)=(TOP(I,J)-TOP(I+1,J))*VCV(I)
    ELSE
     KHV(I,J)=EPS
     KVV(I,J)=EPS
    ENDIF
   ENDDO  
  ENDDO
  
  !## save variables at currrent location
  DO I=1,N
   WRITE(TP(I,2)%IU) (TOP(I,J),J=1,BSIZE)
   WRITE(KH(I,2)%IU) (KHV(I,J),J=1,BSIZE)
   DO J=1,BSIZE; KVA(I,J)=KVV(I,J)/KHV(I,J); ENDDO
   WRITE(VA(I,2)%IU) (KVA(I,J),J=1,BSIZE)
   WRITE(IB(I,2)%IU) (IBD(I,J),J=1,BSIZE)
  ENDDO
  WRITE(TP(N+1,2)%IU) (TOP(N+1,J),J=1,BSIZE)
 
  INODE=INODE+BSIZE

  WRITE(6,'(A,F10.3,A)') '+Progress Aggregation ',REAL(100*INODE)/REAL(NODES),'%     '

 ENDDO
 
 DEALLOCATE(TOP,KHV,KVV,VCV,KDW)
 
 END SUBROUTINE LHM_CONVERTREGIS_AGGREGATE

 !###===========================
 LOGICAL FUNCTION IDF_EXIST(IDF)
 !###===========================
 IMPLICIT NONE
 TYPE(IDFOBJ), INTENT(IN) :: IDF
 LOGICAL :: LEX
 
 IDF_EXIST = .FALSE.
 INQUIRE(FILE=TRIM(IDF%FNAME), EXIST=LEX)
 IDF_EXIST = LEX
 
 END FUNCTION IDF_EXIST

 !###===========================
 SUBROUTINE LHM_CONVERTREGIS_WRITEPRJ()
 !###===========================
 IMPLICIT NONE
 INTEGER :: I,IU,N
 
 N=SIZE(FLIST)
 !CALL IMOD_UTL_CREATEDIR(WORKDIR)
 IU=IMOD_UTL_GETUNIT(); OPEN(IU,FILE=TRIM(OUTPUTDIR)//'\MODEL_v5.PRJ',STATUS='UNKNOWN',ACTION='WRITE')

 WRITE(IU,'(/A)') '0001,(BND),1, BOUNDARY CONDITION,[BND]'
 WRITE(IU,'(A,I10)') '001,',N
 DO I = 1,N
  FNAME = FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:)
  FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-IB_L'//TRIM(ITOS(I))//'.IDF'
  WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
 ENDDO
 
 WRITE(IU,'(/A)') '0001,(TOP),1, TOP ELEVATION,[TOP]'
 WRITE(IU,'(A,I10)') '001,',N
 DO I = 1,N
  FNAME = FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:)
  FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-T_L'//TRIM(ITOS(I))//'.IDF'
  WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
 ENDDO
 
 WRITE(IU,'(/A)') '0001,(BOT),1, BOTTOM ELEVATION,[BOT]'
 WRITE(IU,'(A,I10)') '001,',N
 DO I = 1,N
  IF(I.EQ.N)THEN
   FNAME = FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:)
   FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-B_L'//TRIM(ITOS(I))//'.IDF'
   WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
  ELSE
   FNAME = FLIST(I+1)(INDEX(FLIST(I+1),'\',.TRUE.)+1:)
   FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-T_L'//TRIM(ITOS(I+1))//'.IDF'
   WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
  ENDIF
 ENDDO

 WRITE(IU,'(/A)') '0001,(KHV),1, HORIZONTAL PERMEABILITY,[KHV]'
 WRITE(IU,'(A,I10)') '001,',N
 DO I = 1,N
  FNAME = FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:)
  FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-KH_L'//TRIM(ITOS(I))//'.IDF'
  WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
 ENDDO

 WRITE(IU,'(/A)') '0001,(KVA),1, VERTICAL ANISOTROPY,[KVA]'
 WRITE(IU,'(A,I10)') '001,',N
 
 DO I = 1,N
  FNAME = FLIST(I)(INDEX(FLIST(I),'\',.TRUE.)+1:)
  FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-VA_L'//TRIM(ITOS(I))//'.IDF'
  WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
 ENDDO
 
 WRITE(IU,'(/A)') '0001,(KVV),1, VERTICAL PERMEABILITY,[KVV]'
 WRITE(IU,'(A,I10)') '001,',N
 
 DO I = 1,N
  WRITE(IU,'(A,I3,A)') '1,1,',I,',1.0,0.0,1.0,""'
 ENDDO
 
 WRITE(IU,'(/A)') '0001 ,(SHD),1, STARTING HEADS,[SHD]'
 WRITE(IU,'(A,I10)') '001,',N
 
 DO I = 1,N
  FNAME = FLIST(1)(INDEX(FLIST(1),'\',.TRUE.)+1:)
  FNAME = TRIM(OUTPUTDIR)//'\'//TRIM(FNAME)//'-T_L'//TRIM(ITOS(1))//'.IDF'
  WRITE(IU,'(A,I3,A)') '1,2,',I,',1.0,0.0,-999.99,'//TRIM(FNAME)
 ENDDO

 WRITE(IU,'(/A)') '0001,(RIV),1, RIVER'
 WRITE(IU,'(A)') 'STEADY-STATE'
 WRITE(IU,'(A)') '004,005'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\COND_HL1_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\COND_HL2_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\COND_P_L0.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\COND_S_L0.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\COND_T_L0.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\PEIL_HW_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\PEIL_HW_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\PEIL_P1W_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\PEIL_S1W_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\PEIL_T1W_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\BOTH_W_L1.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\BOTH_W_L2.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\BODH_P1W_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\BODH_S1W_250.IDF' 
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\PEIL_T1W_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\INFMZ_H_250_L1.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\INFMZ_H_250_L2.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\INFMZ_P1_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\INFMZ_S1_250.IDF'
 WRITE(IU,'(A)') '1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\INFMZ_T1_250.IDF'
 
 WRITE(IU,'(/A)') '0001,(DRN),1, DRAINAGE'
 WRITE(IU,'(A)') 'STEADY-STATE'
 WRITE(IU,'(A)') '002,001'
 WRITE(IU,'(A)') ' 1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\COND_B_250.IDF'
 WRITE(IU,'(A)') ' 1,2, -001,   1.000000    ,   0.000000    ,   8.000000    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\BODH_B_250.IDF'
 
 WRITE(IU,'(/A)') '0001,(GHB),1, GENERAL HEAD BOUNDARY,[CON,LVL]'
 WRITE(IU,'(A)') 'STEADY-STATE'
 WRITE(IU,'(A)') '002,001'
 WRITE(IU,'(A)') ' 1,2, -001,   1.000000    ,   0.000000    ,  -999.9900    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\GHB_COND_REGIS_L01.IDF'
 WRITE(IU,'(A)') ' 1,2, -001,   1.000000    ,   0.000000    ,   8.000000    ,'//TRIM(WORKDIR)//'\LHM\LHMV1\GHB_STAGE_REGIS_L01.IDF'
 
 WRITE(IU,'(/A)') '0001,(RCH),1, RECHARGE,[RCH]'
 WRITE(IU,'(A)') 'STEADY-STATE'
 WRITE(IU,'(A)') '001,001'
 WRITE(IU,'(A)') ' 1,1, -001,   1.000000    ,   0.000000    ,   0.700000    ,""'
 
 WRITE(IU,'(/A)') ' 0001,(PCG),1, PRECONDITION CONJUGATE-GRADIENT []'
 WRITE(IU,'(A)')  '50,150,0.10000E-02,10.00000,0.98000,1,0,1,1.0000,1.0000,1,5.00000'
 
 WRITE(IU,'(/A)') 'PERIODS'

 END SUBROUTINE LHM_CONVERTREGIS_WRITEPRJ
 
 
END MODULE MOD_TOOLS

!###===========================
PROGRAM TOOLS
!###===========================
USE MOD_TOOLS
IMPLICIT NONE
CHARACTER(100) :: STRING

IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN; WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENTS REQUIRED, STOPPING'; STOP; ENDIF
CALL GET_COMMAND_ARGUMENT(1,STRING)
READ(STRING,*) CMIN

!## initialize
CALL LHM_CONVERTREGIS_INIT()
!## lhm
!CALL LHM_CONVERTREGIS()
!## aggregate
!CALL LHM_CONVERTREGIS_AGGREGATE()
!!## write prj file
!CALL LHM_CONVERTREGIS_WRITEPRJ()

STOP

END PROGRAM TOOLS