!!  Copyright (C) Stichting Deltares, 2005-2014.
!!
!!  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_SCENTOOL_WELLS

USE WINTERACTER
USE RESOURCE
USE MODPLOT, ONLY : MPW
USE MOD_SCENTOOL_PAR
USE MOD_SCENTOOL_UTL, ONLY : ST_DRAWPNTS,ST1FIELDS,ST_SYMBOLDRAW,ST1_PROPUPDATEGRID,ST1_PROPFIELDS, &
                    ST1_GETCOLOUR,ST1_PROPGRAPH,ST1_PROPOPENSAVE
USE MOD_UTL, ONLY : UTL_HIDESHOWDIALOG,JDATETOGDATE,ITOS,RTOS,GDATETOJDATE,UTL_GETUNIT,UTL_JDATETOIDATE, &
                    UTL_CREATEDIR,UTL_IDATETOJDATE,UTL_GETMED
USE MOD_IDF, ONLY : IDFIROWICOL,IDFGETVAL,IDFREAD,IDFDEALLOCATE
USE MOD_OSD, ONLY : OSD_OPEN
USE MOD_IPF, ONLY : IPFPLOTLABEL,IPFREAD2,IPF,IPFWRITE
USE MOD_PMANAGER, ONLY : PMANAGER_SAVEMF2005_PCK_GETTLP,PMANAGER_SAVEMF2005_PCK_READTXT

INTEGER,PRIVATE :: ITIME  !## fill in dates (1) fill in durations (2)
CHARACTER(LEN=10),PRIVATE :: CDATE
INTEGER,PRIVATE :: IMODE

