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

USE WINTERACTER
USE RESOURCE
USE OPENGL
USE MOD_COLOURS
USE IMODVAR, ONLY : DP_KIND,SP_KIND,PI,REPLACESTRING
USE MOD_IDF_PAR
USE MODPLOT, ONLY : MPW,MP
USE MOD_QKSORT, ONLY : UTL_WSORT
USE MOD_ASC2IDF_PAR, ONLY : IDFFILE,IGRIDFUNC,CS,NODATA,XYZFNAMES,IXCOL,IYCOL,IZCOL,TRIMDEPTH_IDF,ELLIPS_IDF
USE MOD_PREF_PAR, ONLY : PREFVAL
USE MOD_UTL, ONLY : ITOS,UTL_GETUNIT,RTOS,JD,UTL_IDATETOJDATE,UTL_WAITMESSAGE,UTL_IDFGETDATE,UTL_CAP, &
           UTL_GETUNIQUE_CHAR,UTL_CREATEDIR,UTL_IDFSNAPTOGRID,UTL_WSELECTFILE,UTL_SUBST,UTL_DIALOGSHOW
USE MOD_IDF, ONLY : IDFREAD,IDFDEALLOCATEX,IDFALLOCATEX,IDFWRITE,IDFIROWICOL,IDFREADSCALE,IDFGETLOC,IDFCOPY,IDFDEALLOCATESX
USE MOD_PLINES_PAR
USE MOD_PLINES_SP, ONLY : TRACEPREPARESP,TRACEREADSP
USE MOD_PLINES_FLOLIN, ONLY : FLOLIN
USE MOD_OSD, ONLY : OSD_OPEN,OSD_IOSTAT_MSG
USE MOD_MANAGER_UTL, ONLY : MANAGER_UTL_DELETE,MANAGER_UTL_ADDFILE
USE MOD_DEMO_PAR
USE MOD_3D_PAR, ONLY : PLLISTINDEX,PLLISTCLR,PLLISTAGE,STPLISTINDEX,MIDPOS,TOP,BOT,IDFPLOT,IPATHLINE_3D
USE MOD_IPF, ONLY : IPFALLOCATE,IPFDEALLOCATE,IPFREAD2,NIPF,IPFFILE=>IPF

REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE,PRIVATE :: XSP,YSP,ZSP
INTEGER,DIMENSION(:),ALLOCATABLE,PRIVATE :: ILAYERS,JLAYERS,NSGRP
CHARACTER(LEN=15),DIMENSION(:),ALLOCATABLE,PRIVATE :: CLAYERS

CONTAINS

 !###======================================================================
 LOGICAL FUNCTION TRACE_3D_INIT(FNAME)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: FNAME
 CHARACTER(LEN=256) :: RUNFILE
 INTEGER :: I,IBATCH,ILAY,IPLOT,NG,IDEF,NCONS

 TRACE_3D_INIT=.FALSE.

 IDEF=WMENUGETSTATE(ID_INTERACTIVEPATHLINES_DEF,2)

 IF(TRIM(FNAME).EQ.'')THEN
  RUNFILE=''
  IF(.NOT.UTL_WSELECTFILE('iMODPATH runfile (*.run)|*.run|', &
                   LOADDIALOG+PROMPTON+DIRCHANGE+APPENDEXT, &
                   RUNFILE,'Select iMODPATH runfile (*.run)'))RETURN
 ELSE
  RUNFILE=FNAME
 ENDIF 
  
 CALL WDIALOGLOAD(ID_D3DSETTINGS_LAYER,ID_D3DSETTINGS_LAYER)
 CALL WDIALOGLOAD(ID_D3DSETTINGS_SINKS,ID_D3DSETTINGS_SINKS)
 CALL WDIALOGLOAD(ID_D3DSETTINGS_LOCATION,ID_D3DSETTINGS_LOCATION)
 CALL WDIALOGLOAD(ID_D3DSETTINGS_IDF,ID_D3DSETTINGS_IDF)
 CALL WDIALOGLOAD(ID_D3DSETTINGS_PART,ID_D3DSETTINGS_PART)
 
 !## deallocate all memory
 CALL TRACEDEALLOCATE(1)

 !## read runfile
 IBATCH=0; IF(.NOT.TRACEREADRUNFILE(RUNFILE,IBATCH))RETURN

 !## take head for determination model-dimensions
 IF(.NOT.IDFREAD(IDF,HFFNAME(1,1,1),0))RETURN
 !## define zoom area for particle tracking - overrule area
 IDF%XMIN=MAX(IDF%XMIN,MPW%XMIN); IDF%XMAX=MIN(IDF%XMAX,MPW%XMAX)
 IDF%YMIN=MAX(IDF%YMIN,MPW%YMIN); IDF%YMAX=MIN(IDF%YMAX,MPW%YMAX)
 CALL UTL_IDFSNAPTOGRID(IDF%XMIN,IDF%XMAX,IDF%YMIN,IDF%YMAX,IDF%DX,IDF%NCOL,IDF%NROW)

 IF(IDF%NCOL.LE.0.OR.IDF%NROW.LE.0)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'The selected results of the iMODPATH runfile are not within'//CHAR(13)// &
      'the limits of the current graphical zoom extent','Error')
  CALL TRACE_3D_CLOSE()
  RETURN
 ENDIF

 !## read data 
 IF(.NOT.TRACECALC_INIT(0,NCONS))RETURN
 !## read fluxes, default forward simulation
 IF(.NOT.TRACEREADBUDGET(1,0))RETURN 

 !## initiate tracing variables
 PL%IREV =0     !## 0=forward; 1=backward simulation (operates as a switch)
 PL%NPART=0     !## number of particles active
 PL%IPLOTSP=0   !## plotsp
 PL%SPCOLOR=WRGB(200,10,10) !## initial colour of start-locations
 NSPG=0
 
 ALLOCATE(IVISIT(IDF%NCOL*IDF%NROW*NLAY)); IVISIT=INT(0,1)
 ALLOCATE(LVISIT(IDF%NCOL*IDF%NROW*NLAY)); LVISIT=0

 IF(IDEF.EQ.1)THEN
 
  !## fill top/bot to imod manager
  !## select files in the imod manager
  MP%ISEL=.FALSE.
  DO ILAY=1,NLAY
   DO IPLOT=1,SIZE(MP) 
    IF(TRIM(UTL_CAP(MP(IPLOT)%IDFNAME,'U')).EQ.TRIM(UTL_CAP(ITBFNAME(2,ILAY),'U')))MP(IPLOT)%ISEL=.TRUE.
    IF(TRIM(UTL_CAP(MP(IPLOT)%IDFNAME,'U')).EQ.TRIM(UTL_CAP(ITBFNAME(3,ILAY),'U')))MP(IPLOT)%ISEL=.TRUE.
   ENDDO       
  END DO

  !## delete them all from manager
  CALL MANAGER_UTL_DELETE(IQ=0)

  DO ILAY=1,NLAY
   !## top
   IF(INDEX(UTL_CAP(ITBFNAME(2,ILAY),'U'),'.IDF').GT.0)CALL MANAGER_UTL_ADDFILE(IDFNAMEGIVEN=ITBFNAME(2,ILAY),LDEACTIVATE=.FALSE.)
   !## bot
   IF(INDEX(UTL_CAP(ITBFNAME(3,ILAY),'U'),'.IDF').GT.0)CALL MANAGER_UTL_ADDFILE(IDFNAMEGIVEN=ITBFNAME(3,ILAY),LDEACTIVATE=.FALSE.)
  ENDDO
 
 ENDIF
 
 ALLOCATE(ILAYERS(NLAY),CLAYERS(NLAY))
 DO ILAY=1,NLAY; CLAYERS(ILAY)='Layer '//TRIM(ITOS(ILAY)); ENDDO; ILAYERS=1
 CALL WDIALOGSELECT(ID_D3DSETTINGS_LAYER)
 CALL WDIALOGPUTMENU(IDF_MENU1,CLAYERS,NLAY,ILAYERS)
 CALL WDIALOGSELECT(ID_D3DSETTINGS_SINKS)
 CALL WDIALOGPUTMENU(IDF_MENU1,CLAYERS,NLAY,ILAYERS)
 DEALLOCATE(ILAYERS,CLAYERS)
 CALL WDIALOGSELECT(ID_D3DSETTINGS_LOCATION)
 CALL WDIALOGSELECT(ID_D3DSETTINGS_IDF)

 CALL WDIALOGSELECT(ID_D3DSETTINGS_PART)
 CALL WDIALOGPUTIMAGE(ID_OPEN,ID_ICONOPEN,1)
 CALL WDIALOGPUTIMAGE(ID_SAVEAS,ID_ICONSAVEAS,1)
 CALL WDIALOGPUTIMAGE(ID_DELETE,ID_ICONDELETE,1)
 NG=WINFOGRID(IDF_GRID1,GRIDROWSMAX)

 ALLOCATE(STPLISTINDEX(NG)); STPLISTINDEX=0

 !## initiate particles
 CALL TRACE_INIT_SP(NG)
 !## maximal number of particles to be traced
 DO I=1,NG; CALL TRACE_AL_SP(1000,I); ENDDO

 IF(IDEF.EQ.1)THEN
  DEMO%IDEMO=2; DEMO%CONFLAG=4; DEMO%ACCFLAG=3; DEMO%IFILL=1
 ELSE
  DEMO%IDEMO=0
 ENDIF
 
 CALL WMENUSETSTATE(ID_INTERACTIVEPATHLINES,2,1)
 CALL WMENUSETSTATE(ID_INTERACTIVEPATHLINES,1,0)
 
 TRACE_3D_INIT=.TRUE.
 
 END FUNCTION TRACE_3D_INIT
 
 !###======================================================================
 LOGICAL FUNCTION TRACE_3D_RESET(ICODE) 
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ICODE
 INTEGER :: I,J,N,M
 
 TRACE_3D_RESET=.FALSE.
 
 IF(ALLOCATED(PLLISTINDEX))THEN
  DO I=1,SIZE(PLLISTINDEX,1)
   DO J=1,SIZE(PLLISTINDEX,2)
    IF(PLLISTINDEX(I,J).NE.0)THEN
     CALL GLDELETELISTS(PLLISTINDEX(I,J),1_GLSIZEI)
     CALL WINDOWOUTSTATUSBAR(4,'Clearing memory for Pathines '//TRIM(ITOS(I))//'...')
    ENDIF
   ENDDO
  ENDDO
  DEALLOCATE(PLLISTINDEX)
 ENDIF
 IF(ALLOCATED(PLLISTAGE))  DEALLOCATE(PLLISTAGE)
 IF(ALLOCATED(PLLISTCLR))  DEALLOCATE(PLLISTCLR)
  
 CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 !## length of the particles tail
 CALL WDIALOGGETINTEGER(IDF_INTEGER4,N)
 !## number of groups active
 CALL WDIALOGGETINTEGER(IDF_INTEGER11,M)
 !## no group active
 IF(M.LE.0)THEN
  IF(ICODE.EQ.1)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Create a set of particles first!','Information')
  RETURN
 ENDIF
 
 !## drawinglist index
 ALLOCATE(PLLISTINDEX(N,M)); PLLISTINDEX=0
 !## colour fraction
 ALLOCATE(PLLISTCLR(N)); PLLISTCLR=0
 !## age
 ALLOCATE(PLLISTAGE(N)); PLLISTAGE=0.0D0
 
 PL%TCUR=0.0D0
 PL%NPER=0
 PL%NTIME=0
 
 CALL WDIALOGPUTDOUBLE(IDF_REAL7,PL%TCUR) 
 
 TRACE_3D_RESET=.TRUE.
     
 END FUNCTION TRACE_3D_RESET

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_MAIN()
 !###======================================================================
 IMPLICIT NONE
 TYPE(WIN_MESSAGE) :: MESSAGE
 INTEGER :: ITYPE,ICOL,IROW,IRGB
 
 CALL WDIALOGSELECT(ID_D3DSETTINGS_PART)
 
 ITYPE=0; IF(PL%IRUN.EQ.0)ITYPE=1
 CALL WDIALOGFIELDSTATE(ID_OPEN,ITYPE)
 CALL WDIALOGFIELDSTATE(ID_SAVEAS,ITYPE)
 CALL WDIALOGFIELDSTATE(ID_DELETE,ITYPE)
 
 CALL TRACE_3D_STARTPOINTS_GRID()
 
 CALL UTL_DIALOGSHOW(-1,-1,0,3)
 DO
  CALL WMESSAGE(ITYPE,MESSAGE)
  SELECT CASE (ITYPE)
   CASE (FIELDCHANGED)
    SELECT CASE (MESSAGE%VALUE2)
     CASE (IDF_GRID1)
      CALL WGRIDPOS(MESSAGE%Y,ICOL,IROW)
      IF(ICOL.EQ.2)THEN
       IRGB=SP(IROW)%ICLR; CALL WSELECTCOLOUR(IRGB)
       IF(WINFODIALOG(4).EQ.1)THEN
        CALL WGRIDCOLOURCELL(IDF_GRID1,2,IROW,IRGB,IRGB)
        CALL WGRIDPUTCELLINTEGER(IDF_GRID1,2,IROW,IRGB)
       ENDIF
       CALL WGRIDSETCELL(IDF_GRID1,1,IROW)
      ENDIF
    END SELECT
   CASE (PUSHBUTTON)
    SELECT CASE (MESSAGE%VALUE1)
     CASE (ID_OPEN)
      CALL TRACE_3D_STARTPOINTS_SAVE(MESSAGE%VALUE1)
      CALL TRACE_3D_STARTPOINTS_GRID()
     CASE (ID_SAVEAS)
      CALL WGRIDGETCHECKBOX(IDF_GRID1,1,SP%IACT ,NSPG)
      CALL WGRIDGETINTEGER (IDF_GRID1,2,SP%ICLR ,NSPG)
      CALL WGRIDGETMENU    (IDF_GRID1,4,SP%IREV ,NSPG); SP%IREV=SP%IREV-1
      CALL WGRIDGETDOUBLE  (IDF_GRID1,5,SP%SPWIDTH,NSPG)
      CALL WGRIDGETDOUBLE  (IDF_GRID1,6,SP%PWIDTH,NSPG)
      CALL WGRIDGETINTEGER (IDF_GRID1,7,SP%IRSTRT,NSPG); SP%IRSTRT=MAX(SP%IRSTRT,0)
      CALL TRACE_3D_STARTPOINTS_SAVE(MESSAGE%VALUE1)
     CASE (IDOK)
      CALL WGRIDGETCHECKBOX(IDF_GRID1,1,SP%IACT ,NSPG)
      CALL WGRIDGETINTEGER (IDF_GRID1,2,SP%ICLR ,NSPG)
      CALL WGRIDGETMENU    (IDF_GRID1,4,SP%IREV ,NSPG); SP%IREV=SP%IREV-1
      CALL WGRIDGETDOUBLE    (IDF_GRID1,5,SP%SPWIDTH,NSPG)
      CALL WGRIDGETDOUBLE    (IDF_GRID1,6,SP%PWIDTH,NSPG)
      CALL WGRIDGETINTEGER (IDF_GRID1,7,SP%IRSTRT,NSPG); SP%IRSTRT=MAX(SP%IRSTRT,0)
      EXIT
     CASE (IDCANCEL)
      EXIT
    END SELECT
  END SELECT
 ENDDO
 CALL WDIALOGHIDE(); CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 
 !## recompute start interval particles
 CALL TRACE_3D_STARTPOINTS_STARTINTERVAL()
 
 IF(NSPG.GT.SIZE(SP))THEN
  CALL WDIALOGFIELDSTATE(ID_ADD,0)
 ELSE
  !## only if no particle-simulation is running
  IF(PL%IRUN.EQ.0)CALL WDIALOGFIELDSTATE(ID_ADD,1)
 ENDIF
 
 END SUBROUTINE TRACE_3D_STARTPOINTS_MAIN
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_GRID()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I
 
 CALL WDIALOGSELECT(ID_D3DSETTINGS_PART)
 CALL WGRIDCLEAR(IDF_GRID1)

 IF(NSPG.GT.0)THEN
 
  CALL WGRIDROWS(IDF_GRID1,NSPG)
  CALL WDIALOGFIELDSTATE(IDF_GRID1,1)
  !## fill in dialog
  CALL WGRIDPUTCHECKBOX(IDF_GRID1,1,SP%IACT,NSPG)
  !## number of particles in each group
  CALL WGRIDPUTINTEGER(IDF_GRID1,3,SP%NPART,NSPG)
  CALL WGRIDSTATE(IDF_GRID1,3,2)
  !## assign current colour to group
  CALL WGRIDPUTINTEGER(IDF_GRID1,2,SP%ICLR,NSPG)
  DO I=1,NSPG; CALL WGRIDCOLOURCELL(IDF_GRID1,2,I,SP(I)%ICLR,SP(I)%ICLR); ENDDO
  !## direction of particle
  SP%IREV=SP%IREV+1; CALL WGRIDPUTOPTION(IDF_GRID1,4,SP%IREV,NSPG); SP%IREV=SP%IREV-1
  !## size of the startpoints
  CALL WGRIDPUTDOUBLE(IDF_GRID1,5,SP%SPWIDTH,NSPG)
  !## size of the particles
  CALL WGRIDPUTDOUBLE(IDF_GRID1,6,SP%PWIDTH,NSPG)
  !## start interval particles
  CALL WGRIDPUTINTEGER(IDF_GRID1,7,SP%IRSTRT,NSPG)
 
 ELSE
 
  CALL WGRIDROWS(IDF_GRID1,1)
  CALL WDIALOGFIELDSTATE(IDF_GRID1,2)
 
 ENDIF

 END SUBROUTINE TRACE_3D_STARTPOINTS_GRID
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_SAVE(ID)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ID
 CHARACTER(LEN=256) :: FNAME
 INTEGER :: IU,I,J,IG2
 LOGICAL :: LIPF
 
 IU=UTL_GETUNIT()
 
 IF(ID.EQ.ID_OPEN)THEN
 
  !## all existing particles will be removed
  IF(PL%NPART.GT.0)THEN
   CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONCANCEL,'Do you want to remove existing particles first?','Question')
   IF(WINFODIALOG(4).NE.1)RETURN
  ENDIF
  IF(.NOT.UTL_WSELECTFILE('iMOD Particle File (*.ptf)|*.ptf|',&
           LOADDIALOG+MUSTEXIST+PROMPTON+DIRCHANGE+APPENDEXT,FNAME,'Load iMOD Particle File (*.ptf)'))RETURN
  CALL OSD_OPEN(IU,FILE=FNAME,STATUS='OLD',ACTION='READ')
 
  CALL TRACE_3D_STARTPOINTS_RESET()

  READ(IU,'(I10)') NSPG

  DO I=1,MIN(SIZE(SP),NSPG)
   READ(IU,'(4(I10,1X),2(F10.2,1X),I10)') SP(I)%NPART,SP(I)%IACT,SP(I)%ICLR,SP(I)%IREV,SP(I)%SPWIDTH,SP(I)%PWIDTH,SP(I)%IRSTRT
   SP(I)%IRSTRT=MAX(SP(I)%IRSTRT,0)
   ALLOCATE(XSP(SP(I)%NPART),YSP(SP(I)%NPART),ZSP(SP(I)%NPART))
   DO J=1,SP(I)%NPART; READ(IU,*) XSP(J),YSP(J),ZSP(J); ENDDO
   !## get new location/particle - click from the 3d tool
   IG2=SP(I)%NPART
   CALL TRACE_3D_STARTPOINTS_ASSIGN(I,1,IG2)
   DEALLOCATE(XSP,YSP,ZSP)
   !## total particles for this group
   PL%NPART=PL%NPART+SP(I)%NPART
  ENDDO
  
  CLOSE(IU)
  CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'Succesfully read '//TRIM(ITOS(PL%NPART))//' particles from'//CHAR(13)//TRIM(FNAME),'Information')

 ELSE
 
  IF(.NOT.UTL_WSELECTFILE('iMOD Map (*.ipf)|*.ipf|iMOD Particle File (*.ptf)|*.ptf|',&
           SAVEDIALOG+DIRCHANGE+APPENDEXT,FNAME,'Save iMOD Particle File (*.ipf,*.ptf)'))RETURN
  CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',ACTION='WRITE')
  LIPF=.FALSE.; IF(INDEX(UTL_CAP(FNAME,'U'),'.IPF').GT.0)LIPF=.TRUE.
  IF(LIPF)THEN
   WRITE(IU,'(I10)') SUM(SP(1:NSPG)%NPART)
   WRITE(IU,'(A)') '4'; WRITE(IU,'(A)') 'X'; WRITE(IU,'(A)') 'Y'; WRITE(IU,'(A)') 'Z'; WRITE(IU,'(A)') 'GROUP'; WRITE(IU,'(A)') '0,TXT'
   DO I=1,NSPG
    DO J=1,SP(I)%NPART
     WRITE(IU,'(2(F10.2,1X),F10.3,1X,I10)') IDF%XMIN+SP(I)%XLC(J),IDF%YMIN+SP(I)%YLC(J),SP(I)%ZLC(J),I
    ENDDO
   ENDDO
  ELSE
   WRITE(IU,'(I10)') NSPG
   DO I=1,NSPG
    WRITE(IU,'(4(I10,1X),2(F10.2,1X),I10)') SP(I)%NPART,SP(I)%IACT,SP(I)%ICLR,SP(I)%IREV,SP(I)%SPWIDTH,SP(I)%PWIDTH,SP(I)%IRSTRT
    DO J=1,SP(I)%NPART
     WRITE(IU,'(2(F10.2,1X),F10.3)') IDF%XMIN+SP(I)%XLC(J),IDF%YMIN+SP(I)%YLC(J),SP(I)%ZLC(J)
    ENDDO
   ENDDO
  ENDIF
  CLOSE(IU)
  CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'Succesfully saved particles to'//CHAR(13)//TRIM(FNAME),'Information')
 ENDIF
 
 END SUBROUTINE TRACE_3D_STARTPOINTS_SAVE
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS() 
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,ISP,NP,NG,IG,IG1,IG2
 
 !## get location of mouse
 CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO5,ISP)

 ALLOCATE(NSGRP(NLAY)); NSGRP=0
 ALLOCATE(ILAYERS(NLAY),JLAYERS(NLAY)); ILAYERS=0; JLAYERS=0
 
 SELECT CASE (ISP)
  !## spatial distributed
  CASE (1)
   CALL TRACE_3D_STARTPOINTS_LOCATION(NSGRP,NG)
  !## spatial distributed via an IDF (2) or an external IPF/IDF-file (5)
  CASE (2,5)
   CALL TRACE_3D_STARTPOINTS_IDF(NSGRP,ISP,NG)
  !## sinks
  CASE (3)
   CALL TRACE_3D_STARTPOINTS_SINK(NSGRP,NG)
  !## border - layers
  CASE (4)
   CALL TRACE_3D_STARTPOINTS_LAYER(NSGRP,NG)
 END SELECT
 
 !## check if really many particles are supposed to be added
 IF(SUM(NSGRP).GT.10000)THEN
  CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'Are you sure to add '//TRIM(ITOS(SUM(NSGRP)))//' particles','Question')
  IF(WINFODIALOG(4).NE.1)NSGRP=-1
 ENDIF
 
 IF(SUM(NSGRP).GT.0)THEN
  
  IF(NSPG+NG.GT.SIZE(SP))THEN
   CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'Total groups '//TRIM(ITOS(NSPG+NG))//' overwrite'//CHAR(13)// &
      'maximal groups sustained ('//TRIM(ITOS(SIZE(SP)))//')'//CHAR(13)//CHAR(13)//'Continue?','Question')
   IF(WINFODIALOG(4).NE.1)NG=0
  ENDIF
  
  DO IG=1,NG
   IG1=1; IF(IG.GT.1)IG1=NSGRP(IG-1)+1; IG2=NSGRP(IG); NP=(IG2-IG1)+1   

   i=nspg
   i=i+1
   nspg=i
! write(*,*) nspg
!   NSPG=NSPG+1
! write(*,*) nspg
   
   IF(NSPG.GT.SIZE(SP))EXIT
   
   !## get new location/particle - click from the 3d tool
   CALL TRACE_3D_STARTPOINTS_ASSIGN(NSPG,IG1,IG2)
   !## total particles for this group
   PL%NPART=PL%NPART+SP(NSPG)%NPART
   !## assign current colour to group
   IF(IG.EQ.1)THEN
    SP(NSPG)%ICLR=PL%SPCOLOR
   ELSE
    !## apply other colour for rest of group
    I=MOD(NSPG,MAXCOLOUR); I=MAX(1,I)
    SP(NSPG)%ICLR=ICOLOR(I)
   ENDIF
   !## activate current group
   SP(NSPG)%IACT=1
   !## flow direction current group = current 
   SP(NSPG)%IREV=PL%IREV
   !## plot size
   SP(NSPG)%SPWIDTH=3.0
   SP(NSPG)%PWIDTH=1.0D0
   !## start times directly
   SP(NSPG)%IRSTRT=0
  
  ENDDO
  
 ELSE
  IF(SUM(NSGRP).EQ.0)CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'iMOD cannot place particles for current configuration.','Information')
 ENDIF
 
 !## deallocate memory for temporary storage of locations
 IF(ALLOCATED(XSP))    DEALLOCATE(XSP)
 IF(ALLOCATED(YSP))    DEALLOCATE(YSP)
 IF(ALLOCATED(ZSP))    DEALLOCATE(ZSP)
 IF(ALLOCATED(NSGRP))  DEALLOCATE(NSGRP)
 IF(ALLOCATED(ILAYERS))DEALLOCATE(ILAYERS)
 IF(ALLOCATED(JLAYERS))DEALLOCATE(JLAYERS)

 CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 IF(NSPG.GE.SIZE(SP))CALL WDIALOGFIELDSTATE(ID_ADD,0)
 NSPG=MIN(NSPG,SIZE(SP))
 
 END SUBROUTINE TRACE_3D_STARTPOINTS

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_SINK(NSGRP,NG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(OUT) :: NG
 INTEGER,DIMENSION(:),INTENT(INOUT) :: NSGRP
 INTEGER,DIMENSION(2) :: IOPT
 INTEGER :: I,II,J,JJ,NC,N,IRANDOM,ISUBSQ,ISGRP,ISHAPE, &
            INCLD,IL1,IL2,ILAY,IROW,ICOL,ISTRONG,NZ,JLAY
 REAL(KIND=DP_KIND) :: R,G2R,Q,OR,QERROR,XC,YC,DZ,X,Y,Z,OFR
 LOGICAL,DIMENSION(2) :: LEX

 CALL WDIALOGGETDOUBLE(IDF_REAL10,Q)
 CALL WDIALOGGETMENU(IDF_MENU2,IOPT(1))
 CALL WDIALOGGETCHECKBOX(IDF_CHECK5,ISTRONG)
   
 CALL WDIALOGSELECT(ID_D3DSETTINGS_SINKS)
 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO9,ISHAPE) !## square/circle 
 CALL WDIALOGGETINTEGER(IDF_INTEGER6,NC)       !## number on circle
 CALL WDIALOGGETINTEGER(IDF_INTEGER7,NZ)       !## number vertially
 CALL WDIALOGGETDOUBLE(IDF_REAL14,R)             !## radius
 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,IRANDOM)   !## randomize positions
 CALL WDIALOGGETMENU(IDF_MENU4,INCLD)
 CALL WDIALOGGETMENU(IDF_MENU3,IOPT(2))

 CALL WDIALOGGETMENU(IDF_MENU1,ILAYERS)
 DO I=1,NLAY; IF(ILAYERS(I).EQ.1)ILAYERS(I)=I; ENDDO; JLAYERS=ILAYERS
 CALL WDIALOGGETCHECKBOX(IDF_CHECK2,ISUBSQ)    !## subsequent layers
 CALL WDIALOGGETCHECKBOX(IDF_CHECK3,ISGRP)     !## new group per (group of) layer(s)

  IF(ISUBSQ.EQ.1)CALL TRACE_3D_STARTPOINTS_SUBQ()
      
 !## from degrees to radians
 G2R=360.0D0/(2.0*PI); OR=360.0D0/REAL(NC); OR=OR/G2R

 DO I=1,2
  N=0; NG=0
  DO JLAY=1,NLAY 
    
   IF(ILAYERS(JLAY).EQ.0.OR.JLAYERS(JLAY).EQ.0)CYCLE
   IL1=ILAYERS(JLAY); IL2=JLAYERS(JLAY)

   DO ILAY=IL1,IL2; DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
       
    !## correct, proven otherwise
    LEX=.TRUE.

    !## check whether strong-sink
    IF(ISTRONG.EQ.1)THEN
     !## forward simulation
     IF(PL%IREV.EQ.0)THEN
      IF(QX(ICOL,IROW,ILAY).LT.0.0D0.OR.QX(ICOL+1,IROW,ILAY).GT.0.0D0.OR. &
         QY(ICOL,IROW,ILAY).GT.0.0D0.OR.QY(ICOL,IROW+1,ILAY).LT.0.0D0.OR. &
         QZ(ICOL,IROW,ILAY).GT.0.0D0.OR.QZ(ICOL,IROW,ILAY+1).LT.0.0D0)LEX(1)=.FALSE.
     !## backward simulation
     ELSE
      IF(QX(ICOL,IROW,ILAY).GT.0.0D0.OR.QX(ICOL+1,IROW,ILAY).LT.0.0D0.OR. &
         QY(ICOL,IROW,ILAY).LT.0.0D0.OR.QY(ICOL,IROW+1,ILAY).GT.0.0D0.OR. &
         QZ(ICOL,IROW,ILAY).LT.0.0D0.OR.QZ(ICOL,IROW,ILAY+1).GT.0.0D0)LEX(1)=.FALSE.
     ENDIF
    ENDIF
    !## no a strong sink
    IF(.NOT.LEX(1))CYCLE
      
    !## net inflow is negative
    QERROR=QSS(ICOL,IROW,ILAY) 

    !## not correct, proven otherwise
    LEX=.FALSE.
    DO J=1,2
     !## forward simulation
     SELECT CASE (IOPT(J))
      CASE (1) !## =
       IF(QERROR.EQ.Q)LEX(J)=.TRUE.
      CASE (2) !## >
       IF(QERROR.GT.Q)LEX(J)=.TRUE.
      CASE (3) !## <
       IF(QERROR.LT.Q)LEX(J)=.TRUE.
      CASE (4) !## >=
       IF(QERROR.GE.Q)LEX(J)=.TRUE.
      CASE (5) !## <=
       IF(QERROR.LE.Q)LEX(J)=.TRUE.
     END SELECT
    ENDDO
    !## and
    IF(INCLD.EQ.2)LEX(1)=LEX(1).AND.LEX(2)
    !## or
    IF(INCLD.EQ.3)LEX(1)=LEX(1).OR.LEX(2)
          
    !## try next
    IF(.NOT.LEX(1))CYCLE
   
    CALL IDFGETLOC(IDF,IROW,ICOL,XC,YC)
    DZ=ZTOP(ICOL,IROW,ILAY)-ZBOT(ICOL,IROW,ILAY)
    IF(IRANDOM.EQ.0)DZ=DZ/REAL(NZ)
     
    DO II=1,NC  !## number on circle      
     !## get coordinates
     IF(IRANDOM.EQ.0)THEN
      X=XC+R*COS(REAL(II)*OR); Y=YC+R*SIN(REAL(II)*OR)
      Z=ZTOP(ICOL,IROW,ILAY)+(DZ*0.5)
     ENDIF
     DO JJ=1,NZ !## number vertically
      N=N+1
      IF(IRANDOM.EQ.0)THEN
       Z=Z-DZ
      ELSE
       CALL RANDOM_NUMBER(X); OFR=(X*360.0D0)/G2R; X=XC+R*COS(OFR); Y=YC+R*SIN(OFR)
       CALL RANDOM_NUMBER(Z); Z=Z*DZ;            Z=Z +ZBOT(ICOL,IROW,ILAY)
      ENDIF
      IF(I.EQ.2)THEN; XSP(N)=X; YSP(N)=Y; ZSP(N)=Z; ENDIF
     ENDDO
    ENDDO
    
   ENDDO; ENDDO; ENDDO

   IF(NG.EQ.0)THEN
    NG=NG+1
   ELSE
    NG=NG+ISGRP
   ENDIF
   NSGRP(NG)=N
    
   ENDDO
  IF(I.EQ.1)ALLOCATE(XSP(N),YSP(N),ZSP(N))
 ENDDO
 
 END SUBROUTINE TRACE_3D_STARTPOINTS_SINK
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_LAYER(NSGRP,NG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(OUT) :: NG
 INTEGER,DIMENSION(:),INTENT(INOUT) :: NSGRP
 INTEGER,DIMENSION(2) :: IOPT
 INTEGER :: I,II,J,NZ,NX,N,IRANDOM,ISUBSQ,ISGRP,IC1,IC2,IR1,IR2,DC,DR,IL1,IL2,ILAY,IROW,ICOL
 REAL(KIND=DP_KIND) :: DZ,XC,YC,Z,X
 
 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO12,IOPT(1))
 CALL WDIALOGSELECT(ID_D3DSETTINGS_LAYER)
 CALL WDIALOGGETINTEGER(IDF_INTEGER8,NZ)
 CALL WDIALOGGETINTEGER(IDF_INTEGER1,NX)
 CALL WDIALOGGETMENU(IDF_MENU1,ILAYERS)
 DO I=1,NLAY; IF(ILAYERS(I).EQ.1)ILAYERS(I)=I; ENDDO; JLAYERS=ILAYERS
 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,IRANDOM)
 CALL WDIALOGGETCHECKBOX(IDF_CHECK2,ISUBSQ)
 CALL WDIALOGGETCHECKBOX(IDF_CHECK3,ISGRP)

 IF(ISUBSQ.EQ.1)CALL TRACE_3D_STARTPOINTS_SUBQ() 

 !## iopt=1 west
 !## iopt=2 north
 !## iopt=3 east
 !## iopt=4 south
 !## iopt=5 all
   
 IC1=1; IC2=IDF%NCOL; IR1=1; IR2=IDF%NROW; DC=NX; DR=NX
 
 !## east
 IF(IOPT(1).EQ.3)THEN; DC=-NX; IC1=IDF%NCOL; IC2=1; ENDIF
 !## south
 IF(IOPT(1).EQ.4)THEN; DR=-NX; IR1=IDF%NROW; IR2=1; ENDIF

 IF(IOPT(1).EQ.5)THEN
  N=NZ*IC2*IR2
 ELSE
  N=NZ*MAX(IC2,IR2)
 ENDIF
 NX=0; DO I=1,NLAY; IF(JLAYERS(I).GT.0)NX=NX+N; ENDDO; N=NX
 ALLOCATE(XSP(N),YSP(N),ZSP(N))

 N=0; NG=0
 DO ILAY=1,NLAY
  IF(ILAYERS(ILAY).EQ.0.OR.JLAYERS(ILAY).EQ.0)CYCLE
  IL1=ILAYERS(ILAY); IL2=JLAYERS(ILAY)
  SELECT CASE (IOPT(1))
   CASE (1,3) !## west/east
    DO IROW=IR1,IR2,DR; DO ICOL=IC1,IC2,DC
     IF(IBOUND(ICOL,IROW,ILAY).EQ.0)CYCLE
     DZ=ZTOP(ICOL,IROW,IL1)-ZBOT(ICOL,IROW,IL2)
     IF(IRANDOM.EQ.0)DZ=DZ/REAL(NZ)
     Z = ZTOP(ICOL,IROW,IL1)+(0.5*DZ)
     CALL IDFGETLOC(IDF,IROW,ICOL,XC,YC)
     DO J=1,NZ
      N=N+1
      IF(IRANDOM.EQ.0)Z=Z-DZ
      IF(IRANDOM.EQ.1)THEN; CALL RANDOM_NUMBER(Z); Z=ZTOP(ICOL,IROW,IL1)-Z*DZ; ENDIF
      XSP(N)=XC; YSP(N)=YC; ZSP(N)=Z
     ENDDO
     EXIT
    ENDDO; ENDDO
   CASE (2,4) !## north/south
    DO ICOL=IC1,IC2,DC; DO IROW=IR1,IR2,DR
     IF(IBOUND(ICOL,IROW,ILAY).EQ.0)CYCLE
     DZ=ZTOP(ICOL,IROW,IL1)-ZBOT(ICOL,IROW,IL2)
     IF(IRANDOM.EQ.0)DZ=DZ/REAL(NZ)
     Z = ZTOP(ICOL,IROW,IL1)+(0.5*DZ)
     CALL IDFGETLOC(IDF,IROW,ICOL,XC,YC)
     DO J=1,NZ
      N=N+1
      IF(IRANDOM.EQ.0)Z=Z-DZ
      IF(IRANDOM.EQ.1)THEN; CALL RANDOM_NUMBER(Z); Z=ZTOP(ICOL,IROW,IL1)-Z*DZ; ENDIF
      XSP(N)=XC; YSP(N)=YC; ZSP(N)=Z
     ENDDO
     EXIT
    ENDDO; ENDDO
   CASE (5) !## all
    DO I=IC1,IC2,DC; DO J=IR1,IR2,DR
     IF(IRANDOM.EQ.0)THEN
      ICOL=I; IROW=J
     ELSE
      CALL RANDOM_NUMBER(X); ICOL=MAX(1.0D0,IDF%NCOL*X)
      CALL RANDOM_NUMBER(X); IROW=MAX(1.0D0,IDF%NROW*X)
     ENDIF
     IF(IBOUND(ICOL,IROW,ILAY).EQ.0)CYCLE
     DZ=ZTOP(ICOL,IROW,IL1)-ZBOT(ICOL,IROW,IL2)
     IF(IRANDOM.EQ.0)DZ=DZ/REAL(NZ)
     Z = ZTOP(ICOL,IROW,IL1)+(0.5*DZ)
     CALL IDFGETLOC(IDF,IROW,ICOL,XC,YC)
     DO II=1,NZ
      N=N+1
      IF(IRANDOM.EQ.0)Z=Z-DZ
      IF(IRANDOM.EQ.1)THEN; CALL RANDOM_NUMBER(Z); Z=ZTOP(ICOL,IROW,IL1)-Z*DZ; ENDIF
      XSP(N)=XC; YSP(N)=YC; ZSP(N)=Z
     ENDDO
    ENDDO; ENDDO
  END SELECT
  IF(NG.EQ.0)THEN
   NG=NG+1
  ELSE
   NG=NG+ISGRP
  ENDIF
  NSGRP(NG)=N 
 ENDDO
 
 END SUBROUTINE TRACE_3D_STARTPOINTS_LAYER

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_IDF(NSGRP,ISP,NG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(OUT) :: NG
 INTEGER,INTENT(IN) :: ISP
 INTEGER,DIMENSION(:),INTENT(INOUT) :: NSGRP
 LOGICAL :: LIDF,LIPF
 LOGICAL,DIMENSION(2) :: LEX
 INTEGER :: I,J,II,IROW,ICOL,IIDF,INCLD,NX,JLAY,IRANDOM,N,M
 REAL(KIND=DP_KIND) :: DZ,X
 TYPE(IDFOBJ) :: PIDF
 INTEGER,ALLOCATABLE,DIMENSION(:) :: ISC,ISR
 
 !## file
 LIDF=.FALSE.; LIPF=.FALSE.
 IF(ISP.EQ.2)THEN; CALL WDIALOGGETMENU(IDF_MENU1,IIDF); LIDF=.TRUE.; ENDIF
 IF(ISP.EQ.5)THEN
  CALL WDIALOGGETSTRING(IDF_STRING2,PIDF%FNAME)
  IF(INDEX(UTL_CAP(PIDF%FNAME,'U'),'.IDF').GT.0)LIDF=.TRUE.
  IF(INDEX(UTL_CAP(PIDF%FNAME,'U'),'.IPF').GT.0)LIPF=.TRUE.
 ENDIF
 
 CALL WDIALOGSELECT(ID_D3DSETTINGS_IDF)
 CALL WDIALOGGETMENU(IDF_MENU1,INCLD)
 CALL WDIALOGGETINTEGER(IDF_INTEGER1,NX)
 CALL WDIALOGGETINTEGER(IDF_INTEGER2,JLAY)   !## layer for vertical fluxes
 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,IRANDOM)
 CALL WDIALOGGETDOUBLE(IDF_REAL1,DZ)

 N=0
 IF(LIDF)THEN
  CALL IDFCOPY(IDF,PIDF); IF(.NOT.IDFALLOCATEX(PIDF))THEN; ENDIF
  CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
  CALL WDIALOGGETSTRING(IDF_STRING2,PIDF%FNAME)
  IF(ISP.EQ.2)THEN
   IF(.NOT.IDFREADSCALE(IDFPLOT(IIDF)%FNAME,PIDF,2,0,0.0D0,0))RETURN
  ELSEIF(ISP.EQ.5)THEN
   IF(.NOT.IDFREADSCALE(PIDF%FNAME,PIDF,2,0,0.0D0,0))RETURN
  ENDIF

  !## store unique number of col/row
  DO II=1,2
   M=0
   DO I=1,PIDF%NROW,NX; DO J=1,PIDF%NCOL,NX
    M=M+1
    IF(II.EQ.2)THEN
     IF(IRANDOM.EQ.0)THEN
      ICOL=J; IROW=I
     ELSE
      CALL RANDOM_NUMBER(X); ICOL=MAX(1.0D0,PIDF%NCOL*X)
      CALL RANDOM_NUMBER(X); IROW=MAX(1.0D0,PIDF%NROW*X)
     ENDIF
     ISC(M)=ICOL; ISR(M)=IROW
    ENDIF
   ENDDO; ENDDO
   IF(II.EQ.1)ALLOCATE(ISC(M),ISR(M))
  ENDDO
  
  DO II=1,2
   N=0; M=0
   DO I=1,PIDF%NROW,NX; DO J=1,PIDF%NCOL,NX
    M=M+1; ICOL=ISC(M); IROW=ISR(M)
    IF(PIDF%X(ICOL,IROW).NE.PIDF%NODATA)THEN
     IF(ISP.EQ.2)THEN
      !## take fluxes layer beneath as modpath stores the fluxes like that
      X=QZ(ICOL,IROW,JLAY+1)  
     ELSE
      X=PIDF%X(ICOL,IROW)
     ENDIF
     LEX(1)=.FALSE.
     SELECT CASE (INCLD)
      CASE (1) !## both (-/+)
       LEX(1)=.TRUE.
      CASE (2) !## down (-)
       IF(X.GT.0.0D0)LEX(1)=.TRUE.
      CASE (3) !## up (+)
       IF(X.LT.0.0D0)LEX(1)=.TRUE.
     END SELECT
     IF(LEX(1))THEN
      N=N+1
      IF(II.EQ.2)THEN
       CALL IDFGETLOC(PIDF,IROW,ICOL,XSP(N),YSP(N))
       ZSP(N)=PIDF%X(ICOL,IROW)+DZ
      ENDIF
     ENDIF
    ENDIF
   ENDDO; ENDDO
   IF(II.EQ.1)ALLOCATE(XSP(N),YSP(N),ZSP(N))
  ENDDO
  CALL IDFDEALLOCATEX(PIDF); DEALLOCATE(ISC,ISR)
 ENDIF
 
 IF(LIPF)THEN
  NIPF=1; CALL IPFALLOCATE(); IPFFILE(1)%FNAME=PIDF%FNAME
  !## x,y,z
  IPFFILE(1)%XCOL=1; IPFFILE(1)%YCOL=2; IPFFILE(1)%ZCOL=3; IPFFILE(1)%Z2COL=4; IPFFILE(1)%QCOL=1
  IF(.NOT.IPFREAD2(1,1,0))RETURN
  N=IPFFILE(1)%NROW; ALLOCATE(XSP(N),YSP(N),ZSP(N))
  DO I=1,N
   XSP(I)=IPFFILE(1)%XYZ(1,I); YSP(I)=IPFFILE(1)%XYZ(2,I); ZSP(I)=IPFFILE(1)%XYZ(3,I) 
  ENDDO
  CALL IPFDEALLOCATE()
 ENDIF
 
 NG=1; NSGRP(NG)=N
 
 END SUBROUTINE TRACE_3D_STARTPOINTS_IDF
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_LOCATION(NSGRP,NG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(OUT) :: NG
 INTEGER,DIMENSION(:),INTENT(INOUT) :: NSGRP
 REAL(KIND=DP_KIND) :: XC,YC,ZC,X,Y,Z,DX,DY,DZ
 INTEGER :: NX,NZ,I,J,K,N
 
 NSGRP=0
 
 CALL WDIALOGSELECT(ID_D3DSETTINGS_LOCATION)
 CALL WDIALOGGETDOUBLE(IDF_REAL1,XC)
 CALL WDIALOGGETDOUBLE(IDF_REAL2,YC)
 CALL WDIALOGGETDOUBLE(IDF_REAL3,ZC)
 CALL WDIALOGGETDOUBLE(IDF_REAL11,X)
 CALL WDIALOGGETDOUBLE(IDF_REAL12,Y)
 CALL WDIALOGGETDOUBLE(IDF_REAL13,Z)
 DX=X-XC; DY=Y-YC; DZ=Z-ZC
 CALL WDIALOGGETINTEGER(IDF_INTEGER3,NX)
 CALL WDIALOGGETINTEGER(IDF_INTEGER2,NZ)

 !## allocate memory for temporary storage of locations
 N=NZ*NX**2; ALLOCATE(XSP(N),YSP(N),ZSP(N))
 !## start at lower-left/bottom-corner
 XC=XC-0.5*DX; YC=YC-0.5*DX; ZC=ZC-0.5*DZ
 DX=DX/REAL(NX); DZ=DZ/REAL(NZ)

 N=0; Z=ZC-DZ
 DO K=1,NZ
  Z=Z+DZ; Y=YC-DX
  DO I=1,NX
   Y=Y+DX; X=XC-DX
   DO J=1,NX
    X=X+DX
    N=N+1; XSP(N)=X; YSP(N)=Y; ZSP(N)=Z
   ENDDO
  ENDDO
 ENDDO
 NG=1; NSGRP(NG)=N

 END SUBROUTINE TRACE_3D_STARTPOINTS_LOCATION
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_SUBQ() 
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J
 
 JLAYERS=0; J=0
 DO I=1,NLAY
  IF(ILAYERS(I).NE.0)THEN
   IF(J.EQ.0)J=I
  ELSE
   IF(J.NE.0)JLAYERS(J)=I-1; J=0
  ENDIF
 ENDDO
 IF(J.NE.0)JLAYERS(J)=I-1
   
 END SUBROUTINE TRACE_3D_STARTPOINTS_SUBQ

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_ASSIGN(IG,IG1,IG2) 
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IG,IG1,IG2
 INTEGER :: IROW,ICOL,ILAY,I,IZL
 REAL(KIND=DP_KIND) :: XC,YC,ZC,ZL,DZ

 SP(IG)%NPART=0
 
 DO I=IG1,IG2 

  XC=XSP(I); YC=YSP(I); ZC=ZSP(I)
  
  !## location outside model dimensions - skip
  IF(XC.LE.IDF%XMIN.OR.XC.GE.IDF%XMAX.OR. &
     YC.LE.IDF%YMIN.OR.YC.GE.IDF%YMAX)CYCLE

  !## get current irow/icol for model-idf
  CALL IDFIROWICOL(IDF,IROW,ICOL,XC,YC)
 
  !## particle too high - above surfacelevel
  IF(ZC.GT.ZTOP(ICOL,IROW,1))CYCLE
  !## particle too low - below lowest bottom
  IF(ZC.LT.ZBOT(ICOL,IROW,NLAY))CYCLE
    
  !## determine in which layer location is situated, skip positions within clay-layers
  DO ILAY=1,NLAY

   !## skip inactive cells (constant head are strong sinks!)
   IF(IBOUND(ICOL,IROW,ILAY).EQ.0)CYCLE
   
   !## no assigned yet
   IZL=0
   ZL =0.0D0
   
   !## inside current modellayer
   IF(ZC.LE.ZTOP(ICOL,IROW,ILAY).AND.ZC.GE.ZBOT(ICOL,IROW,ILAY))THEN
    !## compute local z: top (zl=1); mid (zl=0.5); bot (zl=0)
    DZ=ZTOP(ICOL,IROW,ILAY)-ZBOT(ICOL,IROW,ILAY)
    IF(DZ.NE.0.0D0)ZL=(ZC-ZBOT(ICOL,IROW,ILAY))/DZ
    IZL=1
   ENDIF
   IF(ZL.EQ.0.0D0.AND.ILAY.LT.NLAY)THEN
    !## skip inactive/constant head cells under layer
    IF(IBOUND(ICOL,IROW,ILAY+1).EQ.0)CYCLE
    !## inside interbed between modellayers
    IF(ZC.LT.ZBOT(ICOL,IROW,ILAY).AND.ZC.GT.ZTOP(ICOL,IROW,ILAY+1))THEN
     !## compute local z: top (zl=1); mid (zl=0.5); bot (zl=0)
     DZ=ZBOT(ICOL,IROW,ILAY)-ZTOP(ICOL,IROW,ILAY+1)
     IF(DZ.NE.0.0D0)ZL=(ZTOP(ICOL,IROW,ILAY+1)-ZC)/DZ
     IZL=1
    ENDIF
   ENDIF

   IF(IZL.EQ.1)THEN 

    !## count number of particles
    SP(IG)%NPART=SP(IG)%NPART+1
    
    CALL TRACE_3D_STARTPOINTS_MEMORY_SP(IG)  

    SP(IG)%XLC(SP(IG)%NPART)=XC-IDF%XMIN
    SP(IG)%YLC(SP(IG)%NPART)=YC-IDF%YMIN
    SP(IG)%ZLC(SP(IG)%NPART)=ZC
    SP(IG)%ZLL(SP(IG)%NPART)=ZL
    SP(IG)%KLC(SP(IG)%NPART)=ILAY
    SP(IG)%ILC(SP(IG)%NPART)=IROW
    SP(IG)%JLC(SP(IG)%NPART)=ICOL
    SP(IG)%TOT(SP(IG)%NPART)=0.0D0
    SP(IG)%MXL(SP(IG)%NPART)=0
    SP(IG)%STM(SP(IG)%NPART)=0
    EXIT
    
   ENDIF
  ENDDO

 ENDDO

 END SUBROUTINE TRACE_3D_STARTPOINTS_ASSIGN

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_STARTINTERVAL()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J
 REAL(KIND=DP_KIND) :: X

 DO I=1,NSPG
  IF(SP(I)%IRSTRT.EQ.0)CYCLE
  DO J=1,SP(I)%NPART
   !## random start squence number - negative to indicate it did not start yet
   CALL RANDOM_NUMBER(X)
   SP(I)%STM(J)=X*SP(I)%IRSTRT
  ENDDO   
 ENDDO

 END SUBROUTINE TRACE_3D_STARTPOINTS_STARTINTERVAL

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_MEMORY_SP(IG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IG
 INTEGER :: I,N
 
 IF(SP(IG)%NPART.LE.SIZE(SP(IG)%XLC))RETURN
 
 N=SP(IG)%NPART+999
 
 ALLOCATE(SP(IG)%XLC_BU(N),SP(IG)%YLC_BU(N),SP(IG)%ZLC_BU(N), &
          SP(IG)%ZLL_BU(N),SP(IG)%KLC_BU(N),SP(IG)%ILC_BU(N), &
          SP(IG)%JLC_BU(N),SP(IG)%TOT_BU(N),SP(IG)%MXL_BU(N), &
          SP(IG)%STM_BU(N))

 DO I=1,SP(IG)%NPART-1
  SP(IG)%XLC_BU(I)=SP(IG)%XLC(I)
  SP(IG)%YLC_BU(I)=SP(IG)%YLC(I)
  SP(IG)%ZLC_BU(I)=SP(IG)%ZLC(I)
  SP(IG)%ZLL_BU(I)=SP(IG)%ZLL(I)
  SP(IG)%KLC_BU(I)=SP(IG)%KLC(I)
  SP(IG)%ILC_BU(I)=SP(IG)%ILC(I)
  SP(IG)%JLC_BU(I)=SP(IG)%JLC(I)
  SP(IG)%TOT_BU(I)=SP(IG)%TOT(I)
  SP(IG)%MXL_BU(I)=SP(IG)%MXL(I)
  SP(IG)%STM_BU(I)=SP(IG)%STM(I)
 ENDDO

 DEALLOCATE(SP(IG)%XLC,SP(IG)%YLC,SP(IG)%ZLC, &
            SP(IG)%ZLL,SP(IG)%KLC,SP(IG)%ILC, &
            SP(IG)%JLC,SP(IG)%TOT,SP(IG)%MXL, &
            SP(IG)%STM)

 SP(IG)%XLC=>SP(IG)%XLC_BU
 SP(IG)%YLC=>SP(IG)%YLC_BU
 SP(IG)%ZLC=>SP(IG)%ZLC_BU
 SP(IG)%ZLL=>SP(IG)%ZLL_BU
 SP(IG)%KLC=>SP(IG)%KLC_BU
 SP(IG)%ILC=>SP(IG)%ILC_BU
 SP(IG)%JLC=>SP(IG)%JLC_BU
 SP(IG)%TOT=>SP(IG)%TOT_BU
 SP(IG)%MXL=>SP(IG)%MXL_BU
 SP(IG)%STM=>SP(IG)%STM_BU

 END SUBROUTINE TRACE_3D_STARTPOINTS_MEMORY_SP

 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_MEMORY_SPR(IG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IG
 INTEGER :: I,N
 
 IF(SPR(IG)%NPART.LE.SIZE(SPR(IG)%XLC))RETURN
 
 N=SPR(IG)%NPART+999
 
 ALLOCATE(SPR(IG)%XLC_BU(N),SPR(IG)%YLC_BU(N),SPR(IG)%ZLC_BU(N), &
          SPR(IG)%ZLL_BU(N),SPR(IG)%KLC_BU(N),SPR(IG)%ILC_BU(N), &
          SPR(IG)%JLC_BU(N),SPR(IG)%TOT_BU(N),SPR(IG)%MXL_BU(N), &
          SPR(IG)%STM_BU(N))

 !## reset layer
 SPR(IG)%KLC_BU=0
 
 DO I=1,SPR(IG)%NPART-1
  SPR(IG)%XLC_BU(I)=SPR(IG)%XLC(I)
  SPR(IG)%YLC_BU(I)=SPR(IG)%YLC(I)
  SPR(IG)%ZLC_BU(I)=SPR(IG)%ZLC(I)
  SPR(IG)%ZLL_BU(I)=SPR(IG)%ZLL(I)
  SPR(IG)%KLC_BU(I)=SPR(IG)%KLC(I)
  SPR(IG)%ILC_BU(I)=SPR(IG)%ILC(I)
  SPR(IG)%JLC_BU(I)=SPR(IG)%JLC(I)
  SPR(IG)%TOT_BU(I)=SPR(IG)%TOT(I)
  SPR(IG)%MXL_BU(I)=SPR(IG)%MXL(I)
  SPR(IG)%STM_BU(I)=SPR(IG)%STM(I)
 ENDDO

 DEALLOCATE(SPR(IG)%XLC,SPR(IG)%YLC,SPR(IG)%ZLC, &
            SPR(IG)%ZLL,SPR(IG)%KLC,SPR(IG)%ILC, &
            SPR(IG)%JLC,SPR(IG)%TOT,SPR(IG)%MXL, &
            SPR(IG)%STM)

 SPR(IG)%XLC=>SPR(IG)%XLC_BU
 SPR(IG)%YLC=>SPR(IG)%YLC_BU
 SPR(IG)%ZLC=>SPR(IG)%ZLC_BU
 SPR(IG)%ZLL=>SPR(IG)%ZLL_BU
 SPR(IG)%KLC=>SPR(IG)%KLC_BU
 SPR(IG)%ILC=>SPR(IG)%ILC_BU
 SPR(IG)%JLC=>SPR(IG)%JLC_BU
 SPR(IG)%TOT=>SPR(IG)%TOT_BU
 SPR(IG)%MXL=>SPR(IG)%MXL_BU
 SPR(IG)%STM=>SPR(IG)%STM_BU
    
 END SUBROUTINE TRACE_3D_STARTPOINTS_MEMORY_SPR
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_RESET()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I
 
 !## clear all
 DO I=1,SIZE(SP)
  !## open drawing list per timestep/call
  IF(STPLISTINDEX(I).NE.0)CALL GLDELETELISTS(STPLISTINDEX(I),1_GLSIZEI); STPLISTINDEX(I)=0
 ENDDO
        
 NSPG=0; SP%NPART=0; PL%NPART=0
 CALL TRACE_DEAL_SP(); CALL TRACE_DEAL_SPR() 

 !## maximal number of particles to be traced
 DO I=1,SIZE(SP); CALL TRACE_AL_SP(1000,I); ENDDO 

 CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER5,PL%NPART)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER11,NSPG)
 !## only if no particle simulation is running
 IF(PL%IRUN.EQ.0)CALL WDIALOGFIELDSTATE(ID_ADD,1)
 
 END SUBROUTINE TRACE_3D_STARTPOINTS_RESET
 
 !###======================================================================
 SUBROUTINE TRACE_3D_STARTPOINTS_SHOW()
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=GLDOUBLE) :: X,Y,Z
 INTEGER :: I,IG
  