CONTAINS

 !###======================================================================
 SUBROUTINE STWEL1INIT(IOPT)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IOPT

 IMODE=IOPT

 IF(IOPT.EQ.ID_ADD)THEN
  IF(NWEL.EQ.MAXNWEL)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You can not add another WEL SYSTEM in the current configuration.'//CHAR(13)// &
     'Currently the maximum WEL SYSTEMS is '//TRIM(ITOS(MAXNWEL))//'.'//CHAR(13)// &
     'Increase/add the keyword MAXNWEL= in the SCENTOOL file to increase the number of WEL SYSTEMS allowed','Information')
   IMODE=0
   CALL STWEL1CLOSE()
   IMODE=IOPT
  ENDIF
  NWEL           =NWEL+1
  IWEL           =NWEL
  WEL(IWEL)%CNAME='Well Field '//TRIM(ITOS(IWEL))
  WEL(IWEL)%NLOC =0
  WEL(IWEL)%NQ   =0
  WEL(IWEL)%ILOCT=1
  WEL(IWEL)%ICLR =WRGB(255,0,0)
  WEL(IWEL)%ISYMBOL=28
  ALLOCATE(WEL(IWEL)%Q(NROWQ))
  ALLOCATE(WEL(IWEL)%LOC(NROWL))
 ELSE
  CALL WDIALOGSELECT(ID_DSCENTOOLTAB1)
  CALL WDIALOGGETMENU(IDF_MENU1,IWEL)
 ENDIF

 CALL UTL_HIDESHOWDIALOG(ID_DSCENTOOL,0)

 CALL WINDOWSELECT(0)
 CALL WMENUSETSTATE(ID_SCENTOOL,1,0)

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB1)
 CALL WDIALOGPUTSTRING(IDF_LABEL1,'Well System name:')
 CALL WDIALOGPUTRADIOBUTTON(IDF_RADIO2) !## start enter durations

 !## put grid columnlabels
 CALL WGRIDLABELCOLUMN(IDF_GRID1,1,'Start Date')
 CALL WGRIDLABELCOLUMN(IDF_GRID1,2,'Duration')
 CALL WGRIDLABELCOLUMN(IDF_GRID1,3,'Rate (m3/hr)')

 CALL STWEL1PUTFIELDS()
 CALL ST1_PROPUPDATEGRID(1)
 CALL ST1_PROPFIELDS()

 !## initialize date in grid
 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB1)
 IF(WEL(IWEL)%NQ.EQ.0)THEN
  CALL WDIALOGGETSTRING(IDF_STRING2,CDATE)
  CALL WGRIDPUTCELLSTRING(IDF_GRID1,1,1,CDATE)
 ENDIF

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROP)
 CALL WDIALOGSETTAB(IDF_TAB,ID_DSCENTOOL_PROPTAB1)
 CALL WDIALOGTITLE('Well Systems')
 CALL WDIALOGSHOW(-1,-1,0,2)

 END SUBROUTINE STWEL1INIT

 !###======================================================================
 SUBROUTINE STWEL1MAIN(ITYPE,MESSAGE)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ITYPE
 TYPE(WIN_MESSAGE),INTENT(IN) :: MESSAGE

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROP)

 SELECT CASE (MESSAGE%WIN)

  !## main dialog
  CASE (ID_DSCENTOOL_PROP)
   SELECT CASE (ITYPE)
    CASE (PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      !## quit well editing
      CASE (IDCANCEL)
       CALL STWEL1CLOSE()
     END SELECT
    CASE (TABCHANGED)
     SELECT CASE (MESSAGE%VALUE2)
      CASE (ID_DSCENTOOL_PROPTAB3)
       CALL ST_SYMBOLDRAW(ID_DSCENTOOL_PROPTAB3,IDF_PICTURE1,IDF_MENU1,IDF_INTEGER1)
     END SELECT
   END SELECT

  !## well strengths
  CASE (ID_DSCENTOOL_PROPTAB1)
   SELECT CASE (ITYPE)
    CASE (FIELDCHANGED)
     CALL ST1_PROPFIELDS()
    CASE (PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      CASE (ID_OPEN,ID_SAVEAS)
       CALL ST1_PROPOPENSAVE(MESSAGE%VALUE1,ID_DSCENTOOL_PROPTAB1,1)
      CASE (ID_GRAPH)
       CALL ST1_PROPGRAPH(1)
      CASE (ID_CALC)
       CALL ST1_PROPUPDATEGRID(1)
     END SELECT
   END SELECT

  !## well locations
  CASE (ID_DSCENTOOL_PROPTAB2)
   SELECT CASE (ITYPE)
    CASE (FIELDCHANGED)
!      CALL ST1_PROPFIELDS()()
    CASE (PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      CASE (ID_OPEN,ID_SAVEAS)
       CALL ST1_PROPOPENSAVE(MESSAGE%VALUE1,ID_DSCENTOOL_PROPTAB2,1)
      CASE (ID_MOVE,ID_DRAW,ID_DELETE)
       CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB2)
       CALL WDIALOGFIELDSTATE(ID_OPEN,2)
       CALL WDIALOGFIELDSTATE(ID_SAVEAS,2)
       IF(MESSAGE%VALUE1.EQ.ID_MOVE)THEN
        CALL WDIALOGFIELDSTATE(ID_DRAW,2)
        CALL WDIALOGFIELDSTATE(ID_DELETE,2)
        CALL ST_DRAWPNTS(ID_DSCENTOOL_PROPTAB2,IDF_GRID1,(/4,5/),2,WEL(IWEL)%ISYMBOL,WEL(IWEL)%ICLR,1)  !## move
       ELSEIF(MESSAGE%VALUE1.EQ.ID_DELETE)THEN
        CALL WDIALOGFIELDSTATE(ID_DRAW,2)
        CALL WDIALOGFIELDSTATE(ID_MOVE,2)
        CALL ST_DRAWPNTS(ID_DSCENTOOL_PROPTAB2,IDF_GRID1,(/4,5/),3,WEL(IWEL)%ISYMBOL,WEL(IWEL)%ICLR,1)  !## delete
       ELSE
        CALL WDIALOGFIELDSTATE(ID_MOVE,2)
        CALL WDIALOGFIELDSTATE(ID_DELETE,2)
        CALL ST_DRAWPNTS(ID_DSCENTOOL_PROPTAB2,IDF_GRID1,(/4,5/),1,WEL(IWEL)%ISYMBOL,WEL(IWEL)%ICLR,1)  !## add
       ENDIF
       CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB2)
       CALL WDIALOGFIELDSTATE(ID_OPEN,1)
       CALL WDIALOGFIELDSTATE(ID_SAVEAS,1)
       CALL WDIALOGFIELDSTATE(ID_MOVE,1)
       CALL WDIALOGFIELDSTATE(ID_DRAW,1)
       CALL WDIALOGFIELDSTATE(ID_DELETE,1)
     END SELECT
   END SELECT

  !## well display options
  CASE (ID_DSCENTOOL_PROPTAB3)
   SELECT CASE (ITYPE)
    CASE (RESIZE,EXPOSE)
     CALL ST_SYMBOLDRAW(ID_DSCENTOOL_PROPTAB3,IDF_PICTURE1,IDF_MENU1,IDF_INTEGER1)
    CASE (FIELDCHANGED)
     SELECT CASE (MESSAGE%VALUE2)
      CASE (IDF_MENU1,IDF_REAL1)
       CALL ST_SYMBOLDRAW(ID_DSCENTOOL_PROPTAB3,IDF_PICTURE1,IDF_MENU1,IDF_INTEGER1)
      END SELECT   !IDD,IDF,IDM,IDC,IDS)
    CASE(PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      CASE (ID_COLOUR)
       CALL ST1_GETCOLOUR(ID_DSCENTOOL_PROPTAB3,IDF_INTEGER1)
       CALL ST_SYMBOLDRAW(ID_DSCENTOOL_PROPTAB3,IDF_PICTURE1,IDF_MENU1,IDF_INTEGER1)
     END SELECT
  END SELECT
 END SELECT

 END SUBROUTINE STWEL1MAIN

 !###======================================================================
 SUBROUTINE STWEL1CLOSE()
 !###======================================================================
 IMPLICIT NONE

 IF(IMODE.NE.0)THEN
  CALL WMESSAGEBOX(YESNOCANCEL,QUESTIONICON,COMMONNO,'Do you want to save any adjustments ?','Question')
  !## save and quit
  IF(WINFODIALOG(4).EQ.1)THEN
   CALL ST1_PROPUPDATEGRID(1)         !## recompute grid and sort data
   IF(.NOT.STWEL1GETFIELDS())RETURN   !## get data from grid, if not correct, return - not closing!
  !## no save and quit
  ELSEIF(WINFODIALOG(4).EQ.2)THEN
   IF(IMODE.EQ.ID_ADD)THEN
    IWEL=IWEL-1
    NWEL=NWEL-1
   ENDIF
  ENDIF
  !## not canceling
  IF(WINFODIALOG(4).EQ.0)RETURN
 ENDIF

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROP)
 CALL WDIALOGHIDE()

 !## change tab-header for number of well systems available
 CALL STWEL1TITLE()

 IWEL=0

 CALL UTL_HIDESHOWDIALOG(ID_DSCENTOOL,2)
 CALL WMENUSETSTATE(ID_SCENTOOL,1,1)
 CALL ST1FIELDS()
 CALL IDFPLOT(1)

 END SUBROUTINE STWEL1CLOSE

 !###======================================================================
 SUBROUTINE STWEL1TITLE()
 !###======================================================================
 IMPLICIT NONE

 CALL WDIALOGSELECT(ID_DSCENTOOLTAB1)
 CALL WDIALOGTITLE('Well Systems ('//TRIM(ITOS(NWEL))//')')
 IF(IWEL.NE.0)CALL WDIALOGPUTMENU(IDF_MENU1,WEL%CNAME,NWEL,IWEL)

 END SUBROUTINE STWEL1TITLE

 !###======================================================================
 SUBROUTINE ST1SAVELOADWELLS(IU,ICODE)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IU,ICODE
 CHARACTER(LEN=256) :: LINE
 INTEGER :: I,J,NDAY,IDATE

 IF(ICODE.EQ.ID_SAVEAS.OR.ICODE.EQ.ID_SAVE)THEN

  LINE=TRIM(ITOS(NWEL))
  WRITE(IU,'(A)') TRIM(LINE)//' !## nwells[i]'
  DO I=1,NWEL
   WRITE(IU,'(A)') '"'//TRIM(WEL(I)%CNAME)//'" !## name of well system[c]'
   LINE=TRIM(ITOS(WEL(I)%NQ-1))//','//TRIM(ITOS(WEL(I)%ITYPE))//','// &
      TRIM(ITOS(UTL_JDATETOIDATE(WEL(I)%Q(1)%IDATE)))
   WRITE(IU,'(A)') TRIM(LINE)//' !## nrates[i],dtype[i],sdate[c]'
   DO J=2,WEL(I)%NQ
    LINE=TRIM(ITOS( WEL(I)%Q(J)%IDATE-WEL(I)%Q(J-1)%IDATE ))//','//TRIM(RTOS(WEL(I)%Q(J-1)%QRATE/24.0,'F',2))  !## m3/day -> m3/hr
    WRITE(IU,'(A)') TRIM(LINE)//' !## delt[i],qrate[r]'
   END DO
   LINE=TRIM(ITOS(WEL(I)%NLOC))//','//TRIM(ITOS(WEL(I)%ILOCT))//','//  &
        TRIM(ITOS(WEL(I)%ISYMBOL))//','//TRIM(ITOS(WEL(I)%ICLR)) !## number of locations, iloct:1=msl 2=surface
   WRITE(IU,'(A)') TRIM(LINE)//' !## nloc[i],ltype[i],isymbol[i],iclr[i]'
   DO J=1,WEL(I)%NLOC
    LINE='"'//TRIM(WEL(I)%LOC(J)%ID)//'",'//TRIM(RTOS(WEL(I)%LOC(J)%Z1,'F',2))//','// &
                                            TRIM(RTOS(WEL(I)%LOC(J)%Z2,'F',2))//','// &
                                            TRIM(RTOS(WEL(I)%LOC(J)%X ,'F',2))//','// &
                                            TRIM(RTOS(WEL(I)%LOC(J)%Y ,'F',2))
    WRITE(IU,'(A)') TRIM(LINE)//' !## id[c],z1[r],z2[r],x[r],y[r]'
   END DO
  END DO

 ELSEIF(ICODE.EQ.ID_OPEN)THEN

  READ(IU,*) NWEL
  DO I=1,NWEL
   READ(IU,*) WEL(I)%CNAME
   READ(IU,*) WEL(I)%NQ,WEL(I)%ITYPE,IDATE
   WEL(I)%NQ=WEL(I)%NQ+1
   IF(.NOT.ASSOCIATED(WEL(I)%Q)) ALLOCATE(WEL(I)%Q(NROWQ))
   WEL(I)%Q(1)%IDATE=UTL_IDATETOJDATE(IDATE)
   DO J=2,WEL(I)%NQ
    READ(IU,*) NDAY,WEL(I)%Q(J-1)%QRATE
    WEL(I)%Q(J)%IDATE  =WEL(I)%Q(J-1)%IDATE+NDAY
    WEL(I)%Q(J-1)%QRATE=WEL(I)%Q(J-1)%QRATE*24.0  !## m3/hr -> m3/day
   END DO
   READ(IU,*) WEL(I)%NLOC,WEL(I)%ILOCT,WEL(I)%ISYMBOL,WEL(I)%ICLR
   IF(.NOT.ASSOCIATED(WEL(I)%LOC))ALLOCATE(WEL(I)%LOC(NROWL))
   DO J=1,WEL(I)%NLOC
    READ(IU,*) WEL(I)%LOC(J)%ID,WEL(I)%LOC(J)%Z1,WEL(I)%LOC(J)%Z2,WEL(I)%LOC(J)%X,WEL(I)%LOC(J)%Y
   ENDDO
  END DO

 ENDIF

 END SUBROUTINE ST1SAVELOADWELLS

 !###======================================================================
 SUBROUTINE ST1PLOTWELLS()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,IROW,ICLR,IMARKER
 REAL,PARAMETER :: TSIZE=0.01
 REAL :: X,Y,Z,TWIDTH,THEIGTH
 CHARACTER(LEN=50),DIMENSION(5) :: CLABEL,CVALUE
 CHARACTER(LEN=10),DIMENSION(2) :: CUNIT
 INTEGER,DIMENSION(5) :: IOPTION
 INTEGER :: IALL,JWEL,ILABEL
 DATA CUNIT/'m+MSL','m+SLevel'/

 CALL WDIALOGSELECT(ID_DSCENTOOLTAB1)
 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,IALL)
 CALL WDIALOGGETCHECKBOX(IDF_CHECK2,ILABEL)
 CALL WDIALOGGETMENU(IDF_MENU1,JWEL)

 CLABEL=''
 IF(ILABEL.EQ.1)THEN
  CLABEL(1)='ID'
  CLABEL(2)='Top Screen'
  CLABEL(3)='Bot Screen'
  CLABEL(4)='X-coord.'
  CLABEL(5)='Y-coord.'
 ENDIF

 !## overrule in case iwel.ne.0
 IF(IWEL.NE.0)THEN
  IALL=0
  JWEL=0
 ENDIF

 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO1,I)
 IOPTION=1
 SELECT CASE (I)
  CASE (1) !## no labeling
   IOPTION=0
  CASE (2) !## id
   IOPTION(2:5)=0
  CASE (3) !## id+filter screens
   IOPTION(4:5)=0
  CASE (4) !## all
   IOPTION=1
 END SELECT

 DO I=1,NWEL

  TWIDTH =TSIZE/2.0
  THEIGTH=TWIDTH*(0.03333333/0.013333)
  THEIGTH=THEIGTH*WINFOGRREAL(GRAPHICSRATIO)

  !## currently selected well-system - draw from grid-field
  IF(I.EQ.IWEL)THEN
   CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB3)
   CALL WDIALOGGETINTEGER(IDF_INTEGER1,ICLR)
   CALL WDIALOGGETMENU(IDF_MENU1,IMARKER)
   CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB2)
   DO IROW=1,NROWL
    X=0.0
    Y=0.0
    CVALUE=''
    IF(WINFOGRIDCELL(IDF_GRID1,1,IROW,GRIDCELLDEFINED).EQ.1)CALL WGRIDGETCELLSTRING(IDF_GRID1,1,IROW,CVALUE(1)) !## id
    IF(WINFOGRIDCELL(IDF_GRID1,2,IROW,GRIDCELLDEFINED).EQ.1)THEN
     CALL WGRIDGETCELLREAL(IDF_GRID1,2,IROW,Z)   !## z1
     CVALUE(2)=TRIM(RTOS(Z,'F',2))//' '//TRIM(CUNIT(WEL(I)%ILOCT))
    ENDIF
    IF(WINFOGRIDCELL(IDF_GRID1,3,IROW,GRIDCELLDEFINED).EQ.1)THEN
     CALL WGRIDGETCELLREAL(IDF_GRID1,3,IROW,Z)   !## z2
     CVALUE(3)=TRIM(RTOS(Z,'F',2))//' '//TRIM(CUNIT(WEL(I)%ILOCT))
    ENDIF
    IF(WINFOGRIDCELL(IDF_GRID1,4,IROW,GRIDCELLDEFINED).EQ.1)THEN
     CALL WGRIDGETCELLREAL(IDF_GRID1,4,IROW,X) !## x
     CVALUE(4)=TRIM(RTOS(X,'F',2))//' m'
    ENDIF
    IF(WINFOGRIDCELL(IDF_GRID1,5,IROW,GRIDCELLDEFINED).EQ.1)THEN
     CALL WGRIDGETCELLREAL(IDF_GRID1,5,IROW,Y) !## y
     CVALUE(5)=TRIM(RTOS(Y,'F',2))//' m'
    ENDIF
    CALL WGRTEXTFONT(WIDTH=TSIZE,HEIGHT=TSIZE*WINFOGRREAL(GRAPHICSRATIO))
    CALL IGRCOLOURN(ICLR)       !WEL(I)%ICLR)
    CALL IGRMARKER(X,Y,IMARKER) !WEL(I)%ISYMBOL)
    CALL IPFPLOTLABEL(X,Y,CVALUE,IOPTION,5,TWIDTH,THEIGTH,CLABEL,.FALSE.,-1)
   ENDDO
  ELSE
   DO IROW=1,WEL(I)%NLOC
    IF(IALL.EQ.1.OR.(IALL.EQ.0.AND.JWEL.EQ.I))THEN
     CVALUE(1)=WEL(I)%LOC(IROW)%ID
     CVALUE(2)=TRIM(RTOS(WEL(I)%LOC(IROW)%Z1,'F',2))//' '//TRIM(CUNIT(WEL(I)%ILOCT))
     CVALUE(3)=TRIM(RTOS(WEL(I)%LOC(IROW)%Z2,'F',2))//' '//TRIM(CUNIT(WEL(I)%ILOCT))
     CVALUE(4)=TRIM(RTOS(WEL(I)%LOC(IROW)%X ,'F',2))//' m'
     CVALUE(5)=TRIM(RTOS(WEL(I)%LOC(IROW)%Y ,'F',2))//' m'
     CALL IPFPLOTLABEL(WEL(I)%LOC(IROW)%X,WEL(I)%LOC(IROW)%Y,CVALUE,IOPTION,5,TWIDTH,THEIGTH,CLABEL,.FALSE.,-1)
    ENDIF
    CALL IGRCOLOURN(WEL(I)%ICLR)
    CALL WGRTEXTFONT(WIDTH=TSIZE,HEIGHT=TSIZE*WINFOGRREAL(GRAPHICSRATIO))
    CALL IGRMARKER(WEL(I)%LOC(IROW)%X,WEL(I)%LOC(IROW)%Y,WEL(I)%ISYMBOL)
   ENDDO
  ENDIF
 END DO

 END SUBROUTINE ST1PLOTWELLS

 !###======================================================================
 SUBROUTINE ST1SIMBOXWELLS(XMIN,YMIN,XMAX,YMAX)
 !###======================================================================
 IMPLICIT NONE
 REAL,INTENT(INOUT) :: XMIN,YMIN,XMAX,YMAX
 INTEGER :: I,J

 DO I=1,NWEL
  DO J=1,WEL(I)%NLOC
   XMIN=MIN(XMIN,WEL(I)%LOC(J)%X)
   XMAX=MAX(XMAX,WEL(I)%LOC(J)%X)
   YMIN=MIN(YMIN,WEL(I)%LOC(J)%Y)
   YMAX=MAX(YMAX,WEL(I)%LOC(J)%Y)
  END DO
 END DO

 END SUBROUTINE ST1SIMBOXWELLS

 !###======================================================================
 LOGICAL FUNCTION ST1CREATEIPF(NLAY,DIR)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NLAY
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER,ALLOCATABLE,DIMENSION(:) :: IU,LAYQ
 REAL,ALLOCATABLE,DIMENSION(:) :: TOP,BOT,L,KD,Q
 REAL,ALLOCATABLE,DIMENSION(:,:) :: TLP
 INTEGER :: I,J,K,IOS,ILAY,IROW,ICOL,JU,IUSUBRUN,IDT
 REAL :: Z1,Z2,ZCOR,MEANQ,SUMQ
 CHARACTER(LEN=256) :: LINE
 CHARACTER(LEN=50) :: QID
 LOGICAL :: LEX

 ST1CREATEIPF=.FALSE.

 CALL UTL_CREATEDIR(TRIM(DIR)//'\wells')

 ALLOCATE(IU(0:NLAY),TOP(STNLAY),BOT(STNLAY),L(STNLAY),KD(STNLAY),Q(STNLAY))
 ALLOCATE(LAYQ(NLAY))

 !## create subrun-file for modifications
 IUSUBRUN=UTL_GETUNIT()
 CALL OSD_OPEN(IUSUBRUN,FILE=TRIM(DIR)//'\scn_wel.run',STATUS='UNKNOWN',IOSTAT=IOS)
 IF(IOS.NE.0)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not create file: '//CHAR(13)// &
          TRIM(DIR)//'\scn_wel.run','Error')
  GOTO 10
 ENDIF

 IU=0

 !## open summary for well-strength
 IU(0)=UTL_GETUNIT()
 CALL OSD_OPEN(IU(0),FILE=TRIM(DIR)//'\wells\wells.ipf',STATUS='UNKNOWN',IOSTAT=IOS)
 IF(IOS.NE.0)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not create '//TRIM(DIR)//'\wells\wells.ipf','Error')
  GOTO 10
 ENDIF

 LINE=TRIM(ITOS(STNLAY))//',(wel)'
 WRITE(IUSUBRUN,'(A)') TRIM(LINE)
 DO ILAY=1,STNLAY
  IU(ILAY)=UTL_GETUNIT()
  CALL OSD_OPEN(IU(ILAY),FILE=TRIM(DIR)//'\wells\wells_l'//TRIM(ITOS(ILAY))//'.ipf',STATUS='UNKNOWN',IOSTAT=IOS)
  IF(IOS.NE.0)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not create '//TRIM(DIR)//'\wells\wells_l'// &
        TRIM(ITOS(ILAY))//'.ipf','Error')
   GOTO 10
  ENDIF
  LINE=TRIM(ITOS(ILAY))//',1.0,0.0,"'//TRIM(DIR)//'\wells\wells_l'//TRIM(ITOS(ILAY))//'.ipf"'
  WRITE(IUSUBRUN,'(A)') TRIM(LINE)
 ENDDO
 CLOSE(IUSUBRUN)

 DO ILAY=1,STNLAY
  IF(.NOT.IDFREAD(TOPIDF(ILAY),TOPIDF(ILAY)%FNAME,0))THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(TOPIDF(ILAY)%FNAME),'Error')
   GOTO 10
  ENDIF
  IF(.NOT.IDFREAD(BOTIDF(ILAY),BOTIDF(ILAY)%FNAME,0))THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(BOTIDF(ILAY)%FNAME),'Error')
   GOTO 10
  ENDIF
  IF(.NOT.IDFREAD(KDIDF(ILAY),KDIDF(ILAY)%FNAME,0))THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(KDIDF(ILAY)%FNAME),'Error')
   GOTO 10
  ENDIF
 END DO

 !## everything opened, write/read information for wells
 LINE=TRIM(ITOS(SUM(WEL(1:NWEL)%NLOC)))
 WRITE(IU(0),'(A)') TRIM(LINE)
 LINE=TRIM(ITOS(7+NLAY))
 WRITE(IU(0),'(A)') TRIM(LINE)
 WRITE(IU(0),'(A)') 'X'
 WRITE(IU(0),'(A)') 'Y'
 WRITE(IU(0),'(A)') 'Z1'
 WRITE(IU(0),'(A)') 'Z2'
 WRITE(IU(0),'(A)') 'Q-ID'
 WRITE(IU(0),'(A)') 'MEAN_Q (M3/DAG)'
 WRITE(IU(0),'(A)') 'SUM_Q (M3/SIMULATION)'
 DO ILAY=1,NLAY
  LINE='ILAY'//TRIM(ITOS(ILAY))
  WRITE(IU(0),'(A)') TRIM(LINE)
 END DO
 WRITE(IU(0),'(A)') '0,TXT'
 DO ILAY=1,NLAY
  LINE=TRIM(ITOS(SUM(WEL(1:NWEL)%NLOC)))
  WRITE(IU(ILAY),'(A)') TRIM(LINE)
  WRITE(IU(ILAY),'(A)') '6'
  WRITE(IU(ILAY),'(A)') 'X'
  WRITE(IU(ILAY),'(A)') 'Y'
  WRITE(IU(ILAY),'(A)') 'Q-ID'
  WRITE(IU(ILAY),'(A)') 'ILAY'
  WRITE(IU(ILAY),'(A)') 'Z1'
  WRITE(IU(ILAY),'(A)') 'Z2'
  WRITE(IU(ILAY),'(A)') '3,TXT'
 END DO

 !## number of well-systems
 DO I=1,NWEL

  !## allocate multiplication factors to be used to compute well-strengths
  ALLOCATE(TLP(WEL(I)%NLOC,STNLAY))
  TLP=0.0

  !## each well inside a wel-system ... get top/bottom,permeabilities
  DO J=1,WEL(I)%NLOC

   !## get top/bottoms
   DO ILAY=1,STNLAY
    TOP(ILAY)=TOPIDF(ILAY)%NODATA
    BOT(ILAY)=BOTIDF(ILAY)%NODATA
    CALL IDFIROWICOL(TOPIDF(ILAY),IROW,ICOL,WEL(I)%LOC(J)%X,WEL(I)%LOC(J)%Y)
    IF(IROW.NE.0.AND.ICOL.NE.0)TOP(ILAY)=IDFGETVAL(TOPIDF(ILAY),IROW,ICOL)
    CALL IDFIROWICOL(BOTIDF(ILAY),IROW,ICOL,WEL(I)%LOC(J)%X,WEL(I)%LOC(J)%Y)
    IF(IROW.NE.0.AND.ICOL.NE.0)BOT(ILAY)=IDFGETVAL(BOTIDF(ILAY),IROW,ICOL)
    !## kd averaging
    IF(WEL(I)%ITYPE.EQ.2)THEN
     KD(ILAY)=KDIDF(ILAY)%NODATA
     CALL IDFIROWICOL(KDIDF(ILAY),IROW,ICOL,WEL(I)%LOC(J)%X,WEL(I)%LOC(J)%Y)
     IF(IROW.NE.0.AND.ICOL.NE.0)KD(ILAY)=IDFGETVAL(KDIDF(ILAY),IROW,ICOL)
    ENDIF
   ENDDO

   !## surfacelevel -> msl
   ZCOR=0.0
   IF(WEL(I)%ILOCT.EQ.2)ZCOR=TOP(1)

   !## fit current position (compute length of well inside each modellayer)
   DO ILAY=1,STNLAY
    Z1=0.0
    Z2=0.0
    L(ILAY)=0.0
    IF(TOP(ILAY).NE.TOPIDF(ILAY)%NODATA.AND.BOT(ILAY).NE.BOTIDF(ILAY)%NODATA)THEN
     Z1=MIN(TOP(ILAY),ZCOR+WEL(I)%LOC(J)%Z1)
     Z2=MAX(BOT(ILAY),ZCOR+WEL(I)%LOC(J)%Z2)
     L(ILAY)=MAX(0.0,Z1-Z2)
    ENDIF
   END DO
   
!   CALL ST1CREATEIPF_GETTLP(STNLAY,IL,TL,L,TLP,KD,C)
   
   !## well within any aquifer(s)
   IF(SUM(L).GT.0.0)THEN
    !## compute percentage and include sumkd, only if itype.eq.2
    IF(WEL(I)%ITYPE.EQ.2)L=L*KD
    TLP(J,:)=(1.0/SUM(L))*L(:)  !## percentage (0-1) L*KD
   ENDIF

  ENDDO

  IF(SUM(TLP).GT.0.0)THEN
   !## overall distribution factors, inluding L*KD
   TLP(:,:)=(1.0/SUM(TLP))*TLP(:,:)
  ENDIF

  !## write results
  DO J=1,WEL(I)%NLOC

   SUMQ=0.0
   LAYQ=0

   !## multiply well strength with appropriate factor
   DO ILAY=1,STNLAY
    JU=UTL_GETUNIT()
    QID='well_sys'//TRIM(ITOS(I))//'_loc'//TRIM(ITOS(J))//'_l'//TRIM(ITOS(ILAY))
    CALL OSD_OPEN(JU,FILE=TRIM(DIR)//'\wells\'//TRIM(QID)//'.txt',STATUS='UNKNOWN',IOSTAT=IOS)
    IF(IOS.NE.0)THEN
     CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not create file: '//CHAR(13)// &
          TRIM(DIR)//'\wells\'//TRIM(QID)//'.txt','Error')
     GOTO 10
    ENDIF
    LINE=TRIM(ITOS(WEL(I)%NQ-1))  !## one less since last date is not connected to well strength
    WRITE(JU,'(A)') TRIM(LINE)
    WRITE(JU,'(A)') '2'
    WRITE(JU,'(A)') 'Date,-999'
    WRITE(JU,'(A)') 'Q(m3/d),-999'
    DO K=1,WEL(I)%NQ-1
     Q(ILAY)=TLP(J,ILAY)*WEL(I)%Q(K)%QRATE
     LINE=TRIM(ITOS(UTL_JDATETOIDATE(WEL(I)%Q(K)%IDATE)))//','//TRIM(RTOS(Q(ILAY),'F',2))
     WRITE(JU,'(A)') TRIM(LINE)
     IF(Q(ILAY).NE.0)THEN
      IDT=WEL(I)%Q(K+1)%IDATE-WEL(I)%Q(K)%IDATE
      SUMQ=SUMQ+Q(ILAY)*REAL(IDT) !## sumq for current system
      LAYQ(ILAY)=LAYQ(ILAY)+IDT    !## count for modellayer
     ENDIF
    ENDDO
    CLOSE(JU)

    LINE=TRIM(RTOS(WEL(I)%LOC(J)%X,'F',2))//','// &
         TRIM(RTOS(WEL(I)%LOC(J)%Y,'F',2))//','// &
         TRIM(QID)//','// &
         TRIM(ITOS(ILAY))//','// &
         TRIM(RTOS(WEL(I)%LOC(J)%Z1,'F',2))//','// &
         TRIM(RTOS(WEL(I)%LOC(J)%Z2,'F',2))
    WRITE(IU(ILAY),'(A)') TRIM(LINE)
   ENDDO

   MEANQ=SUMQ/REAL(SUM(LAYQ))
   LINE=TRIM(RTOS(WEL(I)%LOC(J)%X,'F',2))//','// &
        TRIM(RTOS(WEL(I)%LOC(J)%Y,'F',2))//','// &
        TRIM(RTOS(WEL(I)%LOC(J)%Z1,'F',2))//','// &
        TRIM(RTOS(WEL(I)%LOC(J)%Z2,'F',2))//','// &
        TRIM(QID)//','// &
        TRIM(RTOS(MEANQ,'F',2))//','// &
        TRIM(RTOS(SUMQ,'F',2))
   DO ILAY=1,NLAY; LINE=TRIM(LINE)//','//TRIM(ITOS(LAYQ(ILAY))); END DO
   WRITE(IU(0),'(A)') TRIM(LINE)

  ENDDO
  DEALLOCATE(TLP)

 ENDDO

 ST1CREATEIPF=.TRUE.

10 CONTINUE

 DO ILAY=1,STNLAY
  IF(TOPIDF(ILAY)%IU.GT.0)THEN
   INQUIRE(FILE=TOPIDF(ILAY)%FNAME,OPENED=LEX)
   IF(LEX)CLOSE(TOPIDF(ILAY)%IU)
   TOPIDF(ILAY)%IU=0
  ENDIF
  IF(BOTIDF(ILAY)%IU.GT.0)THEN
   INQUIRE(FILE=BOTIDF(ILAY)%FNAME,OPENED=LEX)
   IF(LEX)CLOSE(BOTIDF(ILAY)%IU)
   BOTIDF(ILAY)%IU=0
  ENDIF
  IF(KDIDF(ILAY)%IU.GT.0)THEN
   INQUIRE(FILE=KDIDF(ILAY)%FNAME,OPENED=LEX)
   IF(LEX)CLOSE(KDIDF(ILAY)%IU)
   KDIDF(ILAY)%IU=0
  ENDIF
 ENDDO

 DO ILAY=0,NLAY
  IF(IU(ILAY).GT.0)THEN
   INQUIRE(IU(ILAY),OPENED=LEX)
   IF(LEX)CLOSE(IU(ILAY))
  ENDIF
 ENDDO

 DEALLOCATE(IU,TOP,BOT,L,KD,Q,LAYQ)

 END FUNCTION ST1CREATEIPF

 !###======================================================================
 LOGICAL FUNCTION ST1CREATEIPF_STEADY(IKD,ICW,MAXC,MINKH,IMIDF,SDATE,EDATE,ISS,ICLAY,IFRAC,HNODATA,FNODATA,MINKD)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IKD,ICW,ISS,IMIDF,ICLAY,IFRAC
 INTEGER,INTENT(INOUT) :: SDATE,EDATE
 REAL,INTENT(IN) :: MAXC,MINKH,HNODATA,FNODATA,MINKD
 INTEGER,ALLOCATABLE,DIMENSION(:) :: IU
 REAL,ALLOCATABLE,DIMENSION(:) :: TOP,BOT,KD,Q,C,KH,TLP
 INTEGER :: I,J,K,IOS,ILAY,JLAY,IROW,ICOL,JU,IUSUBRUN,IDT
 REAL :: Z1,Z2,ZCOR,MEANQ,FMID
 CHARACTER(LEN=521) :: LINE
 CHARACTER(LEN=256) :: DIR,FNAME
 CHARACTER(LEN=50) :: QID
 LOGICAL :: LEX
 !REAL,PARAMETER :: FNODATA=-9999.0 !## nodata for filterstellingen (nog in imodbatch)
   
 ST1CREATEIPF_STEADY=.FALSE.

 SDATE=UTL_IDATETOJDATE(SDATE); EDATE=UTL_IDATETOJDATE(EDATE)
 
 DIR=IPF(1)%FNAME(:INDEX(IPF(1)%FNAME,'.',.TRUE.)-1)
 CALL UTL_CREATEDIR(TRIM(DIR))
 
 WRITE(*,'(/A/)') 'Results saved in folder: ['//TRIM(DIR)//']'

 ALLOCATE(IU(-1:STNLAY),TOP(STNLAY),BOT(STNLAY),KD(STNLAY),KH(STNLAY),Q(STNLAY),C(STNLAY-1))
 
 IU=0

 DO ILAY=-1,STNLAY
  IU(ILAY)=UTL_GETUNIT()
  CALL OSD_OPEN(IU(ILAY),FILE=TRIM(DIR)//'\imod_mkipf_wells_l'//TRIM(ITOS(ILAY))//'_tmp.ipf',STATUS='UNKNOWN',IOSTAT=IOS)
  IF(IOS.NE.0)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not create '//TRIM(DIR)//'\imod_mkipf_wells_l'// &
        TRIM(ITOS(ILAY))//'_tmp.ipf','Error'); GOTO 10
  ENDIF
 ENDDO

 DO ILAY=1,STNLAY
  IF(.NOT.IDFREAD(TOPIDF(ILAY),TOPIDF(ILAY)%FNAME,0))THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(TOPIDF(ILAY)%FNAME),'Error'); GOTO 10
  ENDIF
  IF(.NOT.IDFREAD(BOTIDF(ILAY),BOTIDF(ILAY)%FNAME,0))THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(BOTIDF(ILAY)%FNAME),'Error'); GOTO 10
  ENDIF
  IF(IKD.EQ.1)THEN
   IF(.NOT.IDFREAD(KDIDF(ILAY),KDIDF(ILAY)%FNAME,0))THEN
    CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(KDIDF(ILAY)%FNAME),'Error'); GOTO 10
   ENDIF
  ENDIF
  IF(ILAY.LT.STNLAY)THEN
   IF(ICW.EQ.1)THEN
    IF(.NOT.IDFREAD(CIDF(ILAY),CIDF(ILAY)%FNAME,0))THEN
     CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Can not read '//TRIM(CIDF(ILAY)%FNAME),'Error'); GOTO 10
    ENDIF
   ENDIF
  ENDIF
 END DO

 !## everything opened, write/read information for wells
 DO ILAY=-1,STNLAY
  LINE=TRIM(ITOS(IPF(1)%NROW));   WRITE(IU(ILAY),'(A)') TRIM(LINE)
  LINE=TRIM(ITOS(IPF(1)%NCOL+7)); WRITE(IU(ILAY),'(A)') TRIM(LINE)
  DO I=1,IPF(1)%NCOL
   IF(ISS.EQ.0)THEN
    IF(I.EQ.IPF(1)%QCOL)WRITE(IU(ILAY),'(A)') 'Q_ASSIGNED'
    IF(I.NE.IPF(1)%QCOL)WRITE(IU(ILAY),'(A)') TRIM(IPF(1)%ATTRIB(I))
   ELSEIF(ISS.EQ.1)THEN
    IF(I.EQ.IPF(1)%ACOL)WRITE(IU(ILAY),'(A)') 'Q_ASSIGNED'
    IF(I.NE.IPF(1)%ACOL)WRITE(IU(ILAY),'(A)') TRIM(IPF(1)%ATTRIB(I))
   ENDIF
  ENDDO
  WRITE(IU(ILAY),'(A)') 'Q_ORG' 
  WRITE(IU(ILAY),'(A)') 'FRACTION'
  WRITE(IU(ILAY),'(A)') 'ILAY'
  WRITE(IU(ILAY),'(A)') 'TOP'
  WRITE(IU(ILAY),'(A)') 'BOT'
  WRITE(IU(ILAY),'(A)') 'KD'
  WRITE(IU(ILAY),'(A)') 'ERROR_CODE'
  IF(ISS.EQ.0)THEN
   LINE=TRIM(ITOS(IPF(1)%ACOL))//','//TRIM(IPF(1)%FEXT) 
  ELSE
   LINE=TRIM(ITOS(IPF(1)%NCOL+1))//','//TRIM(IPF(1)%FEXT) 
  ENDIF
  WRITE(IU(ILAY),'(A)') TRIM(LINE)
 END DO

 !## allocate multiplication factors to be used to compute well-strengths
 ALLOCATE(TLP(STNLAY)); TLP=0.0

 !## each well inside a wel-system ... get top/bottom,permeabilities
 DO I=1,IPF(1)%NROW

  !## get top/bottoms
  DO ILAY=1,STNLAY

   !## read current top value
   TOP(ILAY)=TOPIDF(ILAY)%NODATA
   CALL IDFIROWICOL(TOPIDF(ILAY),IROW,ICOL,IPF(1)%XYZ(1,I),IPF(1)%XYZ(2,I))
   IF(IROW.NE.0.AND.ICOL.NE.0)TOP(ILAY)=IDFGETVAL(TOPIDF(ILAY),IROW,ICOL)

   !## read current bot value
   BOT(ILAY)=BOTIDF(ILAY)%NODATA
   CALL IDFIROWICOL(BOTIDF(ILAY),IROW,ICOL,IPF(1)%XYZ(1,I),IPF(1)%XYZ(2,I))
   IF(IROW.NE.0.AND.ICOL.NE.0)BOT(ILAY)=IDFGETVAL(BOTIDF(ILAY),IROW,ICOL)

   !## read current kd value
   KD(ILAY)=1.0
   IF(IKD.EQ.1)THEN
    KD(ILAY)=KDIDF(ILAY)%NODATA
    CALL IDFIROWICOL(KDIDF(ILAY),IROW,ICOL,IPF(1)%XYZ(1,I),IPF(1)%XYZ(2,I))
    IF(IROW.NE.0.AND.ICOL.NE.0)KD(ILAY)=IDFGETVAL(KDIDF(ILAY),IROW,ICOL)
   ENDIF

   !## read current c value
   IF(ILAY.LT.STNLAY)THEN
    C(ILAY)=1.0
    IF(ICW.EQ.1)THEN
     C(ILAY)=CIDF(ILAY)%NODATA
     CALL IDFIROWICOL(CIDF(ILAY),IROW,ICOL,IPF(1)%XYZ(1,I),IPF(1)%XYZ(2,I))
     IF(IROW.NE.0.AND.ICOL.NE.0)C(ILAY)=IDFGETVAL(CIDF(ILAY),IROW,ICOL)
    ENDIF
   ENDIF

  ENDDO

  !## fit current position (compute length of well inside each modellayer)
  DO ILAY=1,STNLAY
   Z1=0.0; Z2=0.0; KH(ILAY)=0.0
   IF(TOP(ILAY).NE.TOPIDF(ILAY)%NODATA.AND.BOT(ILAY).NE.BOTIDF(ILAY)%NODATA)THEN
    IF(TOP(ILAY)-BOT(ILAY).GT.0.0)KH(ILAY)=KD(ILAY)/(TOP(ILAY)-BOT(ILAY))
   ENDIF
  END DO
 
  FMID=FNODATA
  IF(IMIDF.EQ.1)THEN
   IF(IPF(1)%XYZ(3,I).NE.FNODATA.AND.IPF(1)%XYZ(4,I).NE.FNODATA)THEN
    FMID=(IPF(1)%XYZ(3,I)+IPF(1)%XYZ(4,I))/2.0
   ELSEIF(IPF(1)%XYZ(3,I).NE.FNODATA)THEN
    FMID=IPF(1)%XYZ(3,I)
   ELSEIF(IPF(1)%XYZ(4,I).NE.FNODATA)THEN
    FMID=IPF(1)%XYZ(4,I)
   ENDIF
   IPF(1)%XYZ(3,I)=FMID; IPF(1)%XYZ(4,I)=FMID   
  ELSE
   IF(IPF(1)%XYZ(3,I).NE.FNODATA.AND.IPF(1)%XYZ(4,I).NE.FNODATA)THEN
    FMID=(IPF(1)%XYZ(3,I)+IPF(1)%XYZ(4,I))/2.0
   ENDIF
  ENDIF
  TLP=0.0
  IF(FMID.NE.FNODATA)THEN
   IF(IPF(1)%ACOL.GT.0)THEN
    CALL PMANAGER_SAVEMF2005_PCK_GETTLP(STNLAY,TLP,KH,TOP,BOT,IPF(1)%XYZ(3,I),IPF(1)%XYZ(4,I))
   ELSE
    CALL PMANAGER_SAVEMF2005_PCK_GETTLP(STNLAY,TLP,KH,TOP,BOT,IPF(1)%XYZ(3,I),IPF(1)%XYZ(4,I)) 
   ENDIF   
   IF(ISS.GT.0.AND.IPF(1)%ACOL.GT.0)THEN
    FNAME=IPF(1)%FNAME(:INDEX(IPF(1)%FNAME,'\',.TRUE.)-1)//'\'// &
          TRIM(IPF(1)%INFO(IPF(1)%ACOL,I))//'.txt'
    IF(.NOT.PMANAGER_SAVEMF2005_PCK_READTXT(2,SDATE,EDATE,2,IPF(1)%XYZ(5,I),FNAME))THEN
     WRITE(*,'(A)') 'Error IPF-TXT '//TRIM(IPF(1)%FNAME)//' failed'
    ENDIF
   ENDIF
  ELSE
   IF(MINKD.GT.0.0)THEN
    !## overrule, neem eerste wvp met kD groter dan minkh ...
    DO ILAY=1,STNLAY
     IF((TOP(ILAY)-BOT(ILAY))*KH(ILAY).GT.MINKD)EXIT
    ENDDO
    TLP=0.0
    IF(ILAY.LE.STNLAY)THEN
     TLP(ILAY)=1.0
    ELSE
     DO ILAY=1,STNLAY; WRITE(*,*) ILAY,TOP(ILAY),BOT(ILAY),KH(ILAY),(TOP(ILAY)-BOT(ILAY))*KH(ILAY); ENDDO
     IPF(1)%XYZ(5,I)=HNODATA
    ENDIF
   ELSE
    IPF(1)%XYZ(5,I)=HNODATA
   ENDIF
  ENDIF
  
  !## multiply well strength with appropriate factor
  DO ILAY=1,STNLAY
   IF(IPF(1)%XYZ(5,I).EQ.HNODATA)THEN
    Q(ILAY)=HNODATA
   ELSE
    IF(TLP(ILAY).EQ.0.0)THEN
     Q(ILAY)=HNODATA
    ELSE
     Q(ILAY)=ABS(TLP(ILAY))*IPF(1)%XYZ(5,I)
    ENDIF
   ENDIF

   LINE='"'//TRIM(IPF(1)%INFO(1,I))//'"'
   DO J=2,IPF(1)%NCOL
    IF(ISS.EQ.0)THEN
     IF(J.EQ.IPF(1)%QCOL)LINE=TRIM(LINE)//',"'//TRIM(RTOS(Q(ILAY),'F',3))//'"'
     IF(J.NE.IPF(1)%QCOL)LINE=TRIM(LINE)//',"'//TRIM(IPF(1)%INFO(J,I))//'"'
    ELSEIF(ISS.EQ.1)THEN
     IF(J.EQ.IPF(1)%ACOL)LINE=TRIM(LINE)//',"'//TRIM(RTOS(Q(ILAY),'F',3))//'"'
     IF(J.NE.IPF(1)%ACOL)LINE=TRIM(LINE)//',"'//TRIM(IPF(1)%INFO(J,I))//'"'
    ENDIF
   ENDDO
   IF(ISS.EQ.0)THEN
    LINE=TRIM(LINE)//',"'//TRIM(IPF(1)%INFO(IPF(1)%QCOL,I))//'","'//TRIM(RTOS(ABS(TLP(ILAY)),'F',4))//'"'
   ELSEIF(ISS.EQ.1)THEN
    LINE=TRIM(LINE)//',"'//TRIM(IPF(1)%INFO(IPF(1)%ACOL,I))//'","'//TRIM(RTOS(ABS(TLP(ILAY)),'F',4))//'"'
   ENDIF
   LINE=TRIM(LINE)//',"'//TRIM(ITOS(ILAY))//'","'//TRIM(RTOS(TOP(ILAY),'F',2))//'","'// &
                          TRIM(RTOS(BOT(ILAY),'F',2))//'","'//TRIM(RTOS(KD(ILAY),'F',2))//'"'
   IF(SUM(TLP).EQ.0.0)THEN
    LINE=TRIM(LINE)//',"#"'
   ELSEIF(SUM(TLP).LT.0.0)THEN
    LINE=TRIM(LINE)//',"@"'
   ELSE
    LINE=TRIM(LINE)//',"-"'
   ENDIF
   WRITE(IU(ILAY),'(A)') TRIM(LINE)
   IF(SUM(TLP).LE.0.0)WRITE(IU(-1),'(A)') TRIM(LINE)
   WRITE(IU(0),'(A)') TRIM(LINE)
  ENDDO

 ENDDO
 DO ILAY=-1,STNLAY; CLOSE(IU(ILAY)); ENDDO

 DEALLOCATE(TLP)

 IF(ISS.EQ.1)IPF(1)%QCOL=IPF(1)%ACOL
 DO ILAY=1,STNLAY
  IPF(1)%FNAME=TRIM(DIR)//'\imod_mkipf_wells_l'//TRIM(ITOS(ILAY))//'_tmp.ipf'
  !# read entire ipf
  IF(.NOT.IPFREAD2(1,1,0))RETURN
  J=0; DO I=1,IPF(1)%NROW
   !## extraction ne zero
   IF(IPF(1)%XYZ(5,I).NE.HNODATA)THEN !## always on five
    J=J+1
    IF(J.NE.I)THEN; IPF(1)%INFO(:,J)=IPF(1)%INFO(:,I); ENDIF
   ENDIF
  ENDDO
  IPF(1)%NROW=J
  IF(IPF(1)%NROW.GT.0)THEN
   IPF(1)%FNAME=TRIM(DIR)//'\imod_mkipf_wells_l'//TRIM(ITOS(ILAY))//'.ipf'
   IF(.NOT.IPFWRITE(1))THEN
   ENDIF
  ENDIF
  CALL IOSDELETEFILE(TRIM(DIR)//'\imod_mkipf_wells_l'//TRIM(ITOS(ILAY))//'_tmp.ipf')
 ENDDO

 ST1CREATEIPF_STEADY=.TRUE.

10 CONTINUE

 CALL IDFDEALLOCATE(TOPIDF,SIZE(TOPIDF)); CALL IDFDEALLOCATE(BOTIDF,SIZE(BOTIDF))
 CALL IDFDEALLOCATE(KDIDF,SIZE(KDIDF));   CALL IDFDEALLOCATE(CIDF,SIZE(CIDF))
 DO ILAY=0,STNLAY; IF(IU(ILAY).GT.0)THEN; INQUIRE(IU(ILAY),OPENED=LEX); IF(LEX)CLOSE(IU(ILAY)); ENDIF; ENDDO

 DEALLOCATE(IU,TOP,BOT,KD,Q,C,KH)

 END FUNCTION ST1CREATEIPF_STEADY

 !###======================================================================
 SUBROUTINE STWEL1DELETE()
 !###======================================================================
 IMPLICIT NONE

 CALL WDIALOGSELECT(ID_DSCENTOOLTAB1)
 CALL WDIALOGGETMENU(IDF_MENU1,IWEL)

 CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'Are you sure to delete the definition for:'//CHAR(13)// &
  TRIM(WEL(IWEL)%CNAME)//' ?','Question')
 IF(WINFODIALOG(4).NE.1)RETURN

 END SUBROUTINE STWEL1DELETE

 !###======================================================================
 SUBROUTINE STWEL1PUTFIELDS()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J,IY,IM,ID

 !## clear fields
 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB1)
 DO I=1,NROWQ
  DO J=1,4; CALL WGRIDCLEARCELL(IDF_GRID1,J,I); END DO
 END DO
 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB2)
 DO I=1,NROWL
  DO J=1,5; CALL WGRIDCLEARCELL(IDF_GRID1,J,I); END DO
 ENDDO

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB1)
 CALL WDIALOGPUTSTRING(IDF_STRING1,WEL(IWEL)%CNAME)

 IF(WEL(IWEL)%NQ.GT.0)THEN
  CALL WDIALOGPUTSTRING(IDF_STRING2,JDATETOGDATE(WEL(IWEL)%Q(1)%IDATE))
 ELSE
  CALL IOSDATE(IY,IM,ID)
  CALL WDIALOGPUTSTRING(IDF_STRING2,TRIM(ITOS(ID))//'/'//TRIM(ITOS(IM))//'/'//TRIM(ITOS(IY)))
 ENDIF

 DO I=1,WEL(IWEL)%NQ
  CALL WGRIDPUTCELLSTRING(IDF_GRID1,1,I,JDATETOGDATE(WEL(IWEL)%Q(I)%IDATE))
  IF(I.LT.WEL(IWEL)%NQ)CALL WGRIDPUTCELLREAL(IDF_GRID1,3,I,WEL(IWEL)%Q(I)%QRATE/24.0)  !## m3/day -> m3/hr
  CALL WGRIDPUTCELLINTEGER(IDF_GRID1,4,I,WEL(IWEL)%Q(I)%IDATE)
 END DO
 CALL WGRIDSETCELL(IDF_GRID1,1,1)!,IPOS)

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB2)

 DO I=1,WEL(IWEL)%NLOC
  CALL WGRIDPUTCELLSTRING(IDF_GRID1,1,I,WEL(IWEL)%LOC(I)%ID)
  CALL WGRIDPUTCELLREAL(IDF_GRID1,2,I,WEL(IWEL)%LOC(I)%Z1)
  CALL WGRIDPUTCELLREAL(IDF_GRID1,3,I,WEL(IWEL)%LOC(I)%Z2)
  CALL WGRIDPUTCELLREAL(IDF_GRID1,4,I,WEL(IWEL)%LOC(I)%X)
  CALL WGRIDPUTCELLREAL(IDF_GRID1,5,I,WEL(IWEL)%LOC(I)%Y)
 END DO
 CALL WGRIDSETCELL(IDF_GRID1,1,1)!,IPOS)

 IF(WEL(IWEL)%ILOCT.EQ.1)CALL WDIALOGPUTRADIOBUTTON(IDF_RADIO1)
 IF(WEL(IWEL)%ILOCT.EQ.2)CALL WDIALOGPUTRADIOBUTTON(IDF_RADIO2)
 IF(WEL(IWEL)%ITYPE.EQ.1)CALL WDIALOGPUTRADIOBUTTON(IDF_RADIO3)
 IF(WEL(IWEL)%ITYPE.EQ.2)CALL WDIALOGPUTRADIOBUTTON(IDF_RADIO4)

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB3)
 CALL WDIALOGPUTOPTION(IDF_MENU1,WEL(IWEL)%ISYMBOL)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER1,WEL(IWEL)%ICLR)

 END SUBROUTINE STWEL1PUTFIELDS

 !###======================================================================
 LOGICAL FUNCTION STWEL1GETFIELDS()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I

 STWEL1GETFIELDS=.FALSE.

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB1)
 CALL WDIALOGGETSTRING(IDF_STRING1,WEL(IWEL)%CNAME)
 IF(LEN_TRIM(WEL(IWEL)%CNAME).EQ.0)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You should specify a name for the Wel-system at least','Warning')
  RETURN
 ENDIF

 CALL WDIALOGUNDEFINED(RVALUE=NODATAGRID)

 WEL(IWEL)%NQ=0
 DO I=1,NROWQ
  CALL WGRIDGETCELLSTRING(IDF_GRID1,1,I,CDATE)
  !## empty
  IF(LEN_TRIM(CDATE).EQ.0)EXIT
  CALL WGRIDGETCELLREAL(IDF_GRID1,3,I,WEL(IWEL)%Q(I)%QRATE)
  WEL(IWEL)%Q(I)%IDATE=GDATETOJDATE(CDATE)
  !## date conversion went wrong
  IF(WEL(IWEL)%Q(I)%IDATE.EQ.0)EXIT
  WEL(IWEL)%NQ=WEL(IWEL)%NQ+1
  WEL(IWEL)%Q(I)%QRATE=WEL(IWEL)%Q(I)%QRATE*24.0  !## m3/hr -> m3/day
 END DO

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB2)

 WEL(IWEL)%NLOC=0

 DO I=1,NROWL
  CALL WGRIDGETCELLSTRING(IDF_GRID1,1,I,WEL(IWEL)%LOC(I)%ID)
  IF(LEN_TRIM(WEL(IWEL)%LOC(I)%ID).EQ.0)EXIT
  CALL WGRIDGETCELLREAL(IDF_GRID1,2,I,WEL(IWEL)%LOC(I)%Z1)
  IF(WEL(IWEL)%LOC(I)%Z1.EQ.NODATAGRID)EXIT
  CALL WGRIDGETCELLREAL(IDF_GRID1,3,I,WEL(IWEL)%LOC(I)%Z2)
  IF(WEL(IWEL)%LOC(I)%Z2.EQ.NODATAGRID)EXIT

  IF(WEL(IWEL)%LOC(I)%Z1.LT.WEL(IWEL)%LOC(I)%Z2)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You should specify top filter to be above the bottom filter'//CHAR(13)// &
    'Check well number '//TRIM(ITOS(I)),'Warning')
   RETURN
  ENDIF

  CALL WGRIDGETCELLREAL(IDF_GRID1,4,I,WEL(IWEL)%LOC(I)%X)
  IF(WEL(IWEL)%LOC(I)%X.EQ.NODATAGRID)EXIT
  CALL WGRIDGETCELLREAL(IDF_GRID1,5,I,WEL(IWEL)%LOC(I)%Y)
  IF(WEL(IWEL)%LOC(I)%Y.EQ.NODATAGRID)EXIT
  WEL(IWEL)%NLOC=WEL(IWEL)%NLOC+1
 END DO

 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO1,WEL(IWEL)%ILOCT)
 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO3,WEL(IWEL)%ITYPE)

 CALL WDIALOGSELECT(ID_DSCENTOOL_PROPTAB3)
 CALL WDIALOGGETMENU(IDF_MENU1,WEL(IWEL)%ISYMBOL)
 CALL WDIALOGGETINTEGER(IDF_INTEGER1,WEL(IWEL)%ICLR)

 STWEL1GETFIELDS=.TRUE.

 END FUNCTION STWEL1GETFIELDS

END MODULE MOD_SCENTOOL_WELLS