! VIEWDX=(TOP%X-BOT%X)/2.0_GLDOUBLE/XYZAXES(1)
! VIEWDY=(TOP%Y-BOT%Y)/2.0_GLDOUBLE/XYZAXES(2)

 !## clear all
 DO IG=1,SIZE(SP)
  !## open drawing list per timestep/call
  IF(STPLISTINDEX(IG).NE.0)CALL GLDELETELISTS(STPLISTINDEX(IG),1_GLSIZEI); STPLISTINDEX(IG)=0
 ENDDO
  
 !## create seperate list for all of the particles
 DO IG=1,NSPG

  !## list index for,  !## start new drawing list
  STPLISTINDEX(IG)=GLGENLISTS(1); CALL GLNEWLIST(STPLISTINDEX(IG),GL_COMPILE)
 
  CALL GLBEGIN(GL_POINTS)

  DO I=1,SP(IG)%NPART
   !## current position
   X= SP(IG)%XLC(I)+IDF%XMIN; Y=SP(IG)%YLC(I)+IDF%YMIN; Z=SP(IG)%ZLC(I)
!   X=(X-MIDPOS%X)/VIEWDX; Y=(Y-MIDPOS%Y)/VIEWDY
   CALL GLVERTEX3D(X,Y,Z)
  ENDDO

  CALL GLEND()
  CALL GLENDLIST()

 ENDDO

 CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER5,PL%NPART)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER11,NSPG)

 END SUBROUTINE TRACE_3D_STARTPOINTS_SHOW
 
 !###======================================================================
 LOGICAL FUNCTION TRACE_3D_COMPUTE() 
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),PARAMETER :: DYEAR=365.25
 INTEGER :: I,IPART,IDSCH,NPACT,IG,ITRAPPED,NTAIL,IFREQ,ICAPT
 REAL(KIND=DP_KIND) :: TIME,TTMAX,MAXVELOCITY,XTREP
 
 TRACE_3D_COMPUTE=.FALSE.

 !## maximal time of simulation - equal to refreshing rate or something to be trickered by the 3d tool
 MAXVELOCITY=0.0D0

 CALL WDIALOGSELECT(ID_D3DSETTINGS_TAB8)
 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO3,PL%ITYPE)

 CALL WDIALOGGETDOUBLE(IDF_REAL4,PL%TMAX); PL%TMAX=PL%TMAX*DYEAR
 CALL WDIALOGGETDOUBLE(IDF_REAL5,PL%TDEL); PL%TDEL=PL%TDEL*DYEAR
 CALL WDIALOGGETDOUBLE(IDF_REAL6,XTREP)
 CALL WDIALOGGETDOUBLE(IDF_REAL7,PL%TCUR); PL%TCUR=PL%TCUR*DYEAR
 CALL WDIALOGGETINTEGER(IDF_INTEGER4,NTAIL)
 CALL WDIALOGGETCHECKBOX(IDF_CHECK7,ITRAPPED)
 CALL WDIALOGGETCHECKBOX(IDF_CHECK8,IFREQ)
 !## make sure ifreq=0 as itrapped=1
 IF(ITRAPPED.EQ.1)IFREQ=0
 CALL WDIALOGGETCHECKBOX(IDF_CHECK6,ICAPT)
 IF(ICAPT.EQ.1)CALL WDIALOGGETMENU(IDF_MENU4,ICAPT)

 CALL WDIALOGGETMENU(IDF_MENU3,ISNK); ISNK=ISNK-1
 CALL WDIALOGGETDOUBLE(IDF_REAL8,FRAC)
 
 PL%NTREP=INT((XTREP+1.0D0)*REAL(NTAIL))

 !## increase number of timesteps
 PL%NPER =PL%NPER+1
 PL%NTIME=PL%NTIME+1
 
 !## shift all one position
 DO I=1,SIZE(PLLISTCLR); PLLISTCLR(I)=PLLISTCLR(I)+1; ENDDO
 
 !## restart in stack
 IF(PL%NPER.GT.SIZE(PLLISTINDEX,1))PL%NPER=1

 !## get next timestep
 TTMAX=PL%TCUR+PL%TDEL
 !## maximize it for tmax
 TTMAX=MIN(TTMAX,PL%TMAX)

 !## current time/years
 PL%TCUR=TTMAX/DYEAR; CALL WDIALOGPUTDOUBLE(IDF_REAL7,PL%TCUR) 

 !## age of current drawing list
 PLLISTAGE(PL%NPER)=TTMAX
 !## position in colouring of current drawing list
 PLLISTCLR(PL%NPER)=1
 
 !## store original direction
 IREV=PL%IREV
 
 !## count number of active particles
 NPACT=0; IMODE=0; IMODE(1)=1
 
 DO IG=1,NSPG

  !## open drawing list per timestep/call
  IF(PLLISTINDEX(PL%NPER,IG).NE.0)CALL GLDELETELISTS(PLLISTINDEX(PL%NPER,IG),1_GLSIZEI)
  
  !## process only whenever particles are active
  IF(SP(IG)%IACT.NE.1)CYCLE
  
  !## decide direction for current particles
  IF(SP(IG)%IREV.NE.IREV)THEN
   CALL TRACEIREV(); IREV=ABS(IREV-1)
  ENDIF

  !## list index for, start new drawing list
  PLLISTINDEX(PL%NPER,IG)=GLGENLISTS(1); CALL GLNEWLIST(PLLISTINDEX(PL%NPER,IG),GL_COMPILE)
 
  !## points
  IF(PL%ITYPE.EQ.2)CALL GLBEGIN(GL_POINTS)

  DO IPART=1,SPR(IG)%NPART 

!   IF(IPART.EQ.647)THEN
!    WRITE(*,*)
!   ENDIF
   
   !## particle not yet started, see whether to start
   IF(SPR(IG)%STM(IPART).GT.0)THEN
    !## skip this one for now
    IF(SPR(IG)%STM(IPART).GT.PL%NTIME)CYCLE
   ENDIF
   
   !## trace selected particle, not yet discharged!
   IF(SPR(IG)%KLC(IPART).GT.0)THEN

    !## time in days
    TIME=SPR(IG)%TOT(IPART)

    !## compute ttmax per particle
    TTMAX=TIME+PL%TDEL
    !## maximize it for tmax
    TTMAX=MIN(TTMAX,PL%TMAX)   

    !## lines
    IF(PL%ITYPE.EQ.1)CALL GLBEGIN(GL_LINE_STRIP)
     
    NPACT=NPACT+1; SPR(IG)%MXL(IPART)=0
    CALL FLOLIN(IPART,IMODE(1),TIME,TTMAX,IDSCH, &
                SPR(IG)%JLC(IPART),SPR(IG)%ILC(IPART),SPR(IG)%KLC(IPART),   &
                SPR(IG)%XLC(IPART),SPR(IG)%YLC(IPART),SPR(IG)%ZLC(IPART),SPR(IG)%ZLL(IPART), &
                IBOUND,ZBOT,ZTOP,LDELR,LDELC,QX,QY,QZ,QSS,POR,NCON,IDF%NCOL,IDF%NROW,NLAY,  &
                NLPOR,IDF%NCOL*IDF%NROW*NLAY,NCP1,NRP1,NLP1,ISNK,IREV,FRAC,IMODE(1),   &
                MAXVELOCITY,DELX,DELY,SPR(IG)%MXL(IPART),IVISIT,LVISIT,NVISIT) 
    
!    IF(SPR(IG)%XLC(IPART).NE.SPR(IG)%XLC(IPART))THEN
!    WRITE(*,*)
!    ENDIF
    
    !## lines
    IF(PL%ITYPE.EQ.1)CALL GLEND()

    !## clean visited places that have been visited before
    DO I=1,NVISIT; IVISIT(LVISIT(I))=INT(0,1); ENDDO
    !## current time in days
    SPR(IG)%TOT(IPART)=TIME
    !## reached tmax - particle can continue next step
    IF(IDSCH.EQ.7)IDSCH=0
    IF(TIME.GE.PL%TMAX)IDSCH=7
    !## particle discharged whenever idsch.ne.0 or maximal time exceeds
    IF(IDSCH.NE.0)SPR(IG)%KLC(IPART)=-IDSCH 

   ELSE

    !## check whether to restart again
    IF(ITRAPPED.EQ.1)THEN
    
     !## particle stopped - decide whether to restart
     IF(ICAPT.GT.0.AND.SPR(IG)%KLC(IPART).LT.0)THEN
      !## no restart since the exitcode (idsch) is not correct
      IF(ICAPT.NE.ABS(SPR(IG)%KLC(IPART)))CYCLE
     ENDIF

     !## restart particle
     SPR(IG)%ILC(IPART)=SP(IG)%ILC(IPART)
     SPR(IG)%JLC(IPART)=SP(IG)%JLC(IPART)
     SPR(IG)%KLC(IPART)=SP(IG)%KLC(IPART)
     SPR(IG)%XLC(IPART)=SP(IG)%XLC(IPART)
     SPR(IG)%YLC(IPART)=SP(IG)%YLC(IPART)
     SPR(IG)%ZLC(IPART)=SP(IG)%ZLC(IPART)
     SPR(IG)%ZLL(IPART)=SP(IG)%ZLL(IPART)
     SPR(IG)%TOT(IPART)=SP(IG)%TOT(IPART)
     SPR(IG)%MXL(IPART)=SP(IG)%MXL(IPART)
     SPR(IG)%STM(IPART)=SP(IG)%STM(IPART)
     IF(SP(IG)%IRSTRT.GT.0)SPR(IG)%STM(IPART)=SPR(IG)%STM(IPART)+PL%NTIME

    ENDIF
   
   ENDIF
  ENDDO

  !## point strip
  IF(PL%ITYPE.EQ.2)CALL GLEND()
  
  CALL GLENDLIST()
 
  !## create new particles ttmax/pl%trep another integer ...
  IF(IFREQ.EQ.1)THEN
   !## initiate a new series of particles
   IF(MOD(PL%NTIME,PL%NTREP).EQ.0)THEN

    I=0
    !## duplicate
    DO IPART=1,SP(IG)%NPART

     !## particle stopped - decide whether to restart
     IF(ICAPT.GT.0.AND.SPR(IG)%KLC(IPART).LT.0)THEN
      !## still going
      IF(ICAPT.NE.ABS(SPR(IG)%KLC(IPART)))CYCLE
     ENDIF

     !## find "empty" spot of terminated particle
     DO
      I=I+1; IF(I.GT.SPR(IG)%NPART)EXIT
      IF(SPR(IG)%KLC(I).LE.0)EXIT
     ENDDO
     
     IF(I.GT.SPR(IG)%NPART)THEN
      !## count number of particles
      SPR(IG)%NPART=SPR(IG)%NPART+1
      CALL TRACE_3D_STARTPOINTS_MEMORY_SPR(IG)
     ENDIF

     !## add copy of particle 
     SPR(IG)%ILC(I)=SP(IG)%ILC(IPART)
     SPR(IG)%JLC(I)=SP(IG)%JLC(IPART)
     SPR(IG)%KLC(I)=SP(IG)%KLC(IPART)
     SPR(IG)%XLC(I)=SP(IG)%XLC(IPART)
     SPR(IG)%YLC(I)=SP(IG)%YLC(IPART)
     SPR(IG)%ZLC(I)=SP(IG)%ZLC(IPART)
     SPR(IG)%ZLL(I)=SP(IG)%ZLL(IPART)
     SPR(IG)%TOT(I)=SP(IG)%TOT(IPART)
     SPR(IG)%MXL(I)=SP(IG)%MXL(IPART)
     SPR(IG)%STM(I)=SP(IG)%STM(IPART)
     IF(SP(IG)%IRSTRT.GT.0)SPR(IG)%STM(IPART)=SPR(IG)%STM(IPART)+PL%NTIME
  
     !## add active particles
     NPACT=NPACT+1
    
    ENDDO
   ENDIF
   
  ENDIF
 ENDDO
  
 !## return direction, if neccessary
 IF(PL%IREV.NE.IREV)CALL TRACEIREV()
 
 !## clean drawing list, whenever nothing in it
 IF(NPACT.LE.0)THEN
  DO IG=1,NSPG
   CALL GLDELETELISTS(PLLISTINDEX(PL%NPER,IG),1_GLSIZEI); PLLISTINDEX(PL%NPER,IG)=0
  ENDDO
 ENDIF
 
 CALL WDIALOGPUTINTEGER(IDF_INTEGER6,NPACT)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER7,SUM(SPR%NPART))

 !## trace as long as there are active particles available
 IF(SUM(PLLISTINDEX).NE.0)TRACE_3D_COMPUTE=.TRUE.
 
 END FUNCTION TRACE_3D_COMPUTE

 !###======================================================================
 SUBROUTINE TRACE_3D_COMPUTE_STOP()
 !###======================================================================
 IMPLICIT NONE
 
 IF(PL%IRUN.EQ.0)RETURN
 
 CALL WDIALOGPUTSTRING(ID_START,'Start'); CALL WDIALOGFIELDSTATE(ID_STOP,0)
 CALL WDIALOGFIELDSTATE(IDF_RADIO1,1); CALL WDIALOGFIELDSTATE(IDF_RADIO2,1)
 CALL WDIALOGFIELDSTATE(IDF_INTEGER4,1)
 CALL TRACE_DEAL_SPR(); PL%IRUN=0

 CALL WDIALOGFIELDSTATE(ID_NEW,1); CALL WDIALOGFIELDSTATE(ID_ADD,1)

 CALL WDIALOGPUTINTEGER(IDF_INTEGER6,0)
 CALL WDIALOGPUTINTEGER(IDF_INTEGER7,0)
 CALL WDIALOGPUTDOUBLE(IDF_REAL7,0.0D0,'(F15.3)')

 END SUBROUTINE TRACE_3D_COMPUTE_STOP
 
 !###======================================================================
 SUBROUTINE TRACE_3D_CLOSE()
 !###======================================================================
 IMPLICIT NONE

 !## deallocate memory
 CALL TRACEDEALLOCATE(1)

! IF(ALLOCATED(IVISIT))DEALLOCATE(IVISIT)
! IF(ALLOCATED(LVISIT))DEALLOCATE(LVISIT)
 
 CALL WDIALOGSELECT(ID_D3DSETTINGS_PART);     CALL WDIALOGUNLOAD()
 CALL WDIALOGSELECT(ID_D3DSETTINGS_LAYER);    CALL WDIALOGUNLOAD()
 CALL WDIALOGSELECT(ID_D3DSETTINGS_LOCATION); CALL WDIALOGUNLOAD()
 CALL WDIALOGSELECT(ID_D3DSETTINGS_SINKS);    CALL WDIALOGUNLOAD()
 
 CALL WINDOWSELECT(0); CALL WMENUSETSTATE(ID_INTERACTIVEPATHLINES,2,0)
 CALL WMENUSETSTATE(ID_INTERACTIVEPATHLINES,1,1)
 
 END SUBROUTINE TRACE_3D_CLOSE
 
 !###======================================================================
 LOGICAL FUNCTION TRACEMAIN(RUNFILE,IBATCH,ICONVERTGEN,ZOOMIDF)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: ZOOMIDF
 CHARACTER(LEN=*),INTENT(IN) :: RUNFILE
 INTEGER,INTENT(IN) :: IBATCH,ICONVERTGEN
 INTEGER :: IULOG
 LOGICAL :: LEX
 
 TRACEMAIN=.FALSE.

 !## deallocate all memory
 CALL TRACEDEALLOCATE(1)

 !## read runfile
 IF(TRACEREADRUNFILE(RUNFILE,IBATCH))THEN

  IULOG=UTL_GETUNIT()
  CALL OSD_OPEN(IULOG,FILE=RUNFILE(:INDEX(RUNFILE,'.',.TRUE.)-1)//'.LOG',STATUS='UNKNOWN')

  !## trace particles
  IF(TRACE_CALC(IBATCH,IULOG,ICONVERTGEN,ZOOMIDF))TRACEMAIN=.TRUE.
  
  INQUIRE(UNIT=IULOG,OPENED=LEX); IF(LEX)CLOSE(IULOG)

 ENDIF
  
 END FUNCTION TRACEMAIN
 
 !###======================================================================
 LOGICAL FUNCTION TRACECALC_INIT(IBATCH,NCONS)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IBATCH
 INTEGER,INTENT(OUT) :: NCONS
 
 TRACECALC_INIT=.FALSE.

 !## allocate space
 CALL TRACEAL()
 !## create cell-sizes (cell-borders expressed in x,y-coordinates)
 CALL TRACEDELRC()
 !## read information-for particle tracking
 IF(.NOT.TRACEDATIN(NCONS,IBATCH))RETURN
 
 TRACECALC_INIT=.TRUE.
 
 END FUNCTION TRACECALC_INIT

 !###======================================================================
 SUBROUTINE TRACE_BRING_IN_DATA(MODELIDF,N,TOP,BOT)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: N
 TYPE(IDFOBJ),INTENT(INOUT) :: MODELIDF
 TYPE(IDFOBJ),INTENT(INOUT),DIMENSION(N) :: TOP,BOT
 INTEGER :: ILAY,IROW,ICOL,NROW,NCOL
 
 NROW=TOP(1)%NROW
 NCOL=TOP(1)%NCOL
 NLAY=N
 
 CALL IDFCOPY(MODELIDF,IDF)
 ALLOCATE(ZTOP(NCOL,NROW,NLAY),ZBOT(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY)); IBOUND=0
 DO ILAY=1,NLAY
  DO IROW=1,NROW; DO ICOL=1,NCOL
   ZTOP(ICOL,IROW,ILAY)=TOP(ILAY)%X(ICOL,IROW)
   ZBOT(ICOL,IROW,ILAY)=BOT(ILAY)%X(ICOL,IROW)
   IF(ZTOP(ICOL,IROW,ILAY).NE.TOP(ILAY)%NODATA.AND. &
      ZBOT(ICOL,IROW,ILAY).NE.BOT(ILAY)%NODATA)IBOUND(ICOL,IROW,ILAY)=1
  ENDDO; ENDDO
 ENDDO
 
 END SUBROUTINE TRACE_BRING_IN_DATA
 
 !###======================================================================
 LOGICAL FUNCTION TRACE_CALC(IBATCH,IULOG,ICONVERTGEN,ZOOMIDF)
 !###======================================================================
 IMPLICIT NONE
 TYPE(IDFOBJ),INTENT(IN) :: ZOOMIDF
 INTEGER,INTENT(IN) :: IBATCH,ICONVERTGEN,IULOG
 CHARACTER(LEN=50) :: CPERIOD
 CHARACTER(LEN=1000) :: STRING
 REAL(KIND=DP_KIND) :: TIME,TTMAX,DT,MAXVELOCITY
 TYPE(WIN_MESSAGE) :: MESSAGE
 INTEGER :: I,J,ILAY,ITYPE,IPER,IPART,NCONS,IPERIOD,NIDSCH,IDSCH,IP,IRAT, &
            IRAT1,DPER,MPER,SPER 
 LOGICAL :: LEX
 
 TRACE_CALC=.FALSE.

 !## take head for determination model-dimensions
 IF(.NOT.IDFREAD(IDF,HFFNAME(1,1,1),0))RETURN; CLOSE(IDF%IU)

 !## define zoom area for particle tracking - overrule area
 IF(ZOOMIDF%XMAX-ZOOMIDF%XMIN.GT.0.0D0)THEN
  IDF%XMIN=ZOOMIDF%XMIN; IDF%XMAX=ZOOMIDF%XMAX
  IDF%YMIN=ZOOMIDF%YMIN; IDF%YMAX=ZOOMIDF%YMAX
  CALL UTL_IDFSNAPTOGRID(IDF%XMIN,IDF%XMAX,IDF%YMIN,IDF%YMAX,IDF%DX,IDF%NCOL,IDF%NROW)
  WRITE(*,'(/1X,A/)') 'Particle tracking for a part of the model only'
 ENDIF

 IF(.NOT.TRACECALC_INIT(IBATCH,NCONS))RETURN
 
 ALLOCATE(IVISIT(IDF%NCOL*IDF%NROW*NLAY)); IVISIT=INT(0,1)
 ALLOCATE(LVISIT(IDF%NCOL*IDF%NROW*NLAY)); LVISIT=0

 !## initialize startpoints memory
 CALL TRACE_INIT_SP(NSPFNAME) 
 
 ALLOCATE(IUOUT(NSPFNAME,2)); IUOUT=0
 
 DO ISPFNAME=1,NSPFNAME
  
  !## read/process particles towards readable format
  IF(.NOT.TRACEPREPARESP(SP(ISPFNAME)%NPART,IBATCH,0))RETURN 
  !## initialize outputfiles
  IF(.NOT.TRACEINITOUTFILES(SP(ISPFNAME)%NPART))RETURN 
 
  !## read all particle in memory
  CALL TRACE_AL_SP(SP(ISPFNAME)%NPART,ISPFNAME) 
  !## set initial time for each particle to zero
  CALL TRACEREADSP(ISPFNAME,SP(ISPFNAME)%NPART)

 ENDDO
 
 !## copy particles to runtime particles
 CALL TRACE_INIT_SPR(NSPFNAME)
 CALL TRACE_AL_SPR()

 TTMAX  =0.0D0
 IPERIOD=1
 IRAT   =0
 IRAT1  =IRAT
 MAXVELOCITY=0.0D0
 IPATHLINE_3D=0
  
 !## transient simulation - get start timestep
 IF(ISS.EQ.1)THEN
  
  !## get start of window
  DO IPER=1,NPER
   IF(PLIPER(IPER,1).GT.JD1)EXIT
  ENDDO
  SPER=MAX(IPER-1,1)
  !## get end of window
  DO IPER=1,NPER
   IF(PLIPER(IPER,1).GT.JD2)EXIT
  ENDDO
  MPER=IPER-1
  
  !## forwards
  IF(IREV.EQ.0)THEN
   MPER=IPER-1
   DO IPER=1,NPER
    IF(PLIPER(IPER,1).LE.JD0.AND.PLIPER(IPER+1,1).GT.JD0)EXIT
   ENDDO
   IPER=IPER-1; DPER=1 !; SPER=1 !; MPER=NPER
  !## backwards
  ELSEIF(IREV.EQ.1)THEN
   DO IPER=NPER,1,-1
    IF(PLIPER(IPER,1).LE.JD0.AND.PLIPER(IPER+1,1).GE.JD0)EXIT
   ENDDO   
   IPER=IPER+1; DPER=-1; I=MPER; MPER=SPER; SPER=I
  ENDIF
   
 ENDIF
   
 DO

  !## transient simulation
  IF(ISS.EQ.1)THEN
   IPER=IPER+DPER   
   !## what to do after sequence has almost ended
   IF(IPER.EQ.MPER)THEN 
    LEX=.TRUE.; DT=PLIPER(IPER+1,1)-PLIPER(IPER,1)
    IF(ISTOPCRIT.EQ.3)THEN; TTMAX=TMAX; ELSE; TTMAX=TTMAX+DT; ENDIF
   ELSE
    !## within selected time-window
    LEX=.FALSE.
    !## get time difference between stress-periods
    IF(PLIPER(IPER+1,1).LE.JD2.AND.PLIPER(IPER,1).GE.JD1)THEN
     !## length of stress-period (days)
     DT =MIN(JD2,PLIPER(IPER+1,1))-MAX(JD1,PLIPER(IPER,1))
     TTMAX=TTMAX+DT
     IF(DT.GT.0.0D0)LEX=.TRUE.
    ENDIF
   ENDIF

  !## steady-state simulation
  ELSE
   TTMAX=TMAX
   LEX=.TRUE.
   IPER=1
  ENDIF

  !## never exceeds given tmax
  TTMAX=MIN(TTMAX,TMAX)

  IF(LEX)THEN

   CALL WMESSAGEPEEK(ITYPE,MESSAGE)
   IF(ITYPE.EQ.KEYDOWN.AND.MESSAGE%VALUE1.EQ.KEYESCAPE)THEN
    CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'Are you sure to terminate current pathline computation?'//CHAR(13)// &
     'Though, current results upto now, are already stored in an IFF-format','Question')
    IF(WINFODIALOG(4).EQ.1)EXIT
   ENDIF

   IF(ISS.EQ.0)WRITE(IULOG,'(1X,A,I10,A,F10.2,A)')    'Using following files for steady-state simulation'
   IF(ISS.EQ.1)WRITE(IULOG,'(1X,A,I10,A,2(F10.2,A))') 'Using following files for current stressperiod',IPER,' time upto: ', &
                             TTMAX,' days (timestep',DT,' days)'
   DO ILAY=1,NLAY
    STRING=TRIM(ITOS(ILAY))
    J=1; IF(ILAY.EQ.NLAY)J=0
    DO I=1,2+J+ISTO 
     IP=INDEX(HFFNAME(I,ILAY,IPER),'\',.TRUE.)+1
     STRING=TRIM(STRING)//','//TRIM(HFFNAME(I,ILAY,IPER)(IP:))
    END DO
    WRITE(IULOG,'(A)') TRIM(STRING)
   END DO
   WRITE(IULOG,*)

   !## read time-dependent data for particle tracking
   IF(.NOT.TRACEREADBUDGET(IPER,IBATCH))RETURN 
   !## backwards tracking
   IF(IREV.EQ.1)CALL TRACEIREV()

   !## start particle loop
   IF(IBATCH.EQ.0)CALL WINDOWSELECT(0)

   IF(ISS.EQ.1)THEN
    CPERIOD=' [Duration '//TRIM(ITOS(INT(DT)))//' day; Cycle '//TRIM(ITOS(IPERIOD))//', max. '//TRIM(ITOS(INT(TTMAX)))//' days]'
    I=0
    DO ISPFNAME=1,NSPFNAME
     DO IPART=1,SP(ISPFNAME)%NPART
      IF(SP(ISPFNAME)%KLC(IPART).NE.0)I=I+1
     ENDDO
    ENDDO
    STRING='Still tracing '//TRIM(ITOS(I))//' out of '//TRIM(ITOS(SUM(SP%NPART)))//' particles for Stress '// &
                             TRIM(ITOS(IPER))//' out of '//TRIM(ITOS(MPER))//TRIM(CPERIOD)
    IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
    IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
   ELSE
    STRING='Tracing '//TRIM(ITOS(SUM(SP%NPART)))//' particles (steady-state) ...'
    IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
    IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
   ENDIF

   !## loop over startpunt groups
   DO ISPFNAME=1,NSPFNAME

!    IF(IBATCH.EQ.1)WRITE(*,'(A)') 'Busy tracking particles ('//TRIM(ITOS(SP(ISPFNAME)%NPART))// &
!      ') startpoints from '//TRIM(SPFNAME(ISPFNAME))//' ...'

    NIDSCH=0
    DO IPART=1,SP(ISPFNAME)%NPART 
     CALL WMESSAGEPEEK(ITYPE,MESSAGE)
     !## time in days!
     TIME=SPR(ISPFNAME)%TOT(IPART) 
     !## trace selected particle, NOT YET discharged!
     IF(SPR(ISPFNAME)%KLC(IPART).NE.0)THEN
      SPR(ISPFNAME)%MXL(IPART)=0
      CALL FLOLIN(IPART,IMODE(1),TIME,TTMAX,IDSCH, &
                  SPR(ISPFNAME)%JLC(IPART),SPR(ISPFNAME)%ILC(IPART),SPR(ISPFNAME)%KLC(IPART),  &
                  SPR(ISPFNAME)%XLC(IPART),SPR(ISPFNAME)%YLC(IPART),SPR(ISPFNAME)%ZLC(IPART),  &
                  SPR(ISPFNAME)%ZLL(IPART), &
                  IBOUND,ZBOT,ZTOP,LDELR,LDELC,QX,QY,QZ,QSS,POR,NCON,IDF%NCOL,IDF%NROW,NLAY,  &
                  NLPOR,IDF%NCOL*IDF%NROW*NLAY,NCP1,NRP1,NLP1,ISNK,IREV,FRAC,IUOUT(ISPFNAME,1),   &
                  MAXVELOCITY,DELX,DELY,SPR(ISPFNAME)%MXL(IPART),IVISIT,LVISIT,NVISIT)
      !## clean visited places that have been visited before
      DO I=1,NVISIT; IVISIT(LVISIT(I))=INT(0,1); ENDDO

      !## time in days
      SPR(ISPFNAME)%TOT(IPART)=TIME !## days
      IF(ISS.EQ.1)THEN
       !## end of current simulation
       IF(IDSCH.EQ.7)IDSCH=0
      ELSE
       !## end of simulation reached!
       IF(SPR(ISPFNAME)%TOT(IPART).GE.TTMAX)IDSCH=7
      ENDIF
      !## particle discharged whenever idsch.ne.0
      IF(IDSCH.NE.0)THEN
       !## write endpoint information current particle that has stopped
       IF(IMODE(2).GT.0)CALL TRACECREATEIPF(IDSCH,IPART)
       SPR(ISPFNAME)%KLC(IPART)=0
      ELSE
       !## particle NOT discharged!
       NIDSCH=NIDSCH+1
      ENDIF
     ENDIF

     IF(IBATCH.EQ.0)CALL UTL_WAITMESSAGE(IRAT,IRAT1,IPART,SP(ISPFNAME)%NPART,'Busy tracking particles ('//TRIM(ITOS(SP(ISPFNAME)%NPART))// &
      ') startpoints from '//TRIM(SPFNAME(ISPFNAME))//' ...')

    ENDDO
   ENDDO

   !## no particles left, stop tracking process!
   IF(NIDSCH.EQ.0)EXIT
  ENDIF

  !## stop while doing a steady-state simulation
  IF(ISS.EQ.0)EXIT 
  !## stop whenever tmax.eq.ttmax
  IF(TTMAX.GE.TMAX)EXIT
  !## determing stopcriterion whenever transient simulation concerned!
  IF(IPER.EQ.MPER)THEN 
   !## stop after last period, assume last period duration is similar to previous one!
   !## continue until particle stops (given tmax, else tmax=10e30)
   IF(ISTOPCRIT.EQ.1.OR.ISTOPCRIT.EQ.3)EXIT
   !## repeat period again tmax
   IF(ISTOPCRIT.EQ.2)THEN
    IPER   =SPER+(-1.0D0*DPER)
    IPERIOD=IPERIOD+1
   ENDIF
  ENDIF

 ENDDO

 !## write remaining non-stopped particles for endpoints (IDSCH=0)
 IF(IMODE(2).GT.0)THEN
  IDSCH=0
  DO ISPFNAME=1,NSPFNAME
   DO IPART=1,SP(ISPFNAME)%NPART; IF(SPR(ISPFNAME)%KLC(IPART).NE.0)CALL TRACECREATEIPF(IDSCH,IPART); ENDDO
  ENDDO
 ENDIF

 IF(IMODE(1).GT.0)THEN
  DO ISPFNAME=1,NSPFNAME; CLOSE(IUOUT(ISPFNAME,1)); ENDDO; IMODE(1)=1
 ENDIF
 IF(IMODE(2).GT.0)THEN
  DO ISPFNAME=1,NSPFNAME; CLOSE(IUOUT(ISPFNAME,2)); ENDDO; IMODE(2)=1
 ENDIF
 
 DO ISPFNAME=1,NSPFNAME
  STRING='Completed particle tracking. Results are stored within:'
  IF(IMODE(1).GT.0)STRING=TRIM(STRING)//' '//TRIM(IFFFNAME(ISPFNAME))//'.IFF'
  IF(IMODE(2).GT.0)STRING=TRIM(STRING)//' '//TRIM(IFFFNAME(ISPFNAME))//'.IPF'
  STRING=TRIM(STRING)//' and added to the iMOD-manager. '// &
    TRIM(ITOS(SP(ISPFNAME)%NPART))//' particles were released out of '//TRIM(ITOS(TPART))//'. '//  &
   'Unreleased particles occured due to inactive/constant head boundary conditions '// &
   'and/or particles positioned above/beneath given thresshold. '// &
    TRIM(ITOS(NCONS))//' inconsequences were removed from top/bottom data! '//CHAR(13)//CHAR(13)// &
    'IMPORTANT: Maximum velocity that occured: '//TRIM(RTOS(MAXVELOCITY,'E',4))//' m/day'
  WRITE(IULOG,'(/A/)') TRIM(STRING)
  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)THEN
   CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,TRIM(STRING),'Information')
  ENDIF
 ENDDO

 !## deallocate memory
 CALL TRACEDEALLOCATE(0)

 DO ISPFNAME=1,NSPFNAME
  !## convert to gen if desired
  IF(ICONVERTGEN.EQ.1)CALL TRACECONVERTTOGEN(TRIM(IFFFNAME(ISPFNAME))//'.iff')
 ENDDO
 
 CLOSE(IULOG)

 TRACE_CALC=.TRUE.

 END FUNCTION TRACE_CALC

 !###======================================================================
 SUBROUTINE TRACECONVERTTOGEN(FNAME) 
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: FNAME
 CHARACTER(LEN=256) :: LINE !,FNAME
 INTEGER :: J,IOS,ID
 INTEGER,DIMENSION(3) :: IU
 TYPE IFFOBJ
  INTEGER :: ID,IL,IROW,ICOL
  REAL(KIND=DP_KIND) :: X,Y,Z,T,V
 END TYPE IFFOBJ
 TYPE(IFFOBJ),DIMENSION(2) :: IFF
 LOGICAL :: LEX
 
 !## nothing to do
 IF(IMODE(1).EQ.0)THEN; WRITE(*,'(A)') 'No export to gen since no iff is created'; RETURN; ENDIF
 
 WRITE(*,'(/A/)') 'Busy with converting to GEN format' 
 WRITE(*,'(A)') '- Converting file '//TRIM(FNAME)
 
 INQUIRE(FILE=FNAME,EXIST=LEX)
 IF(.NOT.LEX)THEN; WRITE(*,'(A)') 'Cannot find specified gen file'; RETURN; ENDIF

 IU(1)=UTL_GETUNIT(); CALL OSD_OPEN(IU(1),FILE=FNAME                                    ,STATUS='OLD'    ,ACTION='READ' ,IOSTAT=IOS,FORM='FORMATTED')
 IU(2)=UTL_GETUNIT(); CALL OSD_OPEN(IU(2),FILE=FNAME(:INDEX(FNAME,'.',.TRUE.)-1)//'.gen',STATUS='UNKNOWN',ACTION='WRITE',IOSTAT=IOS,FORM='FORMATTED')
 IU(3)=UTL_GETUNIT(); CALL OSD_OPEN(IU(3),FILE=FNAME(:INDEX(FNAME,'.',.TRUE.)-1)//'.dat',STATUS='UNKNOWN',ACTION='WRITE',IOSTAT=IOS,FORM='FORMATTED')

 !## skip header
 DO J=1,10; READ(IU(1),*); ENDDO
 READ(IU(1),*,IOSTAT=IOS) IFF(2)%ID,IFF(2)%IL,IFF(2)%X,IFF(2)%Y,IFF(2)%Z,IFF(2)%T,IFF(2)%V,IFF(2)%IROW,IFF(2)%ICOL
 ID=0
 DO
  READ(IU(1),*,IOSTAT=IOS) IFF(1)%ID,IFF(1)%IL,IFF(1)%X,IFF(1)%Y,IFF(1)%Z,IFF(1)%T,IFF(1)%V,IFF(1)%IROW,IFF(1)%ICOL
  IF(IOS.NE.0)EXIT
  IF(IFF(1)%ID.EQ.IFF(2)%ID)THEN
   ID=ID+1
   WRITE(IU(2),'(I10)') ID 
   WRITE(IU(2),'(2(F10.2,1X))') IFF(2)%X,IFF(2)%Y
   WRITE(IU(2),'(2(F10.2,1X))') IFF(1)%X,IFF(1)%Y
   WRITE(IU(2),'(A)') 'END'
   LINE=TRIM(ITOS(ID))//','//TRIM(ITOS(IFF(1)%IL))//','//TRIM(RTOS(IFF(1)%T,'E',5))
   WRITE(IU(3),'(A)') TRIM(LINE)
  ENDIF
  IFF(2)=IFF(1)
 ENDDO
 
 WRITE(IU(2),'(A)') 'END'

 DO J=1,SIZE(IU); CLOSE(IU(J)); ENDDO
 
 END SUBROUTINE TRACECONVERTTOGEN

 !###======================================================================
 SUBROUTINE TRACEPOSTPROCESSING(IFFFLOW,IPFFLOW,IDFFLOW,IPFFNAME,IPFICOL,ICONVERTGEN,TOPFNAME,BOTFNAME,IEXTRACT) 
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: IFFFLOW,IPFFLOW,IDFFLOW,IPFFNAME,TOPFNAME,BOTFNAME
 INTEGER,INTENT(IN),DIMENSION(4) :: IPFICOL
 INTEGER,INTENT(IN) :: ICONVERTGEN,IEXTRACT
 CHARACTER(LEN=256) :: LINE
 TYPE IFFOBJ
  INTEGER :: ID,ILAY,IROW,ICOL
  REAL(KIND=DP_KIND) :: X,Y,Z,T,V
 END TYPE IFFOBJ
 TYPE(IFFOBJ),DIMENSION(:),POINTER :: IFF,IFF_BU
 TYPE(IDFOBJ) :: TOPIDF,BOTIDF
 INTEGER :: IU,I,J,N,IOS,NROW,NCOL,EPIROW,EPICOL,EPILAY,EPCMTT,EPX,EPY,EPZ,J1,J2
 REAL(KIND=DP_KIND) :: EPS,Z1,Z2
 INTEGER,DIMENSION(:,:),ALLOCATABLE :: JU
 INTEGER,DIMENSION(:,:),ALLOCATABLE :: NP
 REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: XYZ
 CHARACTER(LEN=12) :: CATTRIB  
 LOGICAL :: LEX

 IMODE=0; IF(IFFFLOW.NE.'')IMODE(1)=1; IF(IPFFLOW.NE.'')IMODE(2)=1
 WRITE(*,'(/A/)') 'Busy with postprocessing particles ...'
 IF(IMODE(1).EQ.1)WRITE(*,'(1X,A)') 'Processing '//TRIM(IFFFLOW)
 IF(IMODE(2).EQ.1)WRITE(*,'(1X,A)') 'Processing '//TRIM(IPFFLOW)

 IF(TRIM(IPFFNAME).NE.'')THEN
  IF(.NOT.TRACEPOSTPROCESSING_INIT(IPFFNAME,IPFICOL))RETURN
  IF(TRIM(IDFFLOW).NE.'')THEN 
   IF(.NOT.IDFREAD(IDF,IDFFLOW,0))RETURN; EPS=IDF%DX/10.0D0
  ENDIF
  !## get ipf-information - equal to imodflow to assign well to locations
  DO I=1,SIZE(IPF); CALL IDFIROWICOL(IDF,IPF(I)%IROW,IPF(I)%ICOL,IPF(I)%X,IPF(I)%Y); ENDDO
 ELSE
  IF(TRIM(TOPFNAME).NE.'')THEN; IF(.NOT.IDFREAD(TOPIDF,TOPFNAME,1))RETURN; ENDIF
  IF(TRIM(BOTFNAME).NE.'')THEN; IF(.NOT.IDFREAD(BOTIDF,BOTFNAME,1))RETURN; ENDIF
  NUNQ=1; ALLOCATE(IPF(NUNQ)); IPF(1)%UNQLABEL='3D_iextract'//TRIM(ITOS(IEXTRACT)); IPF(1)%INQ=1
 ENDIF 

 !## open all files
 ALLOCATE(JU(NUNQ,2),NP(NUNQ,2)); JU=0
 DO I=1,NUNQ
  IF(IMODE(1).EQ.1)THEN
   JU(I,1)=UTL_GETUNIT(); J=INDEX(IFFFLOW,'.',.TRUE.)-1; IF(J.LE.0)J=LEN_TRIM(IFFFLOW)
   CALL OSD_OPEN(JU(I,1),FILE=IFFFLOW(:J)//'_'//TRIM(IPF(I)%UNQLABEL)//'.iff',STATUS='REPLACE', &
                         ACTION='WRITE',IOSTAT=IOS,FORM='FORMATTED')
   CALL TRACEINITOUTFILES_IFF(JU(I,1))
  ENDIF
  IF(IMODE(2).EQ.1)THEN
   JU(I,2)=UTL_GETUNIT(); J=INDEX(IPFFLOW,'.',.TRUE.)-1; IF(J.LE.0)J=LEN_TRIM(IPFFLOW)
   CALL OSD_OPEN(JU(I,2),FILE=IPFFLOW(:J)//'_'//TRIM(IPF(I)%UNQLABEL)//'_.ipf',STATUS='REPLACE', &
                         ACTION='WRITE',IOSTAT=IOS,FORM='FORMATTED')
   CALL TRACEINITOUTFILES_IPF(JU(I,2),0)
  ENDIF
 ENDDO
 
 NP=0; IU=UTL_GETUNIT()
 !## split pathlines to appropriate files
 IF(IMODE(1).EQ.1)THEN
  IMODE(1)=UTL_GETUNIT()
  CALL OSD_OPEN(IMODE(1),FILE=IFFFLOW,STATUS='OLD',ACTION='READ',IOSTAT=IOS,FORM='FORMATTED')
  READ(IMODE(1),*) NCOL; DO I=1,NCOL; READ(IMODE(1),*); ENDDO
  ALLOCATE(IFF(100))
  READ(IMODE(1),*,IOSTAT=IOS) IFF(1)%ID,IFF(1)%ILAY,IFF(1)%X,IFF(1)%Y,IFF(1)%Z,IFF(1)%T,IFF(1)%V,IFF(1)%IROW,IFF(1)%ICOL
  
  !## process entire particle
  N=1; DO
   N=N+1
   IF(N.GT.SIZE(IFF))THEN
    IF(N.GT.10000)THEN
     LINE='### Currently particle '//TRIM(ITOS(IFF(1)%ID))//' probably went wrong, will be skipped ###'
     WRITE(*,'(/A/)') TRIM(LINE)
     !## skip particle ... too big to process
     ALLOCATE(IFF(100))
     DO
      READ(IMODE(1),*) IFF(2)%ID,IFF(2)%ILAY,IFF(2)%X,IFF(2)%Y,IFF(2)%Z,IFF(2)%T,IFF(2)%V
      IF(IFF(2)%ID.NE.IFF(1)%ID)THEN; IFF(1)=IFF(2); N=2; EXIT; ENDIF
     ENDDO
    ELSE
     ALLOCATE(IFF_BU(N+100),STAT=IOS)
     DO I=1,SIZE(IFF); IFF_BU(I)=IFF(I); ENDDO; DEALLOCATE(IFF); IFF=>IFF_BU
    ENDIF
   ENDIF
   READ(IMODE(1),*,IOSTAT=IOS) IFF(N)%ID,IFF(N)%ILAY,IFF(N)%X,IFF(N)%Y,IFF(N)%Z,IFF(N)%T,IFF(N)%V,IFF(N)%IROW,IFF(N)%ICOL
   
   !## new particle starts - evaluate particle
   IF(IFF(1)%ID.NE.IFF(N)%ID.OR.IOS.NE.0)THEN

    !## see what ipf gets the data
IFFLOOP: DO I=1,SIZE(IPF)
    
     IF(TRIM(IPFFNAME).NE.'')THEN

      LEX=IFF(N-1)%ILAY.EQ.IPF(I)%ILAY.AND. &
          IFF(N-1)%IROW.EQ.IPF(I)%IROW.AND. &
          IFF(N-1)%ICOL.EQ.IPF(I)%ICOL
      J1=1; J2=N-1

     ELSE

      LEX=.FALSE.
      SELECT CASE (IEXTRACT)
       !## test if any of the fragments pass the 3d zone
       CASE (1,2,3)
        J1=0; J2=0; DO J=1,N-1
         IF(TRACEPOSTPROCESSING_Z1Z2(IFF(J)%X,IFF(J)%Y,TOPIDF,BOTIDF,Z1,Z2))THEN
          IF(IFF(J)%Z.LE.Z1.AND.IFF(J)%Z.GE.Z2)THEN
           !## start of the (first) zone
           IF(J1.EQ.0)J1=J
          ELSE
           !## end of the (first) zone
           IF(J1.NE.0.AND.J2.EQ.0)J2=J
          ENDIF
         ENDIF
        ENDDO
        !## okay somewhere crosses the 3d zone - determine what part to write
        IF(J1.NE.0)THEN
         LEX=.TRUE.
         IF(IEXTRACT.EQ.1)THEN
          J1=1; J2=N-1
         ELSEIF(IEXTRACT.EQ.2)THEN
          J2=J1; J1=1
         ELSEIF(IEXTRACT.EQ.3)THEN
          IF(J2.GT.0)THEN
           J1=J2; J2=N-1; IF(J1.EQ.0)J1=1
          ELSE
           LEX=.FALSE.
          ENDIF
         ENDIF
        ENDIF
       !## test end of pathline
       CASE (4)
        J=N-1
        IF(TRACEPOSTPROCESSING_Z1Z2(IFF(J)%X,IFF(J)%Y,TOPIDF,BOTIDF,Z1,Z2))THEN
         IF(IFF(J)%Z.LE.Z1.AND.IFF(J)%Z.GE.Z2)LEX=.TRUE.; J1=1; J2=J
        ENDIF
      END SELECT

     ENDIF
      
     IF(LEX)THEN
      NP(IPF(I)%INQ,1)=NP(IPF(I)%INQ,1)+1
      DO J=J1,J2
       WRITE(JU(IPF(I)%INQ,1),'(2(I10,1X),5(E15.7,1X),2(I10,1X))') IFF(J)%ID,IFF(J)%ILAY,IFF(J)%X,IFF(J)%Y,IFF(J)%Z,IFF(J)%T,IFF(J)%V,IFF(J)%IROW,IFF(J)%ICOL
      ENDDO
      EXIT IFFLOOP
     ENDIF   
     
    ENDDO IFFLOOP
    
    IF(IOS.NE.0)EXIT
    IFF(1)=IFF(N); N=1
 
   ENDIF
  
  ENDDO

  DO I=1,NUNQ; IF(JU(I,1).NE.0)THEN
   IF(NP(I,1).EQ.0)THEN; CLOSE(JU(I,1),STATUS='DELETE'); JU(I,1)=0; ENDIF
   IF(NP(I,1).GT.0)CLOSE(JU(I,1))
  ENDIF; ENDDO
  CLOSE(IMODE(1)); IMODE(1)=1; DEALLOCATE(IFF)
  
 ENDIF

 !## open output channel endpoints
 IF(IMODE(2).EQ.1)THEN
  IMODE(2)=UTL_GETUNIT()
  CALL OSD_OPEN(IMODE(2),FILE=IPFFLOW,STATUS='OLD',ACTION='READ',IOSTAT=IOS,FORM='FORMATTED')
  READ(IMODE(2),*) NROW; READ(IMODE(2),*) NCOL
  !## get appropriate column due to version compatability
  DO I=1,NCOL;
   READ(IMODE(2),'(A)') CATTRIB
   IF(TRIM(CATTRIB).EQ.'EP_IROW')    EPIROW=I
   IF(TRIM(CATTRIB).EQ.'END_IROW')   EPIROW=I
   IF(TRIM(CATTRIB).EQ.'EP_ICOL')    EPICOL=I
   IF(TRIM(CATTRIB).EQ.'END_ICOL')   EPICOL=I
   IF(TRIM(CATTRIB).EQ.'EP_ILAY')    EPILAY=I
   IF(TRIM(CATTRIB).EQ.'END_ILAY')   EPILAY=I
   IF(TRIM(CATTRIB).EQ.'TIME(YEARS)')EPCMTT=I
   IF(TRIM(CATTRIB).EQ.'EP_XCRD.')   EPX   =I
   IF(TRIM(CATTRIB).EQ.'END_XCRD.')  EPX   =I
   IF(TRIM(CATTRIB).EQ.'EP_YCRD.')   EPY   =I
   IF(TRIM(CATTRIB).EQ.'END_YCRD.')  EPY   =I
   IF(TRIM(CATTRIB).EQ.'EP_ZCRD.')   EPZ   =I
   IF(TRIM(CATTRIB).EQ.'END_ZCRD.')  EPZ   =I
  ENDDO
  READ(IMODE(2),*)
!  WRITE(*,'(A,4I3)') 'EPICOL,EPIROW,EPILAY,EPCMTT ',EPICOL,EPIROW,EPILAY,EPCMTT
 
  ALLOCATE(XYZ(NCOL))
  !## split endpoint to appropriate files
  DO J=1,NROW

   READ(IMODE(2),'(A256)') LINE
   READ(LINE,*) (XYZ(I),I=1,NCOL) 

   IF(TRIM(IPFFNAME).NE.'')THEN

    !## see what ipf gets the data
IPFLOOP: DO I=1,SIZE(IPF)
     IF(INT(XYZ(EPILAY)).EQ.IPF(I)%ILAY)THEN
      IF(IPF(I)%IROW.EQ.INT(XYZ(EPIROW)).AND.IPF(I)%ICOL.EQ.INT(XYZ(EPICOL)))THEN
       NP(IPF(I)%INQ,2)=NP(IPF(I)%INQ,2)+1
       WRITE(JU(IPF(I)%INQ,2),'(A)') TRIM(LINE)
       EXIT IPFLOOP
      ENDIF
     ENDIF
    ENDDO IPFLOOP
   
   ELSE

    IF(TRACEPOSTPROCESSING_Z1Z2(XYZ(EPX),XYZ(EPY),TOPIDF,BOTIDF,Z1,Z2))THEN
     IF(XYZ(EPZ).LE.Z1.AND.XYZ(EPZ).GE.Z2)THEN  
      NP(IPF(1)%INQ,2)=NP(IPF(1)%INQ,2)+1
      WRITE(JU(IPF(1)%INQ,2),'(A)') TRIM(LINE)
     ENDIF
    ENDIF
     
   ENDIF
   
  ENDDO

  DEALLOCATE(XYZ)
  DO I=1,NUNQ; IF(JU(I,2).NE.0)THEN
   IF(NP(I,2).EQ.0)THEN; CLOSE(JU(I,2),STATUS='DELETE'); JU(I,2)=0; ENDIF
   IF(NP(I,2).GT.0)CLOSE(JU(I,2))
  ENDIF; ENDDO
  CLOSE(IMODE(2)); IMODE(2)=1
 ENDIF

 WRITE(*,'(/3A10)') 'UNIQUE','N_IFF','N_IPF'
 !## construct header for ipf files
 DO I=1,NUNQ
  WRITE(*,'(3I10)') I,NP(I,1),NP(I,2)
  IF(IMODE(2).EQ.1)THEN
   IF(NP(I,2).GT.0)THEN
    JU(I,2)=UTL_GETUNIT()
    J=INDEX(IPFFLOW,'.',.TRUE.)-1; IF(J.LE.0)J=LEN_TRIM(IPFFLOW)
    CALL OSD_OPEN(JU(I,2),FILE=IPFFLOW(:J)//'_'//TRIM(IPF(I)%UNQLABEL)//'_.ipf',STATUS='OLD', &
                          ACTION='READ',IOSTAT=IOS,FORM='FORMATTED')
    JU(I,1)=UTL_GETUNIT()
    CALL OSD_OPEN(JU(I,1),FILE=IPFFLOW(:J)//'_'//TRIM(IPF(I)%UNQLABEL)//'.ipf',STATUS='REPLACE', &
                          ACTION='WRITE',IOSTAT=IOS,FORM='FORMATTED')  
    READ(JU(I,2),*); WRITE(JU(I,1),*) NP(I,2)
    DO
     READ(JU(I,2),'(A256)',IOSTAT=IOS) LINE
     IF(IOS.NE.0)EXIT; WRITE(JU(I,1),'(A)') TRIM(LINE)
    ENDDO
    CLOSE(JU(I,2),STATUS='DELETE'); CLOSE(JU(I,1))

   ENDIF
  ENDIF
 ENDDO
 
 !## convert to gen?
 IF(ICONVERTGEN.EQ.1)THEN
  IF(IMODE(1).EQ.1)THEN
   DO I=1,NUNQ
    CALL TRACECONVERTTOGEN(IFFFLOW(:J)//'_'//TRIM(IPF(I)%UNQLABEL)//'.iff')
   ENDDO
  ENDIF
 ENDIF

 END SUBROUTINE TRACEPOSTPROCESSING
 
 !###======================================================================
 LOGICAL FUNCTION TRACEPOSTPROCESSING_Z1Z2(X,Y,TOPIDF,BOTIDF,Z1,Z2)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(IN) :: X,Y
 REAL(KIND=DP_KIND),INTENT(OUT) :: Z1,Z2
 TYPE(IDFOBJ),INTENT(IN) :: TOPIDF,BOTIDF
 INTEGER :: IROW,ICOL
 
 TRACEPOSTPROCESSING_Z1Z2=.FALSE.
 
 Z1=HUGE(1.0D0)
 IF(ASSOCIATED(TOPIDF%X))THEN
  CALL IDFIROWICOL(TOPIDF,IROW,ICOL,X,Y)
  IF(ICOL.GT.0.AND.ICOL.LE.TOPIDF%NCOL.AND. &
     IROW.GT.0.AND.IROW.LT.TOPIDF%NROW)THEN
   Z1=TOPIDF%X(ICOL,IROW); IF(Z1.EQ.TOPIDF%NODATA)RETURN
  ELSE
   !## outside idf extent
   RETURN
  ENDIF
 ENDIF

 Z2=-HUGE(1.0D0)
 IF(ASSOCIATED(BOTIDF%X))THEN
  CALL IDFIROWICOL(BOTIDF,IROW,ICOL,X,Y)
  IF(ICOL.GT.0.AND.ICOL.LE.BOTIDF%NCOL.AND. &
     IROW.GT.0.AND.IROW.LT.BOTIDF%NROW)THEN
   Z2=BOTIDF%X(ICOL,IROW); IF(Z2.EQ.BOTIDF%NODATA)RETURN
  ELSE
   !## outside idf extent
   RETURN
  ENDIF
 ENDIF

 TRACEPOSTPROCESSING_Z1Z2=.TRUE.
 
 END FUNCTION TRACEPOSTPROCESSING_Z1Z2

 !###======================================================================
 LOGICAL FUNCTION TRACEPOSTPROCESSING_INIT(IPFFNAME,IPFICOL)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: IPFFNAME
 INTEGER,INTENT(IN),DIMENSION(4) :: IPFICOL
 CHARACTER(LEN=52),ALLOCATABLE,DIMENSION(:) :: STRING
 INTEGER :: I,J,K,M,IOS,IU,NJ
 
 TRACEPOSTPROCESSING_INIT=.FALSE.

 IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=IPFFNAME,STATUS='OLD',FORM='FORMATTED',IOSTAT=IOS)
 READ(IU,*) M
 ALLOCATE(IPF(M))
 READ(IU,*) NJ; DO J=1,NJ+1; READ(IU,*); ENDDO
 ALLOCATE(STRING(NJ))
 DO J=1,M
  READ(IU,*) (STRING(K),K=1,NJ)
  READ(STRING(IPFICOL(1)),*) IPF(J)%X;     READ(STRING(IPFICOL(2)),*) IPF(J)%Y
  READ(STRING(IPFICOL(3)),*) IPF(J)%LABEL; READ(STRING(IPFICOL(4)),*) IPF(J)%ILAY
  IPF(J)%LABEL=UTL_CAP(IPF(J)%LABEL,'U')
 ENDDO
 DEALLOCATE(STRING); CLOSE(IU)

 !## get unique items
 IPF%UNQLABEL=IPF%LABEL
 CALL UTL_GETUNIQUE_CHAR(IPF%UNQLABEL,SIZE(IPF),NUNQ)
 IPF%INQ=0
 !## search for unique labels
 DO J=1,SIZE(IPF)
  DO I=1,NUNQ
   IF(IPF(J)%LABEL.EQ.IPF(I)%UNQLABEL)THEN; IPF(J)%INQ=I; EXIT; ENDIF
  ENDDO
 ENDDO
  
 !## found following unique labels
 WRITE(*,'(/A)') 'Found following unique labels'
 WRITE(*,'(A10,A52/)') 'Number','Label'
 DO I=1,NUNQ
  K=0; DO J=1,SIZE(IPF)
   IF(IPF(J)%INQ.EQ.I)K=K+1
  ENDDO
  WRITE(*,'(I10,A52)') K,TRIM(IPF(I)%UNQLABEL)
 ENDDO
 
 TRACEPOSTPROCESSING_INIT=.TRUE.
 
 END FUNCTION TRACEPOSTPROCESSING_INIT
 
 !###======================================================================
 LOGICAL FUNCTION TRACEINITOUTFILES(NPART)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NPART
 INTEGER :: IOS,I,J
 CHARACTER(LEN=50) :: ERRORMESSAGE

 TRACEINITOUTFILES=.FALSE.

 CALL UTL_CREATEDIR(IFFFNAME(ISPFNAME)(:INDEX(IFFFNAME(ISPFNAME),'\',.TRUE.)-1))

 I=INDEX(IFFFNAME(ISPFNAME),'.',.TRUE.)-1
 J=INDEX(IFFFNAME(ISPFNAME),'\',.TRUE.)
 !## no filename given ...
 IF(I.LT.J)I=LEN_TRIM(IFFFNAME(ISPFNAME))
 
 !## open output channel pathlines
 IF(IMODE(1).GT.0)THEN
  IUOUT(ISPFNAME,1)=UTL_GETUNIT()
  CALL OSD_OPEN(IUOUT(ISPFNAME,1),FILE=IFFFNAME(ISPFNAME)(:I)//'.iff',STATUS='REPLACE',ACTION='WRITE',IOSTAT=IOS,FORM='FORMATTED')
  CALL TRACEINITOUTFILES_IFF(IUOUT(ISPFNAME,1))
 ENDIF
 IF(IMODE(2).GT.0)THEN
  IUOUT(ISPFNAME,2)=UTL_GETUNIT()
  CALL OSD_OPEN(IUOUT(ISPFNAME,2),FILE=IFFFNAME(ISPFNAME)(:I)//'.ipf',STATUS='REPLACE',ACTION='WRITE',IOSTAT=IOS)
  CALL TRACEINITOUTFILES_IPF(IUOUT(ISPFNAME,2),NPART)
 ENDIF
 IF(IOS.NE.0)THEN
  CALL OSD_IOSTAT_MSG(IOS,ERRORMESSAGE)
  IF(IMODE(1).GT.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK, &
                     TRIM(ERRORMESSAGE)//CHAR(13)//TRIM(IFFFNAME(ISPFNAME)),'Error')
  IF(IMODE(2).GT.1)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK, &
                     TRIM(ERRORMESSAGE)//CHAR(13)//TRIM(IFFFNAME(ISPFNAME)),'Error')
  RETURN
 ENDIF

 TRACEINITOUTFILES=.TRUE.

 END FUNCTION TRACEINITOUTFILES

 !###======================================================================
 SUBROUTINE TRACEINITOUTFILES_IFF(IU)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IU

 WRITE(IU,'(A)') '9'
 WRITE(IU,'(A)') 'PARTICLE_NUMBER'
 WRITE(IU,'(A)') 'ILAY'
 WRITE(IU,'(A)') 'XCRD.'
 WRITE(IU,'(A)') 'YCRD.'
 WRITE(IU,'(A)') 'ZCRD.'
 WRITE(IU,'(A)') 'TIME(YEARS)'
 WRITE(IU,'(A)') 'VELOCITY(M/DAY)'
 WRITE(IU,'(A)') 'IROW'
 WRITE(IU,'(A)') 'ICOL'

 END SUBROUTINE TRACEINITOUTFILES_IFF

!###======================================================================
 SUBROUTINE TRACEINITOUTFILES_IPF(IU,NPART)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IU,NPART

 WRITE(IU,'(I10)') NPART
 WRITE(IU,'(A)') '17'
 WRITE(IU,'(A)') 'SP_XCRD.'
 WRITE(IU,'(A)') 'SP_YCRD.'
 WRITE(IU,'(A)') 'SP_ZCRD.'
 WRITE(IU,'(A)') 'SP_ILAY'
 WRITE(IU,'(A)') 'SP_IROW'
 WRITE(IU,'(A)') 'SP_ICOL'
 WRITE(IU,'(A)') 'EP_XCRD.'
 WRITE(IU,'(A)') 'EP_YCRD.'
 WRITE(IU,'(A)') 'EP_ZCRD.'
 WRITE(IU,'(A)') 'EP_ILAY'
 WRITE(IU,'(A)') 'EP_IROW'
 WRITE(IU,'(A)') 'EP_ICOL'
 WRITE(IU,'(A)') 'TIME(YEARS)'
 WRITE(IU,'(A)') 'MAXLAYER'
 WRITE(IU,'(A)') 'DISTANCE'
 WRITE(IU,'(A)') 'IDENT.NO.'
 WRITE(IU,'(A)') 'CAPTURED_BY'
 WRITE(IU,'(A)') '0,TXT'

 END SUBROUTINE TRACEINITOUTFILES_IPF

 !###======================================================================
 LOGICAL FUNCTION TRACEREADRUNFILE(RUNFILE,IBATCH)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IBATCH
 CHARACTER(LEN=*),INTENT(IN) :: RUNFILE
 INTEGER :: IU,J,IDATE,IPER,ILAY,NJ,IOS,I
 LOGICAL :: LEX
 REAL(KIND=DP_KIND) :: X
 CHARACTER(LEN=256) :: LINE
 
 TRACEREADRUNFILE=.FALSE.

 IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=RUNFILE,STATUS='OLD',FORM='FORMATTED',IOSTAT=IOS)
 IF(.NOT.TRACECHECKRUN(IOS,'Cannot open runfile:'//TRIM(RUNFILE),IU))RETURN
 !## optional output types
 READ(IU,*,IOSTAT=IOS) NLAY; IF(.NOT.TRACECHECKRUN(IOS,'NLAY',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I4)') 'NLAY: ',NLAY
 READ(IU,'(A256)',IOSTAT=IOS) LINE; IF(.NOT.TRACECHECKRUN(IOS,'NPER',IU))RETURN
 READ(LINE,*,IOSTAT=IOS) NPER,ISTO
 IF(IOS.NE.0)THEN
  ISTO=0; READ(LINE,*,IOSTAT=IOS) NPER; IF(.NOT.TRACECHECKRUN(IOS,'NPER',IU))RETURN
  IF(IBATCH.EQ.1)WRITE(*,'(A10,I4)') 'NPER: ',NPER
 ELSE
  IF(IBATCH.EQ.1)WRITE(*,'(A10,2I4)') 'NPER,ISTO: ',NPER,ISTO
 ENDIF

 IF(ALLOCATED(HFFNAME)) DEALLOCATE(HFFNAME)
 IF(ALLOCATED(ITBFNAME))DEALLOCATE(ITBFNAME)
 IF(ALLOCATED(SPFNAME)) DEALLOCATE(SPFNAME)
 IF(ALLOCATED(IFFFNAME))DEALLOCATE(IFFFNAME)

 IF(ISTO.EQ.0)THEN
  ALLOCATE(HFFNAME(3,NLAY,NPER))  !## bdgfrf,bdgfff,bdgflf
 ELSE
  ALLOCATE(HFFNAME(4,NLAY,NPER))  !## bdgfrf,bdgfff,bdgflf,bdgsto
 ENDIF
 ALLOCATE(ITBFNAME(5,NLAY))       !## ibound,top,bot,por_aqf,por_aqt
 HFFNAME=''; ITBFNAME=''
 
 NSPFNAME=1; ALLOCATE(SPFNAME(NSPFNAME),IFFFNAME(NSPFNAME))
 
 READ(IU,'(A256)',IOSTAT=IOS) SPFNAME(1); IF(.NOT.TRACECHECKRUN(IOS,'SPFNAME',IU))RETURN
 !## try to read sequence number
 READ(SPFNAME(1),*,IOSTAT=IOS) NSPFNAME 
 IF(IOS.EQ.0)THEN

  DEALLOCATE(SPFNAME,IFFFNAME)
  ALLOCATE(SPFNAME(NSPFNAME),IFFFNAME(NSPFNAME))
  DO I=1,NSPFNAME
   READ(IU,*,IOSTAT=IOS) SPFNAME(I)
   SPFNAME(I)=UTL_SUBST(SPFNAME(I),TRIM(REPLACESTRING),PREFVAL(5))
   IF(.NOT.TRACECHECKRUN(IOS,'SPFNAME('//TRIM(ITOS(I))//')',IU))RETURN
   LINE='SPFNAME: '//TRIM(ITOS(I))
   IF(IBATCH.EQ.1)WRITE(*,'(A10,A)') TRIM(LINE),TRIM(SPFNAME(I))
   READ(IU,*,IOSTAT=IOS) IFFFNAME(I)
   IFFFNAME(I)=UTL_SUBST(IFFFNAME(I),TRIM(REPLACESTRING),PREFVAL(5))
   IF(.NOT.TRACECHECKRUN(IOS,'IFFFNAME('//TRIM(ITOS(I))//')',IU))RETURN
   LINE='IFFNAME: '//TRIM(ITOS(I))
   IF(IBATCH.EQ.1)WRITE(*,'(A10,A)') TRIM(LINE),TRIM(IFFFNAME(I))
  ENDDO

 ELSE

  NSPFNAME=1; READ(SPFNAME(1),*,IOSTAT=IOS) SPFNAME(1)
  SPFNAME(1)=UTL_SUBST(SPFNAME(1),TRIM(REPLACESTRING),PREFVAL(5))
  IF(IBATCH.EQ.1)WRITE(*,'(A)') 'SPFNAME: '//TRIM(SPFNAME(1))
  READ(IU,*,IOSTAT=IOS) IFFFNAME(1); IF(.NOT.TRACECHECKRUN(IOS,'IFFFNAME',IU))RETURN
  IFFFNAME(1)=UTL_SUBST(IFFFNAME(1),TRIM(REPLACESTRING),PREFVAL(5))
  IF(IBATCH.EQ.1)WRITE(*,'(A)') 'IFFNAME: '//TRIM(IFFFNAME(1))

 ENDIF

 READ(IU,*,IOSTAT=IOS) (IMODE(I),I=1,2); IF(.NOT.TRACECHECKRUN(IOS,'IMODE(.)',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,2I2)') 'IMODE(:): ',IMODE
 READ(IU,*,IOSTAT=IOS) IREV ; IF(.NOT.TRACECHECKRUN(IOS,'IREV',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I2)') 'IREV: ',IREV
 READ(IU,*,IOSTAT=IOS) ISNK;  IF(.NOT.TRACECHECKRUN(IOS,'ISNK',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I2)') 'ISNK: ',ISNK
 ISNK=ISNK-1
 READ(IU,*,IOSTAT=IOS) FRAC;  IF(.NOT.TRACECHECKRUN(IOS,'FRAC',IU))RETURN
 IF(ISNK.NE.2)FRAC=1.0D0; FRAC=MIN(1.0D0,MAX(0.0D0,FRAC))
 IF(IBATCH.EQ.1)WRITE(*,'(A10,F5.2)') 'FRAC: ',FRAC
 READ(IU,*,IOSTAT=IOS) ISTOPCRIT; IF(.NOT.TRACECHECKRUN(IOS,'ISTOPCRIT',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I2)') 'ISTOPCRIT: ',ISTOPCRIT
 READ(IU,*,IOSTAT=IOS) TMAX;  IF(.NOT.TRACECHECKRUN(IOS,'TMAX',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,E10.3)') 'TMAX: ',TMAX
 JD0=0; JD1=0; JD2=0
 READ(IU,*,IOSTAT=IOS) IDATE; IF(.NOT.TRACECHECKRUN(IOS,'START DATE',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I10)') 'START DATE: ',IDATE
 IF(NPER.GT.1)JD0=UTL_IDATETOJDATE(IDATE)
 READ(IU,*,IOSTAT=IOS) IDATE; IF(.NOT.TRACECHECKRUN(IOS,'START WINDOW',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I10)') 'START WINDOW: ',IDATE
 IF(NPER.GT.1)JD1=UTL_IDATETOJDATE(IDATE)
 READ(IU,*,IOSTAT=IOS) IDATE; IF(.NOT.TRACECHECKRUN(IOS,'END WINDOW',IU))RETURN
 IF(IBATCH.EQ.1)WRITE(*,'(A10,I10)') 'END WINDOW: ',IDATE
 IF(NPER.GT.1)JD2=UTL_IDATETOJDATE(IDATE)

 !## ib,top,bot,por_aqf,por_aqt
 DO ILAY=1,NLAY
  DO J=1,5
   IF(J.NE.5.OR.ILAY.NE.NLAY)THEN
    READ(IU,*,IOSTAT=IOS) ITBFNAME(J,ILAY)

    ITBFNAME(J,ILAY)=UTL_CAP(ITBFNAME(J,ILAY),'U')
    ITBFNAME(J,ILAY)=UTL_SUBST(ITBFNAME(J,ILAY),TRIM(REPLACESTRING),PREFVAL(5))

    IF(J.EQ.1)THEN

     IF(.NOT.TRACECHECKRUN(IOS,'IBOUND FOR LAYER '//TRIM(ITOS(ILAY)),IU))RETURN
     LINE='IBOUND'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY))
     IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
     IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'IBOUND'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY)))

    ELSEIF(J.EQ.2)THEN

     IF(.NOT.TRACECHECKRUN(IOS,'TOP FOR LAYER '//TRIM(ITOS(ILAY)),IU))RETURN
     LINE='TOP'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY))
     IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
     IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'TOP'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY)))

    ELSEIF(J.EQ.3)THEN

     IF(.NOT.TRACECHECKRUN(IOS,'BOT FOR LAYER '//TRIM(ITOS(ILAY)),IU))RETURN
     LINE='BOT'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY))
     IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
     IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'BOT'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY)))

    ELSEIF(J.EQ.4)THEN

     IF(.NOT.TRACECHECKRUN(IOS,'POR AQUIFER FOR LAYER '//TRIM(ITOS(ILAY)),IU))RETURN
     LINE='POR AQUIFER'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY))
     IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
     IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'POR AQUIFER'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY)))

    ELSEIF(J.EQ.5)THEN

     IF(.NOT.TRACECHECKRUN(IOS,'POR AQUITARD FOR LAYER '//TRIM(ITOS(ILAY)),IU))RETURN
     LINE='POR AQUITARD'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY))
     IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
     IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'POR AQUITARD'//TRIM(ITOS(ILAY))//'='//TRIM(ITBFNAME(J,ILAY)))

    ENDIF

    READ(ITBFNAME(J,ILAY),*,IOSTAT=IOS) X
    IF(IOS.NE.0)THEN
     INQUIRE(FILE=ITBFNAME(J,ILAY),EXIST=LEX)
     IF(.NOT.LEX)THEN
      IF(.NOT.TRACECHECKRUN(1,CHAR(13)//'File not found: '//TRIM(ITBFNAME(J,ILAY)),IU))RETURN      
      RETURN
     ENDIF
     IOS=0
    ENDIF

   ENDIF
  END DO
 END DO
 !## number of input budget files (including bdgsto if isto.eq.1)
 NJ=3+ISTO
 !## steady-state (iss.eq.0) or transient/multiple stressperiods (iss.eq.1)
 ISS=0; IF(NPER.GT.1)ISS=1

 !## frf,fff,flf,(sto)
 DO IPER=1,NPER
  DO ILAY=1,NLAY
   DO J=1,NJ
    IF(J.NE.3.OR.ILAY.NE.NLAY)THEN
     READ(IU,*,IOSTAT=IOS) HFFNAME(J,ILAY,IPER)

     HFFNAME(J,ILAY,IPER)=UTL_CAP(HFFNAME(J,ILAY,IPER),'U')
     HFFNAME(J,ILAY,IPER)=UTL_SUBST(HFFNAME(J,ILAY,IPER),TRIM(REPLACESTRING),PREFVAL(5))

     IF(J.EQ.1)THEN

      IF(.NOT.TRACECHECKRUN(IOS,'BDGFRF FOR LAYER '//TRIM(ITOS(ILAY))//' IPER '//TRIM(ITOS(IPER)),IU))RETURN
      LINE='BDGFRF FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER))
      IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
      IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'BDGFRF FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER)))

     ELSEIF(J.EQ.2)THEN

      IF(.NOT.TRACECHECKRUN(IOS,'BDGFFF FOR LAYER '//TRIM(ITOS(ILAY))//' IPER '//TRIM(ITOS(IPER)),IU))RETURN
      LINE='BDGFFF FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER))
      IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(LINE)
      IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'BDGFFF FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER)))

     ELSEIF(J.EQ.3)THEN

      IF(.NOT.TRACECHECKRUN(IOS,'BDGFLF FOR LAYER '//TRIM(ITOS(ILAY))//' IPER '//TRIM(ITOS(IPER)),IU))RETURN
      LINE='BDGFLF FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER))
      IF(IBATCH.EQ.1)WRITE(*,'(A)') LINE
      IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'BDGFLF FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER)))

     ELSEIF(J.EQ.4)THEN

      IF(.NOT.TRACECHECKRUN(IOS,'BDGSTO FOR LAYER '//TRIM(ITOS(ILAY))//' IPER '//TRIM(ITOS(IPER)),IU))RETURN
      LINE='BDGSTO FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER))
      IF(IBATCH.EQ.1)WRITE(*,'(A)') LINE
      IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'BDGSTO FOR LAYER'//TRIM(ITOS(ILAY))//'='//TRIM(HFFNAME(J,ILAY,IPER)))

     ENDIF

     INQUIRE(FILE=HFFNAME(J,ILAY,IPER),EXIST=LEX)
     IF(.NOT.LEX)THEN
      IF(.NOT.TRACECHECKRUN(1,'File not found: '//TRIM(HFFNAME(J,ILAY,IPER)),IU))RETURN      
      RETURN
     ENDIF

    ENDIF
   END DO
  END DO
 END DO
 CLOSE(IU)

 !## get time data
 IF(ALLOCATED(PLIPER))DEALLOCATE(PLIPER)
 ALLOCATE(PLIPER(NPER+1,1)); PLIPER=0 
 DO IPER=1,NPER
  PLIPER(IPER,1)=UTL_IDATETOJDATE(UTL_IDFGETDATE(HFFNAME(1,1,IPER)))
 ENDDO
 CALL WSORT(PLIPER(:,1),1,NPER)
 !## artifically extent with one day
 PLIPER(NPER+1,1)=PLIPER(NPER,1)+(JD2-PLIPER(NPER,1))
 
 TRACEREADRUNFILE=.TRUE.

 END FUNCTION TRACEREADRUNFILE

 !###======================================================================
 LOGICAL FUNCTION TRACECHECKRUN(IOS,TXT,IU)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IOS,IU
 CHARACTER(LEN=*),INTENT(IN) :: TXT

 TRACECHECKRUN=.TRUE.

 IF(IOS.EQ.0)RETURN
 CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Error reading parameter in runfile.'//CHAR(13)//TRIM(TXT),'Error')
 CLOSE(IU)
 
 TRACECHECKRUN=.FALSE.

 END FUNCTION TRACECHECKRUN

 !###======================================================================
 SUBROUTINE TRACEDEALLOCATE(ICODE)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ICODE

 CALL TRACE_DEAL_SP();  IF(ALLOCATED(SP)) DEALLOCATE(SP)
 CALL TRACE_DEAL_SPR(); IF(ALLOCATED(SPR))DEALLOCATE(SPR)
  
 IF(ICODE.EQ.0)RETURN

 IF(ALLOCATED(ITBFNAME))DEALLOCATE(ITBFNAME)
 IF(ALLOCATED(HFFNAME)) DEALLOCATE(HFFNAME)
 IF(ALLOCATED(QX))      DEALLOCATE(QX);     IF(ALLOCATED(QY))   DEALLOCATE(QY)
 IF(ALLOCATED(QZ))      DEALLOCATE(QZ);     IF(ALLOCATED(QS))   DEALLOCATE(QS)
 IF(ALLOCATED(POR))     DEALLOCATE(POR)
 IF(ALLOCATED(DELR))    DEALLOCATE(DELR);   IF(ALLOCATED(DELC)) DEALLOCATE(DELC)
 IF(ALLOCATED(DELX))    DEALLOCATE(DELX);   IF(ALLOCATED(DELY)) DEALLOCATE(DELY)
 IF(ALLOCATED(LDELR))   DEALLOCATE(LDELR);  IF(ALLOCATED(LDELC))DEALLOCATE(LDELC)
 IF(ALLOCATED(ZBOT))    DEALLOCATE(ZBOT);   IF(ALLOCATED(ZTOP)) DEALLOCATE(ZTOP)
 IF(ALLOCATED(IBOUND))  DEALLOCATE(IBOUND); IF(ALLOCATED(NCON)) DEALLOCATE(NCON)
 IF(ALLOCATED(QSS))  DEALLOCATE(QSS)
 IF(ALLOCATED(SPFNAME)) DEALLOCATE(SPFNAME)
 IF(ALLOCATED(IFFFNAME))DEALLOCATE(IFFFNAME)
 IF(ALLOCATED(IVISIT))  DEALLOCATE(IVISIT)
 IF(ALLOCATED(LVISIT))  DEALLOCATE(LVISIT)
 IF(ALLOCATED(IUOUT))   DEALLOCATE(IUOUT)
 END SUBROUTINE TRACEDEALLOCATE

 !###======================================================================
 SUBROUTINE TRACECREATEIPF(IDSCH,IPART)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IDSCH,IPART
 INTEGER  :: IU,IG
 REAL(KIND=DP_KIND) :: DIST
 
 IU=IUOUT(ISPFNAME,2)
 IG=ISPFNAME
  
 DIST=SQRT((SP(IG)%XLC(IPART)-SPR(IG)%XLC(IPART))**2.0D0+ &
           (SP(IG)%YLC(IPART)-SPR(IG)%YLC(IPART))**2.0D0+ &
           (SP(IG)%ZLC(IPART)-SPR(IG)%ZLC(IPART))**2.0D0)

 WRITE(IU,'(2(3(F15.3,1X),3(I10,1X)),E15.7,1X,I10,1X,F15.3,1X,2(I10,1X))') &
       SP(IG) %XLC(IPART)+IDF%XMIN, &
       SP(IG) %YLC(IPART)+IDF%YMIN, &
       SP(IG) %ZLC(IPART),          &
       SP(IG) %KLC(IPART),          &
       SP(IG) %ILC(IPART),          &
       SP(IG) %JLC(IPART),          &
       SPR(IG)%XLC(IPART)+IDF%XMIN, &
       SPR(IG)%YLC(IPART)+IDF%YMIN, &
       SPR(IG)%ZLC(IPART),          &
       SPR(IG)%KLC(IPART),          &
       SPR(IG)%ILC(IPART),          &
       SPR(IG)%JLC(IPART),          &
       SPR(IG)%TOT(IPART)/365.25D0, &
       SPR(IG)%MXL(IPART),          &
       DIST,                        &
       IPART,                       &
       IDSCH

 END SUBROUTINE TRACECREATEIPF

 !###======================================================================
 REAL(KIND=DP_KIND) FUNCTION TRACEGETV(XYZT,APOR)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(IN) :: APOR
 REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(4,3) :: XYZT
 REAL(KIND=DP_KIND) :: DX,DY,DZ,DT,D

 DX=XYZT(1,2)-XYZT(1,1)
 DY=XYZT(2,2)-XYZT(2,1)
 DZ=XYZT(3,2)-XYZT(3,1)
 DT=XYZT(4,2)-XYZT(4,1)
 D =SQRT(DX**2.0D0+DY**2.0D0+DZ**2.0D0)
 !#m/year
 TRACEGETV=D/DT*APOR

 RETURN
 END FUNCTION

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

 NCP1 =IDF%NCOL+1
 NRP1 =IDF%NROW+1
 NLP1 =NLAY+1
 !## porosity also for aquitards
 NLPOR=NLAY+NLAY-1

 ALLOCATE(QX(NCP1,IDF%NROW,NLAY))
 ALLOCATE(QY(IDF%NCOL,NRP1,NLAY))
 ALLOCATE(QZ(IDF%NCOL,IDF%NROW,NLP1))
 IF(ISTO.EQ.1)ALLOCATE(QS(IDF%NCOL,IDF%NROW,NLAY))
 ALLOCATE(POR(IDF%NCOL,IDF%NROW,NLPOR))
 ALLOCATE(DELR(0:IDF%NCOL))
 ALLOCATE(DELC(0:IDF%NROW))
 ALLOCATE(LDELR(IDF%NCOL))
 ALLOCATE(LDELC(IDF%NROW))
 ALLOCATE(DELX(IDF%NCOL))
 ALLOCATE(DELY(IDF%NROW))
 ALLOCATE(ZBOT(IDF%NCOL,IDF%NROW,NLAY))
 ALLOCATE(ZTOP(IDF%NCOL,IDF%NROW,NLAY))
 ALLOCATE(QSS(IDF%NCOL,IDF%NROW,NLAY))
 ALLOCATE(IBOUND(IDF%NCOL,IDF%NROW,NLAY))
 ALLOCATE(NCON(NLAY))

 RETURN
 END SUBROUTINE

 !###======================================================================
 SUBROUTINE TRACE_INIT_SP(NG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NG
 INTEGER :: I
 
 !## particle groups
 ALLOCATE(SP(NG)); SP%NPART=0; SP%ICLR=0; SP%IACT=0; SP%IREV=0; SP%SPWIDTH=0.0D0; SP%PWIDTH=0.0D0
 
 DO I=1,NG

  NULLIFY(SP(I)%XLC   ,SP(I)%YLC   ,SP(I)%ZLC   ,SP(I)%ZLL   , &
          SP(I)%ILC   ,SP(I)%JLC   ,SP(I)%KLC   ,SP(I)%TOT   , &
          SP(I)%MXL   ,SP(I)%STM)
  NULLIFY(SP(I)%XLC_BU,SP(I)%YLC_BU,SP(I)%ZLC_BU,SP(I)%ZLL_BU, &
          SP(I)%ILC_BU,SP(I)%JLC_BU,SP(I)%KLC_BU,SP(I)%TOT_BU, &
          SP(I)%MXL_BU,SP(I)%STM_BU)

 ENDDO
 
 END SUBROUTINE TRACE_INIT_SP

 !###======================================================================
 SUBROUTINE TRACE_INIT_SPR(NG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NG
 INTEGER :: I
 
 !## particle groups
 ALLOCATE(SPR(NG)); SPR%NPART=0; SPR%ICLR=0; SPR%IACT=0; SPR%IREV=0

 DO I=1,NG

  NULLIFY(SPR(I)%XLC   ,SPR(I)%YLC   ,SPR(I)%ZLC   ,SPR(I)%ZLL   , &
          SPR(I)%ILC   ,SPR(I)%JLC   ,SPR(I)%KLC   ,SPR(I)%TOT   , &
          SPR(I)%MXL   ,SPR(I)%STM)
  NULLIFY(SPR(I)%XLC_BU,SPR(I)%YLC_BU,SPR(I)%ZLC_BU,SPR(I)%ZLL_BU, &
          SPR(I)%ILC_BU,SPR(I)%JLC_BU,SPR(I)%KLC_BU,SPR(I)%TOT_BU, &
          SPR(I)%MXL_BU,SPR(I)%STM_BU)

 ENDDO
 
 END SUBROUTINE TRACE_INIT_SPR
 
 !###======================================================================
 SUBROUTINE TRACE_AL_SP(NPART,IG)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NPART,IG
 
 !## reals
 ALLOCATE(SP(IG)%XLC(NPART), &
          SP(IG)%YLC(NPART), &
          SP(IG)%ZLC(NPART), &
          SP(IG)%ZLL(NPART), &
          SP(IG)%TOT(NPART), &
          SP(IG)%STM(NPART))

 !## integers
 ALLOCATE(SP(IG)%KLC(NPART), &
          SP(IG)%JLC(NPART), &
          SP(IG)%ILC(NPART), &
          SP(IG)%MXL(NPART))
 
 END SUBROUTINE TRACE_AL_SP
 
 !###======================================================================
 SUBROUTINE TRACE_AL_SPR()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J
 
 CALL TRACE_DEAL_SPR(); IF(ALLOCATED(SPR))DEALLOCATE(SPR)
 
 ALLOCATE(SPR(SIZE(SP)))
 
 !## reset particles
 DO I=1,SIZE(SPR)

  !## reals
  ALLOCATE(SPR(I)%XLC(SP(I)%NPART), &
           SPR(I)%YLC(SP(I)%NPART), &
           SPR(I)%ZLC(SP(I)%NPART), &
           SPR(I)%ZLL(SP(I)%NPART), &
           SPR(I)%TOT(SP(I)%NPART), &
           SPR(I)%STM(SP(I)%NPART))
  !## integers
  ALLOCATE(SPR(I)%KLC(SP(I)%NPART), &
           SPR(I)%JLC(SP(I)%NPART), &
           SPR(I)%ILC(SP(I)%NPART), &
           SPR(I)%MXL(SP(I)%NPART))
  
  SPR(I)%NPART=SP(I)%NPART
  
  DO J=1,SP(I)%NPART
   SPR(I)%XLC(J)=SP(I)%XLC(J)
   SPR(I)%YLC(J)=SP(I)%YLC(J)
   SPR(I)%ZLC(J)=SP(I)%ZLC(J)
   SPR(I)%ZLL(J)=SP(I)%ZLL(J)
   SPR(I)%KLC(J)=SP(I)%KLC(J)
   SPR(I)%ILC(J)=SP(I)%ILC(J)
   SPR(I)%JLC(J)=SP(I)%JLC(J)
   SPR(I)%TOT(J)=SP(I)%TOT(J)
   SPR(I)%MXL(J)=SP(I)%MXL(J)
   SPR(I)%STM(J)=SP(I)%STM(J)
  ENDDO

 ENDDO
     
 END SUBROUTINE TRACE_AL_SPR
 
 !###======================================================================
 SUBROUTINE TRACE_DEAL_SP()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I

 IF(.NOT.ALLOCATED(SP))RETURN
 
 !## reset particles
 DO I=1,SIZE(SP)
 
  IF(ASSOCIATED(SP(I)%XLC))DEALLOCATE(SP(I)%XLC)
  IF(ASSOCIATED(SP(I)%YLC))DEALLOCATE(SP(I)%YLC)
  IF(ASSOCIATED(SP(I)%ZLC))DEALLOCATE(SP(I)%ZLC)
  IF(ASSOCIATED(SP(I)%ILC))DEALLOCATE(SP(I)%ILC)
  IF(ASSOCIATED(SP(I)%JLC))DEALLOCATE(SP(I)%JLC)
  IF(ASSOCIATED(SP(I)%KLC))DEALLOCATE(SP(I)%KLC)
  IF(ASSOCIATED(SP(I)%ZLL))DEALLOCATE(SP(I)%ZLL)
  IF(ASSOCIATED(SP(I)%TOT))DEALLOCATE(SP(I)%TOT)
  IF(ASSOCIATED(SP(I)%MXL))DEALLOCATE(SP(I)%MXL)
  IF(ASSOCIATED(SP(I)%STM))DEALLOCATE(SP(I)%STM)

 ENDDO

 END SUBROUTINE TRACE_DEAL_SP
 
 !###======================================================================
 SUBROUTINE TRACE_DEAL_SPR()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I
 
 IF(.NOT.ALLOCATED(SPR))RETURN

 !## reset particles
 DO I=1,SIZE(SPR)

  IF(ASSOCIATED(SPR(I)%XLC))DEALLOCATE(SPR(I)%XLC)
  IF(ASSOCIATED(SPR(I)%YLC))DEALLOCATE(SPR(I)%YLC)
  IF(ASSOCIATED(SPR(I)%ZLC))DEALLOCATE(SPR(I)%ZLC)
  IF(ASSOCIATED(SPR(I)%ILC))DEALLOCATE(SPR(I)%ILC)
  IF(ASSOCIATED(SPR(I)%JLC))DEALLOCATE(SPR(I)%JLC)
  IF(ASSOCIATED(SPR(I)%KLC))DEALLOCATE(SPR(I)%KLC)
  IF(ASSOCIATED(SPR(I)%ZLL))DEALLOCATE(SPR(I)%ZLL)
  IF(ASSOCIATED(SPR(I)%TOT))DEALLOCATE(SPR(I)%TOT)
  IF(ASSOCIATED(SPR(I)%MXL))DEALLOCATE(SPR(I)%MXL)
  IF(ASSOCIATED(SPR(I)%STM))DEALLOCATE(SPR(I)%STM)
  
 ENDDO
 
 END SUBROUTINE TRACE_DEAL_SPR
 
 !###======================================================================
 LOGICAL FUNCTION TRACEREADBUDGET(IPER,IBATCH)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),PARAMETER :: QCRIT=0.01D0  !## smaller value might be noise
 INTEGER,INTENT(IN) :: IBATCH,IPER
 CHARACTER(LEN=256) :: STRING
 INTEGER :: ILAY,IROW,ICOL,ITYPE
 REAL(KIND=DP_KIND) :: QERROR
 TYPE(WIN_MESSAGE) :: MESSAGE
 TYPE(IDFOBJ) :: IDFTMP
 
 TRACEREADBUDGET=.FALSE.

 CALL IDFCOPY(IDF,IDFTMP)

 IDFTMP%XMIN=IDF%XMIN-IDF%DX 
 IDFTMP%YMAX=IDF%YMAX
 IDFTMP%NCOL=IDF%NCOL+1
 IDFTMP%NROW=IDF%NROW

 !## read flux-right-face
 QX  =0.0D0
 DO ILAY=1,NLAY
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING='Reading BDGFRF L'//TRIM(ITOS(ILAY))//' '//TRIM(HFFNAME(1,ILAY,IPER))//'...'
!  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)

  !## set this to an equidistantial grid otherwise sx/sy will not be recomputed
  IDFTMP%IEQ=0

  IF(.NOT.IDFREADSCALE(HFFNAME(1,ILAY,IPER),IDFTMP,10,1,0.0D0,0))RETURN !## block-value
  DO ICOL=1,IDF%NCOL+1; DO IROW=1,IDF%NROW
   IF(IDFTMP%X(ICOL,IROW).EQ.IDFTMP%NODATA)THEN
    QX(ICOL,IROW,ILAY)=0.0D0
   ELSE
    QX(ICOL,IROW,ILAY)=-IDFTMP%X(ICOL,IROW)  !## moved allready a single cell to the left
   ENDIF
  ENDDO; ENDDO
 ENDDO

 IDFTMP%XMIN=IDF%XMIN
 IDFTMP%YMAX=IDF%YMAX+IDF%DX 
 IDFTMP%NCOL=IDF%NCOL
 IDFTMP%NROW=IDF%NROW+1

 !## read flux-front-face
 QY  =0.0D0
 DO ILAY=1,NLAY
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING='Reading BDGFFF L'//TRIM(ITOS(ILAY))//' '//TRIM(HFFNAME(2,ILAY,IPER))//'...'
!  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)

  !## set this to an equidistantial grid otherwise sx/sy will not be recomputed
  IDFTMP%IEQ=0

  IF(.NOT.IDFREADSCALE(HFFNAME(2,ILAY,IPER),IDFTMP,10,1,0.0D0,0))RETURN !## block-value
  DO ICOL=1,IDF%NCOL; DO IROW=1,IDF%NROW+1
   IF(IDFTMP%X(ICOL,IROW).EQ.IDFTMP%NODATA)THEN
    QY(ICOL,IROW,ILAY)=0.0D0
   ELSE
    QY(ICOL,IROW,ILAY)= IDFTMP%X(ICOL,IROW) 
   ENDIF
  ENDDO; ENDDO
 ENDDO

 !## read flux-lower-face
 QZ  =0.0D0
 IF(NLAY.GT.1)THEN
  DO ILAY=1,NLAY-1
   CALL WMESSAGEPEEK(ITYPE,MESSAGE)
   STRING='Reading BDGFLF L'//TRIM(ITOS(ILAY))//' '//TRIM(HFFNAME(3,ILAY,IPER))//'...'
!   IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
   IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)

   !## set this to an equidistantial grid otherwise sx/sy will not be recomputed
   IDFTMP%IEQ=0

   IF(.NOT.IDFREADSCALE(HFFNAME(3,ILAY,IPER),IDF,10,1,0.0D0,0))RETURN !## block-value
   DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
    IF(IDF%X(ICOL,IROW).EQ.IDF%NODATA)THEN
     QZ(ICOL,IROW,ILAY)=0.0D0
    ELSE
     QZ(ICOL,IROW,ILAY+1)=IDF%X(ICOL,IROW) 
    ENDIF
   ENDDO; ENDDO
  ENDDO
 ENDIF

 !## read storage
 IF(ISTO.EQ.1)THEN
  QS=0.0D0
  DO ILAY=1,NLAY
   CALL WMESSAGEPEEK(ITYPE,MESSAGE)
   STRING='Reading BDGSTO L'//TRIM(ITOS(ILAY))//' '//TRIM(HFFNAME(4,ILAY,IPER))//'...'
!   IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
   IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
   IF(.NOT.IDFREADSCALE(HFFNAME(4,ILAY,IPER),IDF,10,1,0.0D0,0))RETURN !## block-value
   DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
    IF(IDF%X(ICOL,IROW).EQ.IDF%NODATA)THEN
     QS(ICOL,IROW,ILAY)=0.0D0
    ELSE
     QS(ICOL,IROW,ILAY)=IDF%X(ICOL,IROW) 
    ENDIF
   ENDDO; ENDDO
  ENDDO
 ENDIF

 !## initialize qss - as nett-term waterbalance
 QSS=0.0D0
 DO ILAY=1,NLAY
  DO IROW=1,IDF%NROW
   DO ICOL=1,IDF%NCOL
    
    !## reset ibound 
    IF(IBOUND(ICOL,IROW,ILAY).EQ. 2000)IBOUND(ICOL,IROW,ILAY)=1
    IF(IBOUND(ICOL,IROW,ILAY).EQ.-2000)IBOUND(ICOL,IROW,ILAY)=1
    
    QERROR=QX(ICOL,IROW,ILAY)-QX(ICOL+1,IROW,ILAY)- &
           QY(ICOL,IROW,ILAY)+QY(ICOL,IROW+1,ILAY)- &
           QZ(ICOL,IROW,ILAY)+QZ(ICOL,IROW,ILAY+1)

    !## add storage as a known source
    IF(ISTO.EQ.1)QERROR=QERROR+QS(ICOL,IROW,ILAY)

    QSS(ICOL,IROW,ILAY)=-QERROR
    !## original code flags all buff().lt.0 - potential escape location
    IF(QSS(ICOL,IROW,ILAY).LT.-QCRIT.AND.IBOUND(ICOL,IROW,ILAY).NE.0)IBOUND(ICOL,IROW,ILAY)=-2000
    IF(QSS(ICOL,IROW,ILAY).GT. QCRIT.AND.IBOUND(ICOL,IROW,ILAY).NE.0)IBOUND(ICOL,IROW,ILAY)= 2000
   ENDDO
  ENDDO
 ENDDO
 
! idf%fname='d:\qss.idf'
! idf%x=qss(:,:,1)
! if(idfwrite(idf,idf%fname,1))then; endif

 CALL IDFDEALLOCATEX(IDFTMP)

 TRACEREADBUDGET=.TRUE.

 END FUNCTION TRACEREADBUDGET

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

 QX =-1.0D0*QX
 QY =-1.0D0*QY
 QZ =-1.0D0*QZ

 END SUBROUTINE TRACEIREV

 !###======================================================================
 LOGICAL FUNCTION TRACEDATIN(NCONS,IBATCH)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),PARAMETER :: MINTHICKNESS=0.01D0 !## OTHERWISE TRACING IS PROBLEM, DIVIDING BY ZERO AND NO PARTICLE BEHAVIOR IN Z-DIRECTION
 INTEGER,INTENT(IN) :: IBATCH
 INTEGER,INTENT(OUT) :: NCONS
 CHARACTER(LEN=256) :: STRING
 REAL(KIND=DP_KIND) :: APOR,XZTOP,XZBOT
 INTEGER :: ILAY,ICOL,IROW,IOS,IBND,ITYPE
 LOGICAL :: LEX
 TYPE(WIN_MESSAGE) :: MESSAGE
 
 TRACEDATIN=.FALSE.

 !## all modellayers do have confining beds by default
 NCON=1

 IF(IBATCH.EQ.0)CALL WINDOWSELECT(0)

 !## read iboundary idf-files
 DO ILAY=1,NLAY
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING=ITBFNAME(1,ILAY)
  READ(STRING,*,IOSTAT=IOS) IBND
  IF(IOS.EQ.0)THEN
   IBOUND(1:IDF%NCOL,1:IDF%NROW,ILAY)=IBND
   STRING='Constant value ['//TRIM(ITOS(IBND))//'] for IBOUND L'//TRIM(ITOS(ILAY))
  ELSE
   STRING='Reading IBOUND L'//TRIM(ITOS(ILAY))//' '//TRIM(ITBFNAME(1,ILAY))//'...'
   LEX=.TRUE.
   IF(ILAY.GT.1)THEN
    LEX=.TRUE.
    IF(TRIM(ITBFNAME(1,ILAY)).EQ.TRIM(ITBFNAME(1,ILAY-1)))LEX=.FALSE.
   ENDIF
   IF(LEX)THEN
    IF(.NOT.IDFREADSCALE(ITBFNAME(1,ILAY),IDF,1,1,0.0D0,0))RETURN !## scale boundary
    IBOUND(:,:,ILAY)=INT(IDF%X)
   ELSE
    STRING='Reusing '//TRIM(ITBFNAME(1,ILAY-1))
    IBOUND(:,:,ILAY)=IBOUND(:,:,ILAY-1)
   ENDIF
  ENDIF
  IF(IBATCH.EQ.1)WRITE(*,*) TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)

  DO IROW=1,IDF%NROW
   DO ICOL=1,IDF%NCOL

    IF(IBOUND(ICOL,IROW,ILAY).LT.0)THEN
     IBOUND(ICOL,IROW,ILAY)=-1
     !## POTENTIAL DISCHARGE CELLS ARE FLAGGED IN THE IBOUND ARRAY
     !## BY MULTIPLYING IBOUND VALUE BY 1000. THIS LOOP FLAGS CHANGING
     !## CELLS THAT HAVE NET DISCHARGE.
     IBOUND(ICOL,IROW,ILAY)=1000*IBOUND(ICOL,IROW,ILAY)
    ENDIF
    IF(IBOUND(ICOL,IROW,ILAY).GT.0)IBOUND(ICOL,IROW,ILAY)= 1
   END DO
  END DO
  !## FLAG SURROUNDING MODEL AS DISCHARGE AREA!!!!

 ENDDO

 !## read vertical coordinate data (TOP)
 DO ILAY=1,NLAY
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING=ITBFNAME(2,ILAY)
  READ(STRING,*,IOSTAT=IOS) XZTOP
  IF(IOS.EQ.0)THEN
   ZTOP(1:IDF%NCOL,1:IDF%NROW,ILAY)=XZTOP
   STRING='Constant value ['//TRIM(RTOS(XZTOP,'F',2))//'] for ZTOP L'//TRIM(ITOS(ILAY))
  ELSE
   STRING='Reading TOP L'//TRIM(ITOS(ILAY))//' '//TRIM(ITBFNAME(2,ILAY))//'...'
   IF(.NOT.IDFREADSCALE(ITBFNAME(2,ILAY),IDF,2,1,0.0D0,0))RETURN !## scale mean
   ZTOP(:,:,ILAY)=IDF%X
   !## adjust boundary whenever top
   DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
    IF(ZTOP(ICOL,IROW,ILAY).EQ.IDF%NODATA)IBOUND(ICOL,IROW,ILAY)=0
   ENDDO; ENDDO
  ENDIF
  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
 ENDDO

 !## read vertical coordinate data (BOT)
 DO ILAY=1,NLAY
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING=ITBFNAME(3,ILAY)
  READ(STRING,*,IOSTAT=IOS) XZBOT
  IF(IOS.EQ.0)THEN
   ZBOT(1:IDF%NCOL,1:IDF%NROW,ILAY)=XZBOT
   STRING='Constant value ['//TRIM(RTOS(XZBOT,'F',2))//'] for ZBOT L'//TRIM(ITOS(ILAY))
  ELSE
   STRING='Reading BOT L'//TRIM(ITOS(ILAY))//' '//TRIM(ITBFNAME(3,ILAY))//'...'
   IF(.NOT.IDFREADSCALE(ITBFNAME(3,ILAY),IDF,2,1,0.0D0,0))RETURN !## scale mean
   ZBOT(:,:,ILAY)=IDF%X   
   !## adjust boundary whenever top
   DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
    IF(ZBOT(ICOL,IROW,ILAY).EQ.IDF%NODATA)IBOUND(ICOL,IROW,ILAY)=0
   ENDDO; ENDDO
  ENDIF
  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
 ENDDO

 !## assign porosity
 DO ILAY=1,NLAY
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING=ITBFNAME(4,ILAY)
  READ(STRING,*,IOSTAT=IOS) APOR
  IF(IOS.EQ.0)THEN
   POR(1:IDF%NCOL,1:IDF%NROW,ILAY)=APOR
   STRING='Constant value ['//TRIM(RTOS(APOR,'F',2))//'] for POR_AQF L'//TRIM(ITOS(ILAY))
  ELSE
   STRING='Reading POR_AQF L'//TRIM(ITOS(ILAY))//' '//TRIM(ITBFNAME(4,ILAY))//'...'
   IF(.NOT.IDFREADSCALE(ITBFNAME(4,ILAY),IDF,2,1,0.0D0,0))RETURN !## scale mean
   POR(:,:,ILAY)=IDF%X   
   !## adjust boundary whenever por
   DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
    IF(POR(ICOL,IROW,ILAY).EQ.IDF%NODATA)IBOUND(ICOL,IROW,ILAY)=0
    IF(IBOUND(ICOL,IROW,ILAY).NE.0)THEN
     IF(POR(ICOL,IROW,ILAY).GE.1.0D0.OR.POR(ICOL,IROW,ILAY).LE.0.0D0)THEN
      IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You cannot specify a porosity of '//TRIM(RTOS(POR(ICOL,IROW,ILAY),'F',3))//CHAR(13)// &
         'It need to be always 0.0D0 > porosity < 1.0D0','Error')
      IF(IBATCH.EQ.1)THEN; WRITE(*,'(A)') 'You cannot specify a porosity of '//TRIM(RTOS(POR(ICOL,IROW,ILAY),'F',3))
                     WRITE(*,'(A)') 'It need to be always 0.0D0 > porosity < 1.0D0'; ENDIF
      RETURN
     ENDIF
    ENDIF
   ENDDO; ENDDO
  ENDIF
  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
 ENDDO

 DO ILAY=1,NLAY-1
  CALL WMESSAGEPEEK(ITYPE,MESSAGE)
  STRING=ITBFNAME(5,ILAY)
  READ(STRING,*,IOSTAT=IOS) APOR
  IF(IOS.EQ.0)THEN
   POR(1:IDF%NCOL,1:IDF%NROW,NLAY+ILAY)=APOR
   STRING='Constant value ['//TRIM(RTOS(APOR,'F',2))//'] for POR_AQT L'//TRIM(ITOS(ILAY))
  ELSE
   STRING='Reading POR_AQT L'//TRIM(ITOS(ILAY))//' '//TRIM(ITBFNAME(5,ILAY))//'...'
   IF(.NOT.IDFREADSCALE(ITBFNAME(5,ILAY),IDF,2,1,0.0D0,0))RETURN !## scale mean
   POR(:,:,NLAY+ILAY)=IDF%X   
   !## adjust boundary whenever por
   DO IROW=1,IDF%NROW; DO ICOL=1,IDF%NCOL
    IF(POR(ICOL,IROW,NLAY+ILAY).EQ.IDF%NODATA)THEN; IBOUND(ICOL,IROW,ILAY)=0; IBOUND(ICOL,IROW,ILAY+1)=0;ENDIF
   ENDDO; ENDDO
  ENDIF
  IF(IBATCH.EQ.1)WRITE(*,'(A)') TRIM(STRING)
  IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,STRING)
 ENDDO

 !## check top/bot consequentie - minimal thickness=minthickness
 NCONS=0
 DO IROW=1,IDF%NROW
  DO ICOL=1,IDF%NCOL
   DO ILAY=1,NLAY
    IF(IBOUND(ICOL,IROW,ILAY).NE.0)THEN
     IF(ZTOP(ICOL,IROW,ILAY).LT.ZBOT(ICOL,IROW,ILAY)+MINTHICKNESS)THEN
      !## only account for an error if they intersect
      IF(ZTOP(ICOL,IROW,ILAY).LT.ZBOT(ICOL,IROW,ILAY))NCONS=NCONS+1
      ZBOT(ICOL,IROW,ILAY)=ZTOP(ICOL,IROW,ILAY)-MINTHICKNESS
     ENDIF
     IF(ILAY.LT.NLAY)THEN
      IF(IBOUND(ICOL,IROW,ILAY+1).NE.0)THEN
       IF(ZBOT(ICOL,IROW,ILAY).LT.ZTOP(ICOL,IROW,ILAY+1))THEN !+MINTHICKNESS)THEN
        NCONS=NCONS+1
        ZTOP(ICOL,IROW,ILAY+1)=ZBOT(ICOL,IROW,ILAY)-MINTHICKNESS
       ENDIF
      ENDIF
     ENDIF
    ENDIF
   ENDDO
  ENDDO
 ENDDO

 IF(IBATCH.EQ.0)CALL WINDOWOUTSTATUSBAR(4,'')

! IF(.not.IDFALLOCATEX(IDF))stop
! DO ILAY=1,NLAY
!  idf%x=REAL(IBOUND(:,:,ilay))
!  LINE=TRIM(PREFVAL(1))//'\ibound_l'//TRIM(ITOS(ILAY))//'.idf'
!  IF(.not.idfwrite(idf,TRIM(LINE),1))then
!  endif
! END DO
! DO ILAY=1,NLAY
!  idf%x=ZTOP(:,:,ilay)
!  LINE=TRIM(PREFVAL(1))//'\ztop_l'//TRIM(ITOS(ILAY))//'.idf'
!  IF(.not.idfwrite(idf,TRIM(LINE),1))then
!  endif
! END DO
! DO ILAY=1,NLAY
!  idf%x=Zbot(:,:,ilay)
!  LINE=TRIM(PREFVAL(1))//'\zBOT_l'//TRIM(ITOS(ILAY))//'.idf'
!  IF(.not.idfwrite(idf,TRIM(LINE),1))then
!  endif
! END DO
! DO ILAY=1,(NLAY*2)-1
!  idf%x=POR(:,:,ilay)
!  LINE=TRIM(PREFVAL(1))//'\por_l'//TRIM(ITOS(ILAY))//'.idf'
!  IF(.not.idfwrite(idf,LINE,1))then
!  endif
! END DO

 TRACEDATIN=.TRUE.

 END FUNCTION TRACEDATIN

 !###====================================================================
 SUBROUTINE TRACEDELRC()
 !###====================================================================
 IMPLICIT NONE
 INTEGER :: I,II

 IF(IDF%IEQ.EQ.0)THEN

  !#fill global coordinates
  DELR(0)=0.0D0 
  DO I=1,IDF%NCOL
   DELR(I)=REAL(I)*IDF%DX
  END DO
  DO I=0,IDF%NCOL
   DELR(I)=IDF%XMIN+DELR(I)
  ENDDO
  
  DELC(0)=0.0D0 
  DO I=1,IDF%NROW
   DELC(I)=-REAL(I)*IDF%DY
  END DO
  DO I=0,IDF%NROW
   DELC(I)=IDF%YMAX+DELC(I)
  ENDDO
  
  !#fill local coordinates
  DO I=1,IDF%NCOL
   LDELR(I)=REAL(I)*IDF%DX
  END DO
  II=0
  DO I=IDF%NROW,1,-1 
   II=II+1
   LDELC(I)=REAL(II)*IDF%DY
  ENDDO

  DELX=IDF%DX
  DELY=IDF%DY

 ELSE

  DELR(0)=IDF%XMIN
  DO I=1,IDF%NCOL
   !## fill global coordinates
   DELR(I) =IDF%SX(I)
   !## delta x
   DELX(I) =DELR(I)-DELR(I-1)
   !## fill local coordinates
   LDELR(I)=DELR(I)-IDF%XMIN
  END DO
  DELC(0)=IDF%YMAX
  DO I=1,IDF%NROW
   !## fill global coordinates
   DELC(I) =IDF%SY(I)
   !## delta y
   DELY(I) =DELC(I-1)-DELC(I)
   !## fill local coordinates
   LDELC(I)=DELC(I-1)-IDF%YMIN
  END DO

 ENDIF
 
 END SUBROUTINE

END MODULE MOD_PLINES_TRACE