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

USE WINTERACTER
USE RESOURCE
USE IMODVAR
USE MOD_IDFPLOT
USE MODPLOT
USE MOD_UTL
USE MOD_OSD
USE MOD_IDF
USE MOD_SOLID_PCG
USE MOD_POLYGON 
USE MOD_POLYGON_DRAW 
USE MOD_POLYGON_UTL 
USE MOD_POLYGON_PAR
USE MOD_SOLID_PAR
USE MOD_SOLID_UTL
USE MOD_PROFILE 
USE MOD_PROFILE_PAR 
USE MOD_SOLID_PROFILE

USE MOD_3D, ONLY : IMOD3D_INIT
USE MOD_3D_UTL, ONLY : IMOD3D_CLOSE

USE MOD_MANAGER_UTL
USE MOD_COLOURS
USE MOD_MAIN_UTL

CHARACTER(LEN=256),PRIVATE :: SOLDIR
CHARACTER(LEN=256),DIMENSION(:),POINTER,PRIVATE :: REGISFILES

CONTAINS

 !###======================================================================
 SUBROUTINE SOLID_MAIN(ITYPE,MESSAGE)
 !###======================================================================
 IMPLICIT NONE
 TYPE(WIN_MESSAGE),INTENT(IN) :: MESSAGE
 INTEGER,INTENT(IN) :: ITYPE
 INTEGER :: IWIN
 LOGICAL :: LEX
 REAL(KIND=DP_KIND),PARAMETER :: ZOFFSET=0.0D0
 
 SELECT CASE(MESSAGE%WIN)

  CASE (ID_DSOLID)
   SELECT CASE(ITYPE)
    CASE(PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      CASE (ID_FEED)
       CALL SOLID_CLEANMANAGER() !## clean manager
       IF(SOLID_FILLMANAGER())THEN
        !## zoom full extent
        CALL IDFZOOM(ID_ZOOMFULLMAP,(MPW%XMAX+MPW%XMIN)/2.0D0,(MPW%YMAX+MPW%YMIN)/2.0D0,0)
        CALL IDFPLOTFAST(1)
       ENDIF
      CASE (IDCANCEL)
       CALL SOLID_CLOSE()
!       CALL IDFZOOM(ID_DGOTOXY,0.0D0,0.0D0,0)
       CALL IDFPLOT(1)

      CASE (IDHELP)
       CALL UTL_GETHELP('5.4','TMO.SolTool')               
     END SELECT
    CASE(TABCHANGED)
     SELECT CASE (MESSAGE%VALUE2)
      !## new tab is solids
      CASE(ID_DSOLIDTAB1)

       !## read entire sol, incl. spf-files
       IF(SOLIDOPENSOL('W',GETSOLNAME(),IQ=1,TXT='Polygons'))THEN
       ENDIF

       CALL POLYGON1DRAWSHAPE(1,SHP%NPOL)
       !## reread sol file again
       IF(SOLIDOPENSOL('R',GETSOLNAME()))THEN
        CALL WDIALOGSELECT(ID_DSOLIDTAB2)
        CALL WDIALOGCLEARFIELD(IDF_MENU1)
        IF(SHP%NPOL.GT.0)THEN
         SHP%POL%IACT=0
         CALL WDIALOGPUTMENU(IDF_MENU1,SHP%POL%PNAME,SHP%NPOL,SHP%POL%IACT)
        ENDIF
        CALL POLYGON1DRAWSHAPE(1,SHP%NPOL)
        CALL WDIALOGSELECT(ID_DSOLIDTAB4)
        CALL WDIALOGCLEARFIELD(IDF_GRID1)
        IF(NGENSOL.GT.0)THEN
         CALL WGRIDPUTINTEGER(IDF_GRID1,1,GENSOL%ILAY,NGENSOL)
         CALL WGRIDPUTSTRING(IDF_GRID1,2,GENSOL%FNAME,NGENSOL)
        ENDIF
        CALL WDIALOGSELECT(ID_DSOLIDTAB5)
        CALL WDIALOGCLEARFIELD(IDF_GRID1)
        IF(NIPFP.GT.0)THEN
         CALL WGRIDPUTINTEGER(IDF_GRID1,1,IPFP%ILAY,NIPFP)
         CALL WGRIDPUTSTRING(IDF_GRID1,2,IPFP%FNAME,NIPFP)
        ENDIF
        CALL WDIALOGSELECT(ID_DSOLIDTAB3)
        CALL WDIALOGCLEARFIELD(IDF_MENU1)
        IF(NMASK.GT.0)CALL WDIALOGPUTMENU(IDF_MENU1,MASK%ALIAS,NMASK,1)
       ENDIF

     END SELECT
   END SELECT

  CASE (ID_DSOLIDTAB1)
   SELECT CASE(ITYPE)

    CASE(FIELDCHANGED)
     SELECT CASE (MESSAGE%VALUE2)
      !## select different solid in menu
      CASE (IDF_MENU1)
       CALL SOLID_READSOL()

     END SELECT

    CASE(PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      CASE (ID_NEW)
       IF(SOLID_NEW())CALL SOLID_READSOL()
      CASE (ID_PROFILE)
       CALL SOLID_PROFILE()
      CASE (ID_INFO)
       CALL WDIALOGSELECT(ID_DSOLIDTAB1)
       CALL WDIALOGGETMENU(IDF_MENU1,IWIN,FNAME)
       INQUIRE(FILE=GETSOLNAME(),EXIST=LEX)
       IF(LEX)THEN
        IWIN=0
        CALL WINDOWOPENCHILD(IWIN,FLAGS=SYSMENUON,WIDTH=1000,HEIGHT=500)
        CALL WINDOWSELECT(IWIN)
        CALL WEDITFILE(GETSOLNAME(),ITYPE=MODAL,IDMENU=0,IFONT=COURIERNEW,ISIZE=10)
       ELSE
        CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'iMOD cannot find/open:'//CHAR(13)//TRIM(GETSOLNAME()),'Error')
       ENDIF
      CASE (ID_3D)
       CALL SOLID_3D()
      CASE (ID_CALCULATE)
       CALL SOLID_CALC()
       CALL WINDOWSELECT(0); CALL WINDOWOUTSTATUSBAR(4,''); CALL WINDOWOUTSTATUSBAR(3,'')
      !## remove solid (incl. clean manager)
      CASE (ID_DELETE)
       CALL WDIALOGSELECT(ID_DSOLIDTAB1)
       CALL WDIALOGGETMENU(IDF_MENU1,IWIN,FNAME)
       CALL SOLID_CLEANMANAGER() !## clean manager
       IF(UTL_DEL1TREE(TRIM(PREFVAL(1))//'\SOLIDS\'//TRIM(FNAME)))THEN
        CALL SOLID_FIELDS('')
       ENDIF
     END SELECT
   END SELECT

  CASE (ID_DSOLIDTAB2)
   !## check polygon actions
   IACTSHAPES=(/3,3,1,3,3,3/)
   CALL POLYGON1MAIN(ITYPE,MESSAGE)
   IF(ITYPE.EQ.PUSHBUTTON.AND.MESSAGE%VALUE1.EQ.ID_ZOOMSELECT)THEN
    CALL IDFZOOM(ID_DGOTOXY,0.0D0,0.0D0,0)
    CALL IDFPLOT(1)
   ENDIF

  !## check masks
  CASE (ID_DSOLIDTAB3)
   SELECT CASE(ITYPE)

    CASE(FIELDCHANGED)
     SELECT CASE (MESSAGE%VALUE2)
     END SELECT

    CASE(PUSHBUTTON)
     SELECT CASE (MESSAGE%VALUE1)
      CASE (ID_NEW)
       IF(SOLID_NEWMASKS(ZOFFSET))THEN; ENDIF
      CASE (ID_OPENIDF)
      
      CASE (ID_INFO)
       CALL WDIALOGSELECT(ID_DSOLIDTAB3)
       CALL WDIALOGGETMENU(IDF_MENU1,IWIN,FNAME)
      CASE (ID_DELETE)
       CALL WDIALOGSELECT(ID_DSOLIDTAB3)
       CALL WDIALOGGETMENU(IDF_MENU1,IWIN,FNAME)
       CALL SOLID_FIELDS('')
     END SELECT
   END SELECT

 END SELECT

 END SUBROUTINE SOLID_MAIN

 !###======================================================================
 LOGICAL FUNCTION SOLID_WVP()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J,N,IROW,ICOL
 REAL(KIND=DP_KIND) :: T,B
 LOGICAL :: LEX
 
 SOLID_WVP=.FALSE.

 !## use idf settings from first top() idf
 IF(MDLIDF%XMIN.EQ.MDLIDF%XMAX)THEN
  IF(.NOT.IDFREAD(TOPIDF(1),TOPIDF(1)%FNAME,0))THEN; RETURN; ENDIF
  CLOSE(TOPIDF(1)%IU); CALL IDFCOPY(TOPIDF(1),MDLIDF)
  MDLIDF%FNAME=TOPIDF(1)%FNAME
 ENDIF
 
 DO I=1,SIZE(TOPIDF)
  WRITE(*,'(A)') 'Reading '//TRIM(TOPIDF(I)%FNAME)
  CALL IDFCOPY(MDLIDF,TOPIDF(I))
  IF(.NOT.IDFREADSCALE(TOPIDF(I)%FNAME,TOPIDF(I),2,1,0.0D0,0))RETURN !## scale mean
 ENDDO
 DO I=1,SIZE(BOTIDF)
  WRITE(*,'(A)') 'Reading '//TRIM(BOTIDF(I)%FNAME)
  CALL IDFCOPY(MDLIDF,BOTIDF(I))
  IF(.NOT.IDFREADSCALE(BOTIDF(I)%FNAME,BOTIDF(I),2,1,0.0D0,0))RETURN !## scale mean
 ENDDO
 DO I=1,SIZE(CIDF)
  WRITE(*,'(A)') 'Reading '//TRIM(CIDF(I)%FNAME)
  CALL IDFCOPY(MDLIDF,CIDF(I))
  IF(.NOT.IDFREADSCALE(CIDF(I)%FNAME,CIDF(I),6,1,0.0D0,0))RETURN !## scale mean
 ENDDO

 N=SIZE(TOPIDF)
 DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
  J=1; T=TOPIDF(J)%X(ICOL,IROW); B=BOTIDF(J)%X(ICOL,IROW)
  DO I=2,SIZE(TOPIDF)
   !## no thickness available and agrees minimal resistance
   LEX=TOPIDF(I)%X(ICOL,IROW).GE.BOTIDF(I-1)%X(ICOL,IROW)
   !## resistance less than minimal value (minc)
   IF(.NOT.LEX)LEX=CIDF(I-1)%X(ICOL,IROW).LE.MINC
   !## aquifer continues vertically
   IF(LEX)THEN
    B=BOTIDF(I)%X(ICOL,IROW)
   ELSE
    TOPIDF(J)%X(ICOL,IROW)=T; BOTIDF(J)%X(ICOL,IROW)=B
    J=J+1; T=TOPIDF(I)%X(ICOL,IROW); B=BOTIDF(I)%X(ICOL,IROW)
   ENDIF
  ENDDO
  TOPIDF(J)%X(ICOL,IROW)=T; BOTIDF(J)%X(ICOL,IROW)=B
  !## blank out remaining
  DO I=J+1,N; TOPIDF(I)%X(ICOL,IROW)=TOPIDF(I)%NODATA; BOTIDF(I)%X(ICOL,IROW)=BOTIDF(I)%NODATA; ENDDO
 ENDDO; ENDDO
  
 DO I=1,SIZE(TOPIDF)
  TOPIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\top_wvp_l'//TRIM(ITOS(I))//'.idf'
  IF(IDFWRITE(TOPIDF(I),TOPIDF(I)%FNAME,1))THEN; ENDIF
  BOTIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\bot_wvp_l'//TRIM(ITOS(I))//'.idf'
  IF(IDFWRITE(BOTIDF(I),BOTIDF(I)%FNAME,1))THEN; ENDIF
 ENDDO
 
 CALL IDFDEALLOCATE(TOPIDF,SIZE(TOPIDF)); DEALLOCATE(TOPIDF)
 CALL IDFDEALLOCATE(BOTIDF,SIZE(BOTIDF)); DEALLOCATE(BOTIDF)
 CALL IDFDEALLOCATE(CIDF  ,SIZE(CIDF  )); DEALLOCATE(CIDF  )
 
 SOLID_WVP=.TRUE.
 
 END FUNCTION SOLID_WVP

 !###======================================================================
 LOGICAL FUNCTION SOLID_GEOTOP()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J,K,IROW,ICOL,ILAY,JLAY,IL,IOS
 REAL(KIND=DP_KIND) :: C1,C2,C3,X
 
 SOLID_GEOTOP=.FALSE.
 
 CALL UTL_CREATEDIR(RESULTFOLDER)
 CALL UTL_CREATEDIR(TRIM(RESULTFOLDER)//'\top')
 CALL UTL_CREATEDIR(TRIM(RESULTFOLDER)//'\bot')
 CALL UTL_CREATEDIR(TRIM(RESULTFOLDER)//'\kdw')
 CALL UTL_CREATEDIR(TRIM(RESULTFOLDER)//'\vcw')
 !CALL UTL_CREATEDIR(TRIM(RESULTFOLDER)//'\bnd') 
 !CALL UTL_CREATEDIR(TRIM(RESULTFOLDER)//'\shd') 
 
 !## use idf settings from first top() idf
 IF(MDLIDF%XMIN.EQ.MDLIDF%XMAX)THEN
  IF(.NOT.IDFREAD(MDLIDF,TPM(1)%FNAME,0))THEN; RETURN; ENDIF; CLOSE(MDLIDF%IU)
 ELSE
  CALL UTL_IDFSNAPTOGRID(MDLIDF%XMIN,MDLIDF%XMAX,MDLIDF%YMIN,MDLIDF%YMAX,MDLIDF%DX,MDLIDF%NCOL,MDLIDF%NROW)
 ENDIF
 
 !## reading starting head original model
 !!DO I=1,SIZE(SHM)
 !! WRITE(*,'(A)') 'Reading '//TRIM(SHM(I)%FNAME); CALL IDFCOPY(MDLIDF,SHM(I))
 !! READ(SHM(I)%FNAME,*,IOSTAT=IOS) X
 !! IF(IOS.EQ.0)THEN
 !!  SHM(I)%X=X
 !! ELSE
 !!  IF(.NOT.IDFREADSCALE(SHM(I)%FNAME,SHM(I),2,1,0.0D0,0))RETURN !## scale mean
 !! ENDIF
 !!ENDDO
 
 !## reading top original model
 DO I=1,SIZE(TPM)
  WRITE(*,'(A)') 'Reading '//TRIM(TPM(I)%FNAME); CALL IDFCOPY(MDLIDF,TPM(I))
  READ(TPM(I)%FNAME,*,IOSTAT=IOS) X
  IF(IOS.EQ.0)THEN
   TPM(I)%X=X
  ELSE
   IF(.NOT.IDFREADSCALE(TPM(I)%FNAME,TPM(I),2,1,0.0D0,0))RETURN !## scale mean
  ENDIF
 ENDDO
 !## reading bottom original model
 DO I=1,SIZE(BTM)
  WRITE(*,'(A)') 'Reading '//TRIM(BTM(I)%FNAME); CALL IDFCOPY(MDLIDF,BTM(I))
  READ(BTM(I)%FNAME,*,IOSTAT=IOS) X
  IF(IOS.EQ.0)THEN
   BTM(I)%X=X
  ELSE
   IF(.NOT.IDFREADSCALE(BTM(I)%FNAME,BTM(I),2,1,0.0D0,0))RETURN !## scale mean
  ENDIF
 ENDDO
 !## reading boundary condition original model
 !DO I=1,SIZE(IBM)
 ! WRITE(*,'(A)') 'Reading '//TRIM(IBM(I)%FNAME); CALL IDFCOPY(MDLIDF,IBM(I))
 ! READ(IBM(I)%FNAME,*,IOSTAT=IOS) X
 ! IF(IOS.EQ.0)THEN
 !  IBM(I)%X=X
 ! ELSE
 !  IF(.NOT.IDFREADSCALE(IBM(I)%FNAME,IBM(I),1,1,0.0D0,0))RETURN !## scale boundary
 ! ENDIF
 !ENDDO
 !## reading horizontal permeability original model
 DO I=1,SIZE(KHM)
  WRITE(*,'(A)') 'Reading '//TRIM(KHM(I)%FNAME); CALL IDFCOPY(MDLIDF,KHM(I))
  READ(KHM(I)%FNAME,*,IOSTAT=IOS) X
  IF(IOS.EQ.0)THEN
   IF(.NOT.IDFALLOCATEX(KHM(I)))RETURN
   KHM(I)%X=X
  ELSE
   IF(.NOT.IDFREADSCALE(KHM(I)%FNAME,KHM(I),3,1,0.0D0,0))RETURN !## scale geometric
  ENDIF
 ENDDO
 !## reading vertical anisotropy original model
 !DO I=1,SIZE(KAM)
 ! WRITE(*,'(A)') 'Reading '//TRIM(KAM(I)%FNAME); CALL IDFCOPY(MDLIDF,KAM(I))
 ! READ(KAM(I)%FNAME,*,IOSTAT=IOS) X
 ! IF(IOS.EQ.0)THEN
 !  IF(.NOT.IDFALLOCATEX(KAM(I)))RETURN
 !  KAM(I)%X=X
 ! ELSE
 !  IF(.NOT.IDFREADSCALE(KAM(I)%FNAME,KAM(I),2,1,0.0D0,0))RETURN !## scale mean
 ! ENDIF
 !ENDDO
 !## reading vertical permeability original model
 DO I=1,NLAYM-1 
  WRITE(*,'(A)') 'Reading '//TRIM(KVM(I)%FNAME); CALL IDFCOPY(MDLIDF,KVM(I))
  READ(KVM(I)%FNAME,*,IOSTAT=IOS) X
  IF(IOS.EQ.0)THEN
   IF(.NOT.IDFALLOCATEX(KVM(I)))RETURN
   KVM(I)%X=X
  ELSE
   IF(.NOT.IDFREADSCALE(KVM(I)%FNAME,KVM(I),3,1,0.0D0,0))RETURN !## scale geometric
  ENDIF
 ENDDO
 CALL IDFNULLIFY(KVM(NLAYM)); CALL IDFCOPY(KVM(1),KVM(NLAYM))
 
 !## reading KH-GeoTOP
 DO I=1,SIZE(KHG)
  WRITE(*,'(A)') 'Reading '//TRIM(KHG(I)%FNAME); CALL IDFCOPY(MDLIDF,KHG(I))
  READ(KHG(I)%FNAME,*,IOSTAT=IOS) X
  IF(IOS.EQ.0)THEN
   IF(.NOT.IDFALLOCATEX(KHG(I)))RETURN
   KHG(I)%X=X
  ELSE
   IF(.NOT.IDFREADSCALE(KHG(I)%FNAME,KHG(I),3,1,0.0D0,0))RETURN !## scale geometric
  ENDIF
 ENDDO
 !## reading KV-GeoTOP
 DO I=1,SIZE(KVG)
  WRITE(*,'(A)') 'Reading '//TRIM(KVG(I)%FNAME); CALL IDFCOPY(MDLIDF,KVG(I))
  READ(KVG(I)%FNAME,*,IOSTAT=IOS) X
  IF(IOS.EQ.0)THEN
   IF(.NOT.IDFALLOCATEX(KVG(I)))RETURN
   KVG(I)%X=X
  ELSE
   IF(.NOT.IDFREADSCALE(KVG(I)%FNAME,KVG(I),3,1,0.0D0,0))RETURN !!## scale geometric
  ENDIF
 ENDDO
 !ALLOCATE(IBG(NLAYG),SHG(NLAYG))
 !DO I=1,NLAYG; CALL IDFNULLIFY(IBG(I)); CALL IDFCOPY(KVG(1),IBG(I)); ENDDO
 !DO I=1,NLAYG; CALL IDFNULLIFY(SHG(I)); CALL IDFCOPY(KVG(1),SHG(I)); ENDDO

 !## determine lower boundary of GeoTOP (voxels)
 DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
  !## fill in remainder of GeoTOP-voxels on the bottom
  DO I=SIZE(KHG),1,-1
   IF(KHG(I)%X(ICOL,IROW).GT.0.0D0)EXIT
   DO J=1,SIZE(TPM)
    IF(TPM(J)%X(ICOL,IROW).GE.KHG(I)%BOT.AND.BTM(J)%X(ICOL,IROW).LE.KHG(I)%TOP)THEN
     KHG(I)%X(ICOL,IROW)=KHM(J)%X(ICOL,IROW)
     !KVG(I)%X(ICOL,IROW)=KHM(J)%X(ICOL,IROW)*KAM(J)%X(ICOL,IROW)
     KVG(I)%X(ICOL,IROW)=KHM(J)%X(ICOL,IROW)*1.0D0
    ENDIF
   ENDDO
  ENDDO
  !## fill in startingheads/boundary condition for GeoTOP voxels
  !DO I=1,SIZE(KHG)
  ! DO J=1,SIZE(TPM)
  !  TP=TPM(J)%X(ICOL,IROW)
  !  IF(J.LT.SIZE(TPM))THEN; BT=TPM(J+1)%X(ICOL,IROW); ELSE; BT=BTM(J)%X(ICOL,IROW); ENDIF
  !  IF(TP.GE.KHG(I)%BOT.AND.BT.LE.KHG(I)%TOP)THEN
  !   SHG(I)%X(ICOL,IROW)=SHM(J)%X(ICOL,IROW); SHG(I)%NODATA=SHM(J)%NODATA
  !   IBG(I)%X(ICOL,IROW)=IBM(J)%X(ICOL,IROW)
  !  ENDIF
  ! ENDDO
  !ENDDO
  
  !## correct top/bottoms for original model for GeoTOP
  I=NLAYG
  DO J=1,SIZE(TPM)
   !## adjust appropriate top of original model
   IF(TPM(J)%X(ICOL,IROW).GT.KHG(I)%BOT)TPM(J)%X(ICOL,IROW)=KHG(I)%BOT
   IF(BTM(J)%X(ICOL,IROW).GT.KHG(I)%BOT)BTM(J)%X(ICOL,IROW)=KHG(I)%BOT
   !## make khm.eq.0.0D0 for zero thicknesses
   IF(TPM(J)%X(ICOL,IROW)-BTM(J)%X(ICOL,IROW).LE.0.0D0)THEN
    KHM(J)%X(ICOL,IROW)=0.0D0; KVM(J)%X(ICOL,IROW)=0.0D0 !; KAM(J)%X(ICOL,IROW)=0.0D0
   ENDIF
  ENDDO
 ENDDO; ENDDO
 
 IF(.NOT.IDFALLOCATEX(MDLIDF))RETURN; ALLOCATE(IBND(MDLIDF%NCOL,MDLIDF%NROW))
 
 !## search first active modellayer
 ILAY=1
ILLOOP: DO
  DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL; IF(KHG(ILAY)%X(ICOL,IROW).GT.0.0D0)EXIT ILLOOP; ENDDO; ENDDO
  ILAY=ILAY+1; IF(ILAY.GT.NLAYG)EXIT
 ENDDO ILLOOP
 ILAY=ILAY-1

 WRITE(*,'(A,I3,A)') 'Start with layer ',ILAY+1,' from GeoTOP, assigned first active modellayer'
 
 !## continue until finished
 JLAY=0; IL=0; DO
  ILAY=ILAY+1; IL=IL+1
  IF(JLAY.GE.NLAYM)EXIT
  
  !## process geotop files
  IF(ILAY.LE.NLAYG)THEN

   !## write top
   MDLIDF%X=KHG(ILAY)%TOP; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\top\top_l'//TRIM(ITOS(IL))//'.idf'
   MDLIDF%NODATA=KHG(ILAY)%NODATA; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !## write bottom
   MDLIDF%X=KHG(ILAY)%BOT; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\bot\bot_l'//TRIM(ITOS(IL))//'.idf'
   MDLIDF%NODATA=KHG(ILAY)%NODATA; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !## compute transmissivity
   MDLIDF%NODATA=0.0D0
   DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
    MDLIDF%X(ICOL,IROW)=(KHG(ILAY)%TOP-KHG(ILAY)%BOT)*KHG(ILAY)%X(ICOL,IROW)
    IF(KHG(ILAY)%X(ICOL,IROW).EQ.KHG(ILAY)%NODATA)MDLIDF%X(ICOL,IROW)=MDLIDF%NODATA
   ENDDO; ENDDO
   MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\kdw\kdw_l'//TRIM(ITOS(IL))//'.idf'
   MDLIDF%ITB=1; MDLIDF%TOP=KHG(ILAY)%TOP ; MDLIDF%BOT=KHG(ILAY)%BOT
   MDLIDF%NODATA=0.0D0; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !## compute vertical resistance for previous modellayer
   IF(ILAY.GT.1)THEN
    DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
     C1=0.0D0; IF(KVG(ILAY-1)%X(ICOL,IROW).GT.0.0D0)C1=(KHG(ILAY-1)%TOP-KHG(ILAY-1)%BOT)/KVG(ILAY-1)%X(ICOL,IROW)
     C2=0.0D0; IF(KVG(ILAY  )%X(ICOL,IROW).GT.0.0D0)C2=(KHG(ILAY  )%TOP-KHG(ILAY  )%BOT)/KVG(ILAY  )%X(ICOL,IROW)
     MDLIDF%X(ICOL,IROW)=0.5*C1+0.5*C2
    ENDDO; ENDDO
    MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\vcw\vcw_l'//TRIM(ITOS(IL-1))//'.idf'
    MDLIDF%ITB=1; MDLIDF%TOP=(KHG(ILAY-1)%TOP+KHG(ILAY-1)%BOT)/2.0 ; MDLIDF%BOT=(KHG(ILAY)%TOP+KHG(ILAY)%BOT)/2.0
    MDLIDF%NODATA=0.0D0; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   ENDIF

   !## compute boundary based on GeoTop
   !IBND=INT(0,1)
   !DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
   ! IF(KHG(ILAY)%X(ICOL,IROW).GT.0.0D0.AND.IBG(ILAY)%X(ICOL,IROW).GT.0)IBND(ICOL,IROW)=INT(1,1)
   !ENDDO; ENDDO
   !DO ICOL=1,MDLIDF%NCOL
   ! IF(KHG(ILAY)%X(ICOL,1          ).GT.1.0D0.AND.IBND(ICOL,1          ).GT.0)IBND(ICOL,1          )=INT(-1,1)
   ! IF(KHG(ILAY)%X(ICOL,MDLIDF%NROW).GT.1.0D0.AND.IBND(ICOL,MDLIDF%NROW).GT.0)IBND(ICOL,MDLIDF%NROW)=INT(-1,1)
   !ENDDO
   !DO IROW=1,MDLIDF%NROW
   ! IF(KHG(ILAY)%X(1,IROW          ).GT.1.0D0.AND.IBND(1,IROW          ).GT.0)IBND(1,IROW          )=INT(-1,1)
   ! IF(KHG(ILAY)%X(MDLIDF%NCOL,IROW).GT.1.0D0.AND.IBND(MDLIDF%NCOL,IROW).GT.0)IBND(MDLIDF%NCOL,IROW)=INT(-1,1)
   !ENDDO
   !## make sure that for constant heads, true heads are available
   !DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
   ! IF(IBND(ICOL,IROW).LT.INT(0,1).AND.SHG(ILAY)%X(ICOL,IROW).EQ.SHG(ILAY)%NODATA)IBND(ICOL,IROW)=INT(1,1)
   ! IF(IBND(ICOL,IROW).GT.INT(1,1).AND.SHG(ILAY)%X(ICOL,IROW).EQ.SHG(ILAY)%NODATA)SHG(ILAY)%X(ICOL,IROW)=0.0D0
   !ENDDO; ENDDO
   
   !MDLIDF%NODATA=0.0D0; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\bnd\bnd_l'//TRIM(ITOS(IL))//'.idf'
   !MDLIDF%ITB=1; MDLIDF%TOP=KHG(ILAY)%TOP ; MDLIDF%BOT=KHG(ILAY)%BOT
   !MDLIDF%X=REAL(IBND); IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !MDLIDF%NODATA=SHG(ILAY)%NODATA; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\shd\shd_l'//TRIM(ITOS(IL))//'.idf'
   !MDLIDF%X=SHG(ILAY)%X; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
  
  ELSE

   !## search first active modellayer
   IF(JLAY.EQ.0)THEN
    DO; IF(SUM(KHM(JLAY+1)%X).GT.0.0D0)EXIT; JLAY=JLAY+1; ENDDO
    WRITE(*,'(A,I3,A)') 'Start with layer ',JLAY+1,' from ORIGINAL model (previous layers filled in by GeoTOP'
   ENDIF
   JLAY=JLAY+1  

   !## write top
   MDLIDF%X=TPM(JLAY)%X; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\top\top_l'//TRIM(ITOS(IL))//'.idf'
   MDLIDF%ITB=0; MDLIDF%NODATA=TPM(JLAY)%NODATA; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !## write bottom 
   MDLIDF%X=BTM(JLAY)%X; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\bot\bot_l'//TRIM(ITOS(IL))//'.idf'
   MDLIDF%NODATA=BTM(JLAY)%NODATA; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !## compute transmissivity
   DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
    MDLIDF%X(ICOL,IROW)=(TPM(JLAY)%X(ICOL,IROW)-BTM(JLAY)%X(ICOL,IROW))*KHM(JLAY)%X(ICOL,IROW)
   ENDDO; ENDDO
   MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\kdw\kdw_l'//TRIM(ITOS(IL))//'.idf'
   MDLIDF%NODATA=KHM(JLAY)%NODATA; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !## compute vertical resistance
   DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
    IF(ILAY-1.EQ.NLAYG)THEN
     C1=0.0D0; K=KVG(ILAY-1)%X(ICOL,IROW); IF(K.GT.0.0D0)C1=(KHG(ILAY-1)%TOP-KHG(ILAY-1)%BOT)/K
    ELSE
     !C1=0.0D0; K=KHM(JLAY-1)%X(ICOL,IROW)*KAM(JLAY-1)%X(ICOL,IROW)
     C1=0.0D0; K=KHM(JLAY-1)%X(ICOL,IROW)*1.0D0
     IF(K.GT.0.0D0)C1=(TPM(JLAY-1)%X(ICOL,IROW)-BTM(JLAY-1)%X(ICOL,IROW))/K
    ENDIF
    !C2=0.0D0; K=KHM(JLAY)%X(ICOL,IROW)*KAM(JLAY)%X(ICOL,IROW)
    C2=0.0D0; K=KHM(JLAY)%X(ICOL,IROW)*1.0D0
    IF(K.GT.0.0D0)C2=(TPM(JLAY)%X(ICOL,IROW)-BTM(JLAY)%X(ICOL,IROW))/K
    C3=0.0D0; K=KVM(JLAY)%X(ICOL,IROW)
    IF(K.GT.0.0D0)C3=(BTM(JLAY-1)%X(ICOL,IROW)-TPM(JLAY)%X(ICOL,IROW))/K
    MDLIDF%X(ICOL,IROW)=0.5*C1+0.5*C2+C3
   ENDDO; ENDDO
   MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\vcw\vcw_l'//TRIM(ITOS(IL-1))//'.idf'
   MDLIDF%NODATA=0.0D0; IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN

   !## compute boundary, constant head around (optional)
   !IBND=INT(0,1)
   !DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
   ! IF(IBM(JLAY)%X(ICOL,IROW).GT.0)IBND(ICOL,IROW)=INT(1,1)
   !ENDDO; ENDDO
   !DO ICOL=1,MDLIDF%NCOL
   ! IF(KHM(JLAY)%X(ICOL,1          ).GT.1.0D0.AND.IBND(ICOL,1          ).GT.0)IBND(ICOL,1          )=INT(-1,1)
   ! IF(KHM(JLAY)%X(ICOL,MDLIDF%NROW).GT.1.0D0.AND.IBND(ICOL,MDLIDF%NROW).GT.0)IBND(ICOL,MDLIDF%NROW)=INT(-1,1)
   !ENDDO
   !DO IROW=1,MDLIDF%NROW
   ! IF(KHM(JLAY)%X(1,IROW          ).GT.1.0D0.AND.IBND(1,IROW          ).GT.0)IBND(1,IROW          )=INT(-1,1)
   ! IF(KHM(JLAY)%X(MDLIDF%NCOL,IROW).GT.1.0D0.AND.IBND(MDLIDF%NCOL,IROW).GT.0)IBND(MDLIDF%NCOL,IROW)=INT(-1,1)
   !ENDDO
   !!## make sure that for constant heads, true heads are available
   !DO IROW=1,MDLIDF%NROW; DO ICOL=1,MDLIDF%NCOL
   ! IF(IBND(ICOL,IROW).LT.INT(0,1).AND.SHM(JLAY)%X(ICOL,IROW).EQ.SHM(JLAY)%NODATA)IBND(ICOL,IROW)=INT(1,1)
   ! IF(IBND(ICOL,IROW).GT.INT(1,1).AND.SHM(JLAY)%X(ICOL,IROW).EQ.SHM(JLAY)%NODATA)SHM(JLAY)%X(ICOL,IROW)=0.0D0
   !ENDDO; ENDDO

   !!## writing boundary conditions
   !MDLIDF%NODATA=0.0D0; MDLIDF%X=REAL(IBND); MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\bnd\bnd_l'//TRIM(ITOS(IL))//'.idf'
   !IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN
   !!## write starting heads
   !MDLIDF%NODATA=SHM(JLAY)%NODATA; MDLIDF%X=SHM(JLAY)%X; MDLIDF%FNAME=TRIM(RESULTFOLDER)//'\shd\shd_l'//TRIM(ITOS(IL))//'.idf'
   !IF(.NOT.IDFWRITE(MDLIDF,MDLIDF%FNAME,1))RETURN

  ENDIF

 ENDDO 

 !IL=IL-1; IU=UTL_GETUNIT()
 !OPEN(IU,FILE=TRIM(RESULTFOLDER)//'\geotop.run',STATUS='UNKNOWN',ACTION='WRITE')
 !WRITE(IU,'(I10,A)') IL,',(BND)'
 !DO I=1,IL
 ! LINE=',1.0D0,0.0D0,'//TRIM(RESULTFOLDER)//'\BND\BND_L'//TRIM(ITOS(I))//'.IDF'
 ! WRITE(IU,'(I10,A)') I,TRIM(LINE)
 !ENDDO
 !WRITE(IU,'(I10,A)') IL,',(SHD)'
 !DO I=1,IL
 ! LINE=',1.0D0,0.0D0,'//TRIM(RESULTFOLDER)//'\SHD\SHD_L'//TRIM(ITOS(I))//'.IDF'
 ! WRITE(IU,'(I10,A)') I,TRIM(LINE)
 !ENDDO
 !WRITE(IU,'(I10,A)') IL,',(TOP)'
 !DO I=1,IL
 ! LINE=',1.0D0,0.0D0,'//TRIM(RESULTFOLDER)//'\TOP\TOP_L'//TRIM(ITOS(I))//'.IDF'
 ! WRITE(IU,'(I10,A)') I,TRIM(LINE)
 !ENDDO
 !WRITE(IU,'(I10,A)') IL,',(BOT)'
 !DO I=1,IL
 ! LINE=',1.0D0,0.0D0,'//TRIM(RESULTFOLDER)//'\BOT\BOT_L'//TRIM(ITOS(I))//'.IDF'
 ! WRITE(IU,'(I10,A)') I,TRIM(LINE)
 !ENDDO
 !WRITE(IU,'(I10,A)') IL,',(KDW)'
 !DO I=1,IL
 ! LINE=',1.0D0,0.0D0,'//TRIM(RESULTFOLDER)//'\KDW\KDW_L'//TRIM(ITOS(I))//'.IDF'
 ! WRITE(IU,'(I10,A)') I,TRIM(LINE)
 !ENDDO
 !WRITE(IU,'(I10,A)') IL-1,',(VCW)'
 !DO I=1,IL-1
 ! LINE=',1.0D0,0.0D0,'//TRIM(RESULTFOLDER)//'\VCW\VCW_L'//TRIM(ITOS(I))//'.IDF'
 ! WRITE(IU,'(I10,A)') I,TRIM(LINE)
 !ENDDO
 !CLOSE(IU)
 
 SOLID_GEOTOP=.TRUE.
 
 END FUNCTION SOLID_GEOTOP

 !###======================================================================
 SUBROUTINE SOLID_GEOTOP_DEALLOCATE()
 !###======================================================================
 IMPLICIT NONE
 
 IF(ALLOCATED(KHG))THEN; CALL IDFDEALLOCATE(KHG,SIZE(KHG)); DEALLOCATE(KHG); ENDIF
 IF(ALLOCATED(KVG))THEN; CALL IDFDEALLOCATE(KVG,SIZE(KVG)); DEALLOCATE(KVG); ENDIF
 IF(ALLOCATED(SHG))THEN; CALL IDFDEALLOCATE(SHG,SIZE(SHG)); DEALLOCATE(SHG); ENDIF
 IF(ALLOCATED(IBG))THEN; CALL IDFDEALLOCATE(IBG,SIZE(IBG)); DEALLOCATE(IBG); ENDIF
 IF(ALLOCATED(KHM))THEN; CALL IDFDEALLOCATE(KHM,SIZE(KHM)); DEALLOCATE(KHM); ENDIF
 IF(ALLOCATED(KVM))THEN; CALL IDFDEALLOCATE(KVM,SIZE(KVM)); DEALLOCATE(KVM); ENDIF
 IF(ALLOCATED(KAM))THEN; CALL IDFDEALLOCATE(KAM,SIZE(KAM)); DEALLOCATE(KAM); ENDIF
 IF(ALLOCATED(SHM))THEN; CALL IDFDEALLOCATE(SHM,SIZE(SHM)); DEALLOCATE(SHM); ENDIF
 IF(ALLOCATED(IBM))THEN; CALL IDFDEALLOCATE(IBM,SIZE(IBM)); DEALLOCATE(IBM); ENDIF
 IF(ALLOCATED(TPM))THEN; CALL IDFDEALLOCATE(TPM,SIZE(TPM)); DEALLOCATE(TPM); ENDIF
 IF(ALLOCATED(BTM))THEN; CALL IDFDEALLOCATE(BTM,SIZE(BTM)); DEALLOCATE(BTM); ENDIF
 IF(ALLOCATED(IBND))DEALLOCATE(IBND)
 
 END SUBROUTINE SOLID_GEOTOP_DEALLOCATE

 !###======================================================================
 LOGICAL FUNCTION SOLID_NEWMASKS(ZOFFSET)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),INTENT(IN) :: ZOFFSET
 INTEGER :: I,IROW,ICOL,IRAT,IRAT1
 REAL(KIND=DP_KIND) :: T,B
 INTEGER,ALLOCATABLE,DIMENSION(:) :: NMSK
 CHARACTER(LEN=256) :: SOLNAME,STRING

 SOLID_NEWMASKS=.FALSE.

 !## solidtool
 IF(IBATCH.EQ.0)THEN
  !## refresh memory
  CALL SOLID_DEALLOCATE()
  SOLNAME=GETSOLNAME()
  !## read all idf-files in to memory
  IF(.NOT.SOLIDOPENSOL('R',SOLNAME))RETURN
  OUTPUTFOLDER=SOLNAME(:INDEX(SOLNAME,'\',.TRUE.)-1)
 ENDIF

 !## read top/bottom idf files, determine mdlidf()
 IF(.NOT.SOLID_READIDF())RETURN
 
 NMASK=NTBSOL
 IF(ALLOCATED(MASK))THEN
  DO I=1,SIZE(MASK); CALL IDFDEALLOCATEX(MASK(I)%IDF); ENDDO; DEALLOCATE(MASK)
 ENDIF
 ALLOCATE(MASK(NMASK),NMSK(NMASK))
 DO I=1,SIZE(MASK)
  CALL IDFNULLIFY(MASK(I)%IDF)
  !## copy mdlidf setting to masks
  IF(IBATCH.EQ.0)CALL IDFCOPY(SOLIDF(I),MASK(I)%IDF)
  IF(IBATCH.EQ.1)CALL IDFCOPY(MDLIDF,MASK(I)%IDF)
  NMSK(I)=0
 ENDDO

 !## read/scale/cut idf files for mdlidf() extent
 DO I=1,SIZE(SOLIDF)
  IF(IBATCH.EQ.0)THEN
   IF(.NOT.IDFREAD(MASK(I)%IDF,SOLIDF(I)%FNAME,1))RETURN !## scale mean
  ENDIF
  IF(IBATCH.EQ.1)THEN
   WRITE(*,'(A)') 'Reading '//TRIM(SOLIDF(I)%FNAME)//' ...'
   IF(.NOT.IDFREADSCALE(SOLIDF(I)%FNAME,MASK(I)%IDF,2,1,0.0D0,0))RETURN !## scale mean
  ENDIF
 ENDDO
 
 IF(IBATCH.EQ.1)THEN

  !## process all existing aquitards to get extinction
  IRAT=0; IRAT1=IRAT
  DO IROW=1,MDLIDF%NROW
   DO ICOL=1,MDLIDF%NCOL
    DO I=2,SIZE(SOLIDF)-2,2
     !## starting heads --- top/bot sequence
     T=MASK(I)  %IDF%X(ICOL,IROW)
     B=MASK(I+1)%IDF%X(ICOL,IROW)
     MASK(I  )%IDF%X(ICOL,IROW)=1
     MASK(I+1)%IDF%X(ICOL,IROW)=2 !1
     IF(T.NE.MASK(I)%IDF%NODATA.AND.B.NE.MASK(I+1)%IDF%NODATA)THEN
      IF(T-B.GT.ZOFFSET)THEN !## thickness aquitard
       MASK(I  )%IDF%X(ICOL,IROW)=-1
       MASK(I+1)%IDF%X(ICOL,IROW)=-1
       NMSK(I)  =NMSK(I)+1
       NMSK(I+1)=NMSK(I+1)+1
      ENDIF
     ENDIF
    END DO
   END DO
   IF(IBATCH.EQ.0)CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,MDLIDF%NROW,'Constructing Masks ... ')
   IF(IBATCH.EQ.1)WRITE(6,'(A,I10,A)') '+Constructing Masks',(IROW*100)/MDLIDF%NROW,'%'
  ENDDO
  !## get max extinction top/bot system (use layer 1 and layer nlay to distinghuish max. size system)
  DO IROW=1,MDLIDF%NROW
   DO ICOL=1,MDLIDF%NCOL
    !## starting heads --- top/bot sequence
    T=MASK(1)     %IDF%X(ICOL,IROW)
    B=MASK(NTBSOL)%IDF%X(ICOL,IROW)
    MASK(1)%IDF%X(ICOL,IROW)=1.0D0
    MASK(NTBSOL)%IDF%X(ICOL,IROW)=1.0D0
    IF(T.NE.MASK(1)%IDF%NODATA.AND.B.NE.MASK(NTBSOL)%IDF%NODATA)THEN
     IF(T-B.LE.0.0D0)THEN; DO I=1,SIZE(SOLIDF); MASK(I)%IDF%X(ICOL,IROW)=0; ENDDO; ENDIF
    ELSE
     DO I=1,SIZE(SOLIDF); MASK(I)%IDF%X(ICOL,IROW)=0; ENDDO
    ENDIF
   END DO
  END DO
 ENDIF 

 CALL UTL_CREATEDIR(TRIM(OUTPUTFOLDER)//'\MASK')
 
 DO I=1,SIZE(MASK)
  MASK(I)%FNAME=TRIM(OUTPUTFOLDER)//'\MASK\MASK_L'//TRIM(ITOS(I))//'.IDF'
  MASK(I)%ALIAS='MASK\MASK_L'//TRIM(ITOS(I))//'.IDF'
  IF(.NOT.IDFWRITE(MASK(I)%IDF,MASK(I)%FNAME,1))THEN; ENDIF
 END DO

 IF(IBATCH.EQ.0)THEN
  CALL UTL_MESSAGEHANDLE(1)
  CALL WDIALOGSELECT(ID_DSOLIDTAB3)
  CALL WDIALOGCLEARFIELD(IDF_MENU1)
  IF(NMASK.GT.0)CALL WDIALOGPUTMENU(IDF_MENU1,MASK%ALIAS,NMASK,1)
 ENDIF

 IF(IBATCH.EQ.1)THEN
  STRING='Summary of masks (number of cell with thickness of aquitards):'//CHAR(13)
  WRITE(*,'(A)') TRIM(STRING)
  DO I=2,SIZE(NMSK)-1,2
   IF(IBATCH.EQ.0)THEN
    STRING=TRIM(STRING)//'mask '//TRIM(ITOS(I))//'+'//TRIM(ITOS(I+1))//': '//TRIM(ITOS(NMSK(I)))//' points'//CHAR(13)
   ELSEIF(IBATCH.EQ.1)THEN
    STRING='mask '//TRIM(ITOS(I))//'+'//TRIM(ITOS(I+1))//': '//TRIM(ITOS(NMSK(I)))//' points'
    WRITE(*,'(A)') TRIM(STRING)
   ENDIF
  ENDDO
 ELSE
  CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'Successfully created masks','Information')
 ENDIF
 DEALLOCATE(NMSK)

! CALL IDFDEALLOCATE(SOLIDF,SIZE(SOLIDF)); DEALLOCATE(SOLIDF)
! DO I=1,SIZE(MASK); CALL IDFDEALLOCATEX(MASK(I)%IDF); CALL IDFDEALLOCATESX(MASK(I)%IDF); ENDDO; DEALLOCATE(MASK)
 
 CALL UTL_CLOSEUNITS()

 SOLID_NEWMASKS=.TRUE.

 END FUNCTION SOLID_NEWMASKS

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

 CALL WINDOWSELECT(0)
 IF(WMENUGETSTATE(ID_SOLIDS,2).EQ.1)THEN
  CALL SOLID_CLOSE()
  RETURN
 ENDIF

 CALL MAIN_UTL_INACTMODULE(ID_SOLIDS)
 !other module no closed, no approvement given
 IF(IDIAGERROR.EQ.1)RETURN

 CALL WMENUSETSTATE(ID_SOLIDS,2,1)

 IF(.NOT.IOSDIREXISTS(TRIM(PREFVAL(1))//'\SOLIDS'))CALL UTL_CREATEDIR(TRIM(PREFVAL(1))//'\SOLIDS')

 !## initialize sld()-pointer
 CALL SOLID_INITSLD(0)
 !## initialize spf()-pointer
 CALL SOLID_INITSPF(0)
 !## initialize polygon
 CALL POLYGON1INIT()

 CALL WDIALOGLOAD(ID_DSOLID,ID_DSOLID)

 CALL WDIALOGSELECT(ID_DSOLIDTAB1)
 CALL WDIALOGPUTIMAGE(ID_NEW,ID_ICONNEW,1)
 CALL WDIALOGPUTIMAGE(ID_PROFILE,ID_ICONPROFILE,1)
 CALL WDIALOGPUTIMAGE(ID_3D,ID_ICON3DTOOL,1)
 CALL WDIALOGPUTIMAGE(ID_CALCULATE,ID_ICONCALC,1)
 CALL WDIALOGPUTIMAGE(ID_INFO,ID_ICONINFO,1)
 CALL WDIALOGPUTIMAGE(ID_DELETE,ID_ICONDELETE,1)

 CALL WDIALOGSELECT(ID_DSOLIDTAB3)
 CALL WDIALOGPUTIMAGE(ID_NEW,ID_ICONNEW,1)
 CALL WDIALOGPUTIMAGE(ID_OPENIDF,ID_ICONOPENIDF,1)
 CALL WDIALOGPUTIMAGE(ID_INFO,ID_ICONPROPERTIES,1)
 CALL WDIALOGPUTIMAGE(ID_DELETE,ID_ICONDELETE,1)

 CALL WDIALOGSELECT(ID_DSOLID)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB2,0)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB3,0)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB4,0)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB5,0)

 CALL WDIALOGSELECT(ID_DSOLIDTAB1)
 CALL WDIALOGPUTCHECKBOX(IDF_CHECK1,0) 
 
 CALL SOLID_FIELDS('')

 !## fill polygonimages
 CALL POLYGON1IMAGES(ID_DSOLIDTAB2)
 CALL POLYGON1FIELDS(ID_DSOLIDTAB2)

 IBATCH=0
 
 CALL WDIALOGSELECT(ID_DSOLID)
 CALL UTL_DIALOGSHOW(-1,-1,0,2)

 END SUBROUTINE SOLID_INIT

 !###======================================================================
 SUBROUTINE SOLID_READSOL()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I 

 CALL WDIALOGSELECT(ID_DSOLID)
 CALL WDIALOGSELECT(ID_DSOLIDTAB2)
 CALL WDIALOGCLEARFIELD(IDF_MENU1)
 CALL WDIALOGSELECT(ID_DSOLIDTAB3)
 CALL WDIALOGCLEARFIELD(IDF_MENU1)
 CALL WDIALOGSELECT(ID_DSOLIDTAB4)
 CALL WDIALOGCLEARFIELD(IDF_GRID1)
 CALL WDIALOGSELECT(ID_DSOLIDTAB5)
 CALL WDIALOGCLEARFIELD(IDF_GRID1)
 I=0
 IF(SOLIDOPENSOL('R',GETSOLNAME()))THEN
  IF(SHP%NPOL.GT.0)THEN
   CALL WDIALOGSELECT(ID_DSOLIDTAB2)
   CALL WDIALOGPUTMENU(IDF_MENU1,SHP%POL%PNAME,SHP%NPOL,SHP%POL%IACT)
  ENDIF
  IF(NMASK.GT.0)THEN
   CALL WDIALOGSELECT(ID_DSOLIDTAB3)
   CALL WDIALOGPUTMENU(IDF_MENU1,MASK%ALIAS,NMASK,1)
  ENDIF
  IF(NGENSOL.GT.0)THEN
   CALL WDIALOGSELECT(ID_DSOLIDTAB4)
   CALL WGRIDPUTINTEGER(IDF_GRID1,1,GENSOL%ILAY,NGENSOL)
   CALL WGRIDPUTSTRING(IDF_GRID1,2,GENSOL%FNAME,NGENSOL)
  ENDIF
  IF(NIPFP.GT.0)THEN
   CALL WDIALOGSELECT(ID_DSOLIDTAB5)
   CALL WGRIDPUTINTEGER(IDF_GRID1,1,IPFP%ILAY,NIPFP)
   CALL WGRIDPUTSTRING(IDF_GRID1,2,IPFP%FNAME,NIPFP)
  ENDIF
  !## fill iMOD-Manager with appropriate idf files
  I=1
 ENDIF
 
 CALL WDIALOGSELECT(ID_DSOLIDTAB1)
 CALL WDIALOGFIELDSTATE(ID_PROFILE,I)
 CALL WDIALOGFIELDSTATE(ID_3D,I)
 CALL WDIALOGFIELDSTATE(ID_CALCULATE,I)
 CALL WDIALOGSELECT(ID_DSOLID)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB2,I)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB3,I)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB4,I)
 CALL WDIALOGTABSTATE(IDF_TAB1,ID_DSOLIDTAB5,I)

 END SUBROUTINE SOLID_READSOL

 !###======================================================================
 SUBROUTINE SOLID_FIELDS(FNAME)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: FNAME
 INTEGER :: N

 CALL WDIALOGSELECT(ID_DSOLIDTAB1)
 CALL UTL_IMODFILLMENU(IDF_MENU1,TRIM(PREFVAL(1))//'\SOLIDS','*','D',N,0,0,SETNAME=FNAME)
 IF(N.EQ.0)THEN
  CALL WDIALOGFIELDSTATE(ID_PROFILE,0)
  CALL WDIALOGFIELDSTATE(ID_CALCULATE,0)
  CALL WDIALOGFIELDSTATE(ID_INFO,0)
  CALL WDIALOGFIELDSTATE(ID_3D,0)
  CALL WDIALOGFIELDSTATE(ID_DELETE,0)
 ELSE
  CALL WDIALOGFIELDSTATE(ID_INFO,1)
  CALL WDIALOGFIELDSTATE(ID_PROFILE,1)
  CALL WDIALOGFIELDSTATE(ID_CALCULATE,1)
  CALL WDIALOGFIELDSTATE(ID_3D,1)
  CALL WDIALOGFIELDSTATE(ID_DELETE,1)
 ENDIF

 END SUBROUTINE SOLID_FIELDS

 !###======================================================================
 SUBROUTINE SOLID_CLEANMANAGER()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,L

 IF(.NOT.ALLOCATED(SLD))RETURN

 !## nothing selected
 MP%ISEL=.FALSE.
 !## select old occurences in list
 DO I=1,SIZE(SLD)
  DO L=1,MXMPLOT
   IF(TRIM(UTL_CAP(MP(L)%IDFNAME,'U')).EQ.UTL_CAP(TRIM(PREFVAL(1))//'\TMP\'//TRIM(SLD(I)%SNAME)//'.MDF','U'))THEN
    MP(L)%ISEL=.TRUE.
   ENDIF
  ENDDO
 ENDDO

 !## remove existing idf from manager (to be able to put them in the right order)
 CALL MANAGER_UTL_DELETE(IQ=0)

 END SUBROUTINE SOLID_CLEANMANAGER

 !###======================================================================
 LOGICAL FUNCTION SOLID_FILLMANAGER()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J,K,L,II
 CHARACTER(LEN=256) :: LINE

 SOLID_FILLMANAGER=.TRUE.

 IF(.NOT.ALLOCATED(SLD))RETURN

 SOLID_FILLMANAGER=.FALSE.

 CALL SOLID_CLEANMANAGER()

 J=0; DO I=1,SIZE(SLD); J=J+SLD(I)%NINT; ENDDO
 K=0; DO I=1,SIZE(MP); IF(MP(I)%IACT)K=K+1; ENDDO
 IF(J+K.GT.MXMPLOT)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'iMOD wants to store '//TRIM(ITOS(J+K))//' files.'//CHAR(13)// &
    'However, iMOD can store maximal '//TRIM(ITOS(MXMPLOT))//' files'//CHAR(13)//'in the iMOD Manager!','Error')
  RETURN
 ENDIF
 !## fill in manager with idf's
 DO I=1,SIZE(SLD)
  DO J=1,SLD(I)%NINT
   !## add top
   CALL MANAGER_UTL_ADDFILE(SLD(I)%INTNAME(J))
  END DO
 ENDDO
 !## select appropriate occurences in list and group them
 DO I=1,SIZE(SLD)
  !## nothing selected
  MP%ISEL=.FALSE.
  DO J=1,SLD(I)%NINT
   DO L=1,MXMPLOT
    IF(MP(L)%IACT)THEN
     IF(TRIM(UTL_CAP(MP(L)%IDFNAME,'U')).EQ.TRIM(UTL_CAP(SLD(I)%INTNAME(J),'U')))THEN
      MP(L)%ISEL=.TRUE.
      MP(L)%SCOLOR=SLD(I)%INTCLR(J)  !## color number for serie-plotting
      II=INDEX(SLD(I)%INTNAME(J),'\',.TRUE.)+1 !; LINE=SLD(I)%INTNAME(J)(:II); II=INDEX(LINE,'\',.TRUE.)+1
      LINE=SLD(I)%INTNAME(J)(II:)
      MP(L)%ALIAS=TRIM(LINE)
      CALL UTL_READARRAY((/1,0,0,1,0,0,0/),7,MP(L)%PRFTYPE) !## solid
      EXIT
     ENDIF
    ENDIF
   ENDDO
  ENDDO

  !## group them into mdf-file
  LINE=TRIM(PREFVAL(1))//'\TMP\'//TRIM(SLD(I)%SNAME)//'.MDF'
  IF(.NOT.MANAGER_UTL_GROUP(LINE))RETURN
  CALL MANAGER_UTL_ADDFILE(LINE)
  IF(I.EQ.1)THEN
   !## zoom full extent
   CALL IDFZOOM(ID_ZOOMFULLMAP,(MPW%XMAX+MPW%XMIN)/2.0,(MPW%YMAX+MPW%YMIN)/2.0,0)
   CALL IDFPLOTFAST(1)
  ENDIF
 ENDDO

 IF(NMASK.GT.0)THEN
  !## fill in manager with masks (why?) and group them
  DO I=1,NMASK; CALL MANAGER_UTL_ADDFILE(MASK(I)%FNAME); ENDDO
  !call idfplotfast()       
  !## nothing selected
  MP%ISEL=.FALSE.
  !## select appropriate occurences in list and group them
  DO I=1,SIZE(MASK)
   DO L=1,MXMPLOT
    IF(TRIM(UTL_CAP(MP(L)%IDFNAME,'U')).EQ.TRIM(UTL_CAP(MASK(I)%FNAME,'U')))THEN
     MP(L)%ISEL=.TRUE.
    ENDIF
   ENDDO
  ENDDO
  LINE=TRIM(PREFVAL(1))//'\TMP\MASK.MDF'
  !## group them into mdf-file
  IF(.NOT.MANAGER_UTL_GROUP(LINE))RETURN
  CALL MANAGER_UTL_ADDFILE(LINE)
 ENDIF

 !## select all group_solid{i}.mdf
 DO I=1,SIZE(SLD)
  DO L=1,MXMPLOT
   IF(TRIM(UTL_CAP(MP(L)%IDFNAME,'U')).EQ.TRIM(UTL_CAP(TRIM(PREFVAL(1))// &
     '\TMP\'//TRIM(SLD(I)%SNAME)//'.MDF','U')))THEN
    MP(L)%ISEL=.TRUE.
   ENDIF
  ENDDO
 ENDDO

 CALL MANAGER_UTL_FILL()
 CALL MANAGER_UTL_UPDATE()

 SOLID_FILLMANAGER=.TRUE.

 END FUNCTION SOLID_FILLMANAGER

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

 ISOLID=0

 !## read entire sol, incl. spf-files/polygons
 IF(SOLIDOPENSOL('R',GETSOLNAME()))THEN
  CALL UTL_HIDESHOWDIALOG(ID_DSOLID,0)
  IF(SOLID_PROFILEINIT())THEN
   ISOLID=1
   CALL PROFILE_INIT()
   CALL SOLID_PROFILECLOSE()
  ENDIF
  CALL UTL_HIDESHOWDIALOG(ID_DSOLID,2)
 ENDIF

 CALL IDFPLOT(1)

 END SUBROUTINE SOLID_PROFILE

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

 !## read entire sol, incl. spf-files
 IF(SOLIDOPENSOL('R',GETSOLNAME()))THEN
  CALL UTL_HIDESHOWDIALOG(ID_DSOLID,0)
  !## remove 3d tool prior to starting
  IF(WMENUGETSTATE(ID_3DTOOL,2).EQ.1)CALL IMOD3D_CLOSE()
  ISOLID=1
  CALL IMOD3D_INIT(1,0)
  !## save solid as something might has been changed
  
!  CALL SOLID_PROFILEINIT()
!  CALL PROFILE_INIT()
!  CALL SOLID_PROFILECLOSE()
  CALL UTL_HIDESHOWDIALOG(ID_DSOLID,2)
 ENDIF

 END SUBROUTINE SOLID_3D

 !###======================================================================
 LOGICAL FUNCTION SOLID_NEW()
 !###======================================================================
 IMPLICIT NONE
 TYPE(WIN_MESSAGE) :: MESSAGE
 INTEGER :: ITYPE,I,DID

 DID=WINFODIALOG(CURRENTDIALOG)

 SOLID_NEW=.FALSE.

 CALL WDIALOGLOAD(ID_DCREATESOLID,ID_DCREATESOLID)
 CALL WDIALOGTITLE('Create New Solid')
 CALL WDIALOGPUTSTRING(IDF_GROUP2,'Select the files to be included')
 CALL WDIALOGPUTSTRING(IDF_LABEL1,'Give the name for the Solid')

 CALL WDIALOGPUTDOUBLE(IDF_REAL1,MPW%XMIN,'(F15.3)')
 CALL WDIALOGPUTDOUBLE(IDF_REAL2,MPW%YMIN,'(F15.3)')
 CALL WDIALOGPUTDOUBLE(IDF_REAL3,MPW%XMAX,'(F15.3)')
 CALL WDIALOGPUTDOUBLE(IDF_REAL4,MPW%YMAX,'(F15.3)')

 !## free all resource
 CALL SOLID_DEALLOCATE()

 CALL UTL_MESSAGEHANDLE(1)

 !## idf(s) to be selected
 ALLOCATE(SOLID_ILIST(MXMPLOT))
 SOLID_ILIST=0

 DO I=1,MXMPLOT; IF(MP(I)%IPLOT.EQ.1)SOLID_ILIST(I)=1; END DO

 !## nothing selected, outgrey option 1
 IF(SUM(SOLID_ILIST).EQ.0)THEN
  !## else deactivate scratch modelling as sum(SOLID_ILIST)=0
  CALL WDIALOGFIELDSTATE(IDF_RADIO3,2)
  CALL WDIALOGPUTRADIOBUTTON(IDF_RADIO4)
 ENDIF
 CALL WDIALOGPUTMENU(IDF_MENU3,MP%ALIAS,MPW%NACT,SOLID_ILIST)
 CALL WDIALOGPUTMENU(IDF_MENU4,MP%ALIAS,MPW%NACT,1)
 CALL WDIALOGPUTMENU(IDF_MENU5,MP%ALIAS,MPW%NACT,1)
 CALL WDIALOGPUTDOUBLE(IDF_REAL5,100.0D0,'(F15.3)')
 
 CALL SOLID_NEW_FIELDS()
 CALL UTL_DIALOGSHOW(-1,-1,0,3)
 DO

  CALL WMESSAGE(ITYPE,MESSAGE)
  SELECT CASE(ITYPE)

   CASE(FIELDCHANGED)
    SELECT CASE (MESSAGE%VALUE2)
     CASE (IDF_RADIO3,IDF_RADIO4,IDF_CHECK1)
      CALL SOLID_NEW_FIELDS()
    END SELECT

   CASE(PUSHBUTTON)
    SELECT CASE (MESSAGE%VALUE1)
     CASE (IDOK)
      CALL WDIALOGSELECT(ID_DCREATESOLID)
      CALL WDIALOGGETSTRING(IDF_STRING1,SOLDIR)
      IF(LEN_TRIM(SOLDIR).EQ.0)THEN
       CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You should give a name for the Solid!','Warning')
      ELSE
       IF(SOLID_CREATE())THEN
        SOLID_NEW=.TRUE.
        EXIT
       ENDIF
       CALL WMESSAGEBOX(YESNO,EXCLAMATIONICON,COMMONNO,'Solid not created, process terminated.'//CHAR(13)// &
        'Remove progress upto now?','Error')
       IF(WINFODIALOG(4).EQ.1)THEN
        IF(.NOT.UTL_DEL1TREE(TRIM(SOLDIR)))THEN
!## error occured
        ENDIF
       ENDIF
       !## close all files
       CALL UTL_CLOSEUNITS()
       CALL UTL_MESSAGEHANDLE(1)
       CALL WDIALOGSELECT(ID_DCREATESOLID)
      ENDIF
     CASE (IDCANCEL)
      EXIT
     CASE (IDHELP)
       CALL UTL_GETHELP('5.4.1','TMO.ST.Create')
    END SELECT
  END SELECT

 ENDDO

 CALL SOLID_DEALLOCATE()

 CALL WDIALOGSELECT(ID_DCREATESOLID); CALL WDIALOGUNLOAD()
 CALL WDIALOGSELECT(DID)

 END FUNCTION SOLID_NEW

 !###======================================================================
 SUBROUTINE SOLID_NEW_FIELDS()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I,J,K

 !## select idfs or enter value
 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO3,I)
 IF(I.EQ.1)THEN; J=1; K=2; ELSE; J=2; K=1; ENDIF

 CALL WDIALOGFIELDSTATE(IDF_MENU3,J)
 CALL WDIALOGFIELDSTATE(IDF_LABEL11,K)
 CALL WDIALOGFIELDSTATE(IDF_LABEL12,K)
 CALL WDIALOGFIELDSTATE(IDF_MENU4,K)
 CALL WDIALOGFIELDSTATE(IDF_MENU5,K)

 !## coordinates
 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,I)
 CALL WDIALOGFIELDSTATE(IDF_REAL1,I)
 CALL WDIALOGFIELDSTATE(IDF_REAL2,I)
 CALL WDIALOGFIELDSTATE(IDF_REAL3,I)
 CALL WDIALOGFIELDSTATE(IDF_REAL4,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL5,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL6,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL7,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL8,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL9,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL10,I)

 !## cellsize
 I=0; IF(K.EQ.1)I=1
 CALL WDIALOGFIELDSTATE(IDF_REAL5,I)
 CALL WDIALOGFIELDSTATE(IDF_LABEL13,I)

 END SUBROUTINE SOLID_NEW_FIELDS

 !###======================================================================
 LOGICAL FUNCTION SOLID_CREATE()
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=50) :: FNAME
 INTEGER :: I,J,ICLIP,ICLR,ITOPBOT
 REAL(KIND=DP_KIND) :: XMIN,YMIN,XMAX,YMAX,CS
 LOGICAL :: LEX,LEQUAL

 SOLID_CREATE=.FALSE.

 CALL WDIALOGGETRADIOBUTTON(IDF_RADIO3,ITOPBOT)
 IF(ITOPBOT.EQ.1)THEN
  CALL WDIALOGGETMENU(IDF_MENU3,SOLID_ILIST)
  !## check whether all IDF files
  DO I=1,MXMPLOT
   IF(SOLID_ILIST(I).EQ.1)THEN
    IF(MP(I)%IPLOT.NE.1)THEN
     CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You need to select IDF files only!','Warning')
     RETURN
    ENDIF
   ENDIF
  ENDDO
  IF(MOD(SUM(SOLID_ILIST),2).NE.0)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You should give multiplies of two IDFs!','Warning')
   RETURN
  ENDIF
 ELSE
  SOLID_ILIST(1)=1 !## top idf/value
  SOLID_ILIST(2)=1 !## bot idf/value
 ENDIF

 FNAME=TRIM(SOLDIR)//'.SOL'
 INQUIRE(FILE=TRIM(PREFVAL(1))//'\SOLIDS\'//TRIM(SOLDIR)//'\'//TRIM(FNAME),EXIST=LEX)
 IF(LEX)THEN
  CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'Are you sure to overwrite the file:'//CHAR(13)// &
     TRIM(FNAME)//' in the folder:'//CHAR(13)//TRIM(PREFVAL(1))//'\SOLIDS\'//TRIM(SOLDIR)//CHAR(13)// &
     'Press YES to overwrite it and continue'//CHAR(13)//'Press NO to stop the process.','Question')
  IF(WINFODIALOG(4).NE.1)RETURN
  IF(UTL_DEL1TREE(TRIM(PREFVAL(1))//'\SOLIDS\'//TRIM(SOLDIR)))THEN
   INQUIRE(FILE=TRIM(PREFVAL(1))//'\SOLIDS\'//TRIM(SOLDIR)//'\'//TRIM(FNAME),EXIST=LEX)
   !## stil exists, probably terminated by user
   IF(LEX)RETURN
  ENDIF
 ENDIF

 SOLDIR=TRIM(PREFVAL(1))//'\SOLIDS\'//TRIM(SOLDIR)
 CALL UTL_CREATEDIR(SOLDIR)

 CALL WDIALOGGETCHECKBOX(IDF_CHECK1,ICLIP)
 CALL WDIALOGGETDOUBLE(IDF_REAL1,XMIN); CALL WDIALOGGETDOUBLE(IDF_REAL2,YMIN)
 CALL WDIALOGGETDOUBLE(IDF_REAL3,XMAX); CALL WDIALOGGETDOUBLE(IDF_REAL4,YMAX)
 CALL WDIALOGGETDOUBLE(IDF_REAL5,CS)

 IF(ALLOCATED(SOLIDF))THEN
  CALL IDFDEALLOCATE(SOLIDF,SIZE(SOLIDF)); DEALLOCATE(SOLIDF)
 ENDIF
 ALLOCATE(SOLIDF(2)); DO I=1,SIZE(SOLIDF); CALL IDFNULLIFY(SOLIDF(I)); ENDDO

 CALL WINDOWSELECT(0)

 !## ask to create zones between top and bottom
 IF(ITOPBOT.EQ.2)THEN 
  IF(.NOT.SOLID_SPLIT(ICLIP))RETURN
 ELSE
  !## store solid-fnames()
  CALL SOLID_INITSLD(1)
  SLD(1)%NINT=SUM(SOLID_ILIST) 
  CALL SOLID_INITSLDPOINTER(1,SLD(1)%NINT)

  LEQUAL=.FALSE.

  J=0
  ICLR=0
   
  DO I=1,MXMPLOT
   IF(SOLID_ILIST(I).EQ.1)THEN
    J=J+1
     
    IF(.NOT.LEQUAL)THEN
     !## make sure all idf files are equal to eachother
     IF(.NOT.IDFREAD(SOLIDF(2),MP(I)%IDFNAME,0))RETURN
     LEQUAL=.TRUE.
    ENDIF

    FNAME=MP(I)%IDFNAME(INDEX(MP(I)%IDFNAME,'\',.TRUE.)+1:)
    IF(ICLIP.EQ.0)THEN
     CALL WINDOWSELECT(0)
     CALL WINDOWOUTSTATUSBAR(4,'Copying '//TRIM(MP(I)%IDFNAME)//'...')
     IF(.NOT.IDFREAD(SOLIDF(1),MP(I)%IDFNAME,1))RETURN
     IF(.NOT.IDFEQUAL(SOLIDF(1),SOLIDF(2),1))RETURN
     IF(.NOT.IDFWRITE(SOLIDF(1),TRIM(SOLDIR)//'\'//TRIM(FNAME),1))RETURN
    ELSE
     CALL WINDOWSELECT(0)
     CALL WINDOWOUTSTATUSBAR(4,'Clipping '//TRIM(MP(I)%IDFNAME)//'...')
     IF(.NOT.IDFREAD(SOLIDF(1),MP(I)%IDFNAME,0))RETURN
     IF(.NOT.IDFREADPART(SOLIDF(1),XMIN,YMIN,XMAX,YMAX))RETURN
     IF(.NOT.IDFWRITE(SOLIDF(1),TRIM(SOLDIR)//'\'//TRIM(FNAME),1))RETURN
     CALL IDFDEALLOCATEX(SOLIDF(1))
    ENDIF

    !## store filenames
    SLD(1)%INTNAME(J)=TRIM(SOLDIR)//'\'//TRIM(FNAME)
    ICLR=ICLR+1
    IF(ICLR.GT.SIZE(ICOLOR))ICLR=1
    SLD(1)%INTCLR(J)=ICOLOR(ICLR)
    SLD(1)%ICLC(J)=1

   ENDIF
  END DO
 ENDIF

 CALL WINDOWSELECT(0)
 CALL WINDOWOUTSTATUSBAR(4,'')

 !## no polygons
 SHP%NPOL=0
 NMASK=0
 !## try to write the solid
 FNAME=SOLDIR(INDEX(SOLDIR,'\',.TRUE.)+1:)
 IF(.NOT.SOLIDOPENSOL('W',TRIM(SOLDIR)//'\'//TRIM(FNAME)//'.SOL',IQ=0))RETURN

 !## free all resource
 CALL SOLID_DEALLOCATE()

 !## refill menu with solids
 CALL SOLID_FIELDS(FNAME)
 
 CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'iMOD created the solid folder'//CHAR(13)// &
  TRIM(SOLDIR),'Information')

 SOLID_CREATE=.TRUE.

 END FUNCTION SOLID_CREATE

 !###======================================================================
 LOGICAL FUNCTION SOLID_SPLIT(ICLIP)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ICLIP
 INTEGER :: N,I,IROW,ICOL,JROW,JCOL,IOS
 REAL(KIND=DP_KIND) :: X,Y,TOP,BOT,DZZ,XTOP,XBOT
 TYPE(WIN_MESSAGE) :: MESSAGE
 INTEGER :: ITYPE
 CHARACTER(LEN=256) :: CTOP,CBOT

 SOLID_SPLIT=.FALSE.
 
 CALL WDIALOGLOAD(ID_DLEGENDNOCLASSES,ID_DLEGENDNOCLASSES)
 CALL WDIALOGPUTSTRING(IDF_RADIO1,'Number of interfaces')
 CALL WDIALOGPUTSTRING(IDF_GROUP2,'')
 CALL WDIALOGPUTINTEGER(IDF_INTEGER1,5)
 CALL WDIALOGFIELDSTATE(IDF_CHECK2,3)
 CALL WDIALOGFIELDSTATE(IDF_RADIO2,3)
 CALL WDIALOGFIELDSTATE(IDF_REAL1,3)
 CALL WDIALOGFIELDSTATE(IDF_REAL2,3)
 CALL WDIALOGFIELDSTATE(IDF_REAL3,3)
 CALL WDIALOGFIELDSTATE(IDF_LABEL2,3)
 CALL WDIALOGFIELDSTATE(IDF_LABEL3,3)

 CALL WDIALOGTITLE('Define the number of interfaces')
 CALL UTL_DIALOGSHOW(-1,-1,0,3)
 DO
  CALL WMESSAGE(ITYPE,MESSAGE)
  SELECT CASE (ITYPE)
   CASE (PUSHBUTTON)
    SELECT CASE (MESSAGE%VALUE1)
     CASE (IDOK)
      CALL WDIALOGGETINTEGER(IDF_INTEGER1,N)
      EXIT
     CASE (IDCANCEL)
      EXIT
     CASE (IDHELP)
       CALL UTL_GETHELP('5.4.1','TMO.ST.Create')
    END SELECT
  END SELECT
 ENDDO

 CALL WDIALOGUNLOAD(); IF(MESSAGE%VALUE1.EQ.IDCANCEL)RETURN

 CALL UTL_MESSAGEHANDLE(0)

 !## store solid-fnames()
 CALL SOLID_INITSLD(1)
 SLD(1)%NINT=N
 SLD(1)%SNAME=TRIM(SOLDIR(INDEX(SOLDIR,'\',.TRUE.)+1:))
 CALL SOLID_INITSLDPOINTER(1,SLD(1)%NINT); SLD(1)%ICLC=1; SLD(1)%ICHECK=1

 !## number of intervals
 NTBSOL=N 

 IF(ALLOCATED(SOLIDF))THEN; CALL IDFDEALLOCATE(SOLIDF,SIZE(SOLIDF)); DEALLOCATE(SOLIDF); ENDIF
 IF(ALLOCATED(TB))THEN; CALL IDFDEALLOCATE(TB,SIZE(TB)); DEALLOCATE(TB); ENDIF
 ALLOCATE(SOLIDF(NTBSOL),TB(2))
 DO I=1,SIZE(SOLIDF); CALL IDFNULLIFY(SOLIDF(I)); END DO
 DO I=1,SIZE(TB);     CALL IDFNULLIFY(TB(I));     END DO

 CALL WDIALOGSELECT(ID_DCREATESOLID)
 TB%IU=0
 !## top idf/value
 CALL WDIALOGGETMENU(IDF_MENU4,I,CTOP)

 !## existing filename
 IF(I.GT.0)THEN; CTOP=MP(I)%IDFNAME; IOS=1
 ELSE; READ(CTOP,*,IOSTAT=IOS) XTOP; ENDIF
 IF(IOS.NE.0)THEN
  !## check if it is an IDF file
  IF(INDEX(UTL_CAP(CTOP,'U'),'.IDF').EQ.0)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You need to specify an IDF for TOP-level.','Error')
   RETURN
  ENDIF
  IF(.NOT.IDFREAD(TB(1),CTOP,0))THEN; RETURN; ENDIF
 ENDIF
 !## bot idf/value
 CALL WDIALOGGETMENU(IDF_MENU5,I,CBOT)

 !## existing filename
 IF(I.GT.0)THEN; CBOT=MP(I)%IDFNAME; IOS=1
 ELSE; READ(CBOT,*,IOSTAT=IOS) XBOT; ENDIF
 IF(IOS.NE.0)THEN
  !## check if it is an IDF file
  IF(INDEX(UTL_CAP(CBOT,'U'),'.IDF').EQ.0)THEN
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You need to specify an IDF for BOT-level.','Error')
   RETURN
  ENDIF
  IF(.NOT.IDFREAD(TB(2),CBOT,0))THEN; RETURN; ENDIF
 ENDIF
 
 !## adjust first idf to iclip
 IF(ICLIP.EQ.1)THEN
  CALL WDIALOGGETDOUBLE(IDF_REAL1,SOLIDF(1)%XMIN)
  CALL WDIALOGGETDOUBLE(IDF_REAL2,SOLIDF(1)%YMIN)
  CALL WDIALOGGETDOUBLE(IDF_REAL3,SOLIDF(1)%XMAX)
  CALL WDIALOGGETDOUBLE(IDF_REAL4,SOLIDF(1)%YMAX)
  CALL WDIALOGGETDOUBLE(IDF_REAL5,SOLIDF(1)%DX)
  SOLIDF(1)%DY=SOLIDF(1)%DX
!  CALL IDFCOPY(TB(1),SOLIDF(1))
  CALL UTL_IDFSNAPTOGRID(SOLIDF(1)%XMIN,SOLIDF(1)%XMAX,SOLIDF(1)%YMIN,SOLIDF(1)%YMAX,SOLIDF(1)%DX,SOLIDF(1)%NCOL,SOLIDF(1)%NROW) 
 ELSEIF(ICLIP.EQ.0)THEN
  !## try first
  IF(TB(1)%IU.GT.0)THEN
   CALL IDFCOPY(TB(1),SOLIDF(1))
  ELSEIF(TB(2)%IU.GT.0)THEN
   CALL IDFCOPY(TB(2),SOLIDF(1))
  ELSE
   CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You need to specify at least one IDF if you do not enter clip coordinates.','Warning')
   RETURN
  ENDIF
 ENDIF

 !## write new files ... including top/bot given!
 DO I=1,NTBSOL

  !## first given IDF will be used to create all the others ... for now!
  IF(I.GT.1)CALL IDFCOPY(SOLIDF(1),SOLIDF(I))
  SOLIDF(I)%FNAME=TRIM(SOLDIR)//'\INT_L'//TRIM(ITOS(I))//'.IDF'
  !## open new idfs
  IF(.NOT.IDFOPEN(SOLIDF(I)%IU,SOLIDF(I)%FNAME,'WO',SOLIDF(I)%ITYPE,0,IQUESTION=1))THEN
   RETURN
  ENDIF

  !## store filenames
  SLD(1)%INTNAME(I)=SOLIDF(I)%FNAME
  SLD(1)%INTCLR(I)=ICOLOR(I)

 ENDDO

 SOLIDF%DMIN= 10.0D10
 SOLIDF%DMAX=-10.0D10

 DO IROW=1,SOLIDF(1)%NROW
  DO ICOL=1,SOLIDF(1)%NCOL

   !## get x/y coordinate
   CALL IDFGETLOC(SOLIDF(1),IROW,ICOL,X,Y)

   !## get top/bot information for current location
   IF(TB(1)%IU.GT.0)THEN
    CALL IDFIROWICOL(TB(1),JROW,JCOL,X,Y)
    TOP=IDFGETVAL(TB(1),JROW,JCOL)
   ELSE
    TOP=XTOP
   ENDIF
   IF(TB(2)%IU.GT.0)THEN
    CALL IDFIROWICOL(TB(2),JROW,JCOL,X,Y)
    BOT=IDFGETVAL(TB(2),JROW,JCOL)
   ELSE 
    BOT=XBOT
   ENDIF

   IF(TOP.NE.SOLIDF(1)%NODATA.AND. &
      BOT.NE.SOLIDF(2)%NODATA)THEN
    DZZ=(TOP-BOT)/REAL(NTBSOL-1) 
    !## check consistency
    IF(DZZ.LT.0.0D0)THEN
     DO I=1,SIZE(SOLIDF); CLOSE(SOLIDF(I)%IU); ENDDO
     CALL UTL_MESSAGEHANDLE(1)
     CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'iMOD cannot create a Solid whenever the values from the first'//CHAR(13)// &
      'given IDF are lower than the values in for the second given IDF.'//CHAR(13)// &
      'Check whether you should change the order of the selected IDF files.','Information')
     RETURN
    ENDIF
   ENDIF
   
   !## compute all intervals
   DO I=1,NTBSOL
    IF(TOP.NE.SOLIDF(1)%NODATA.AND. &
       BOT.NE.SOLIDF(2)%NODATA)THEN
     X=TOP-(DZZ*(I-1))
     SOLIDF(I)%DMIN=MIN(SOLIDF(I)%DMIN,X)
     SOLIDF(I)%DMAX=MAX(SOLIDF(I)%DMAX,X)
    ELSE
     X=SOLIDF(I)%NODATA
    ENDIF
    CALL IDFPUTVAL(SOLIDF(I),IROW,ICOL,X)
   END DO

  END DO
 END DO

 DO I=1,NTBSOL 
  IF(.NOT.IDFWRITEDIM(0,SOLIDF(I)))THEN
  ENDIF
  CLOSE(SOLIDF(I)%IU)
  CALL MANAGER_UTL_ADDFILE(IDFNAMEGIVEN=SOLIDF(I)%FNAME)
 ENDDO

 NSPF=0

 CALL UTL_MESSAGEHANDLE(1)

 SOLID_SPLIT=.TRUE.

 END FUNCTION SOLID_SPLIT

 !###======================================================================
 LOGICAL FUNCTION SOLID_CALC_KDC(IHYPO)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND) :: MINTHICKNESS=0.0D0
 INTEGER,INTENT(IN) :: IHYPO
 INTEGER :: I,J,K,ILAY,JLAY,IROW,ICOL,IU,IOS,NLAY,N 
 REAL(KIND=DP_KIND) :: TR,BR,Z1,Z2,F,KVAL,XTOP,XBOT,D
 CHARACTER(LEN=256) :: ROOT,FNAME,LINE
 CHARACTER(LEN=52) :: WC,CTYPE
 
 SOLID_CALC_KDC=.FALSE.

 !## try to read all idf's
 NLAY=SLD(1)%NINT/2; ALLOCATE(TOPIDF(NLAY),BOTIDF(NLAY)); J=1
 DO I=1,SLD(1)%NINT
  IF(IHYPO.EQ.0)THEN
   FNAME=SLD(1)%INTNAME(I)
  ELSE
   FNAME=TRIM(OUTPUTFOLDER)//'\'//TRIM(SLD(1)%INTNAME(I)(INDEX(SLD(1)%INTNAME(I),'\',.TRUE.)+1:))
  ENDIF
  WRITE(*,'(A)') 'Reading '//TRIM(FNAME)
  IF(MOD(I,2).NE.0)THEN; IF(.NOT.IDFREAD(TOPIDF(J),FNAME,1))RETURN; ENDIF
  IF(MOD(I,2).EQ.0)THEN; IF(.NOT.IDFREAD(BOTIDF(J),FNAME,1))RETURN; J=J+1; ENDIF
 ENDDO
 
 IF(.NOT.ASSOCIATED(REGISFILES_TOP))THEN
  !## get list of "regis"-files
  I=INDEX(REGISTOP,'\',.TRUE.); ROOT=REGISTOP(:I-1); WC=UTL_CAP(TRIM(REGISTOP(I+1:)),'U')
  IF(.NOT.UTL_DIRINFO_POINTER(ROOT,WC,REGISFILES,'F'))RETURN
  DO I=1,SIZE(REGISFILES); REGISFILES(I)=UTL_CAP(TRIM(ROOT)//'\'//TRIM(REGISFILES(I)),'U'); ENDDO
 ENDIF
  
 IF(.NOT.SOLID_CALC_KDC_INIT(NLAY))RETURN
 
 IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=TRIM(OUTPUTFOLDER)//'\solid_log.txt',STATUS='UNKNOWN',ACTION='WRITE',IOSTAT=IOS)
 IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'Cannot CREATE '//TRIM(OUTPUTFOLDER)//'\solid_log.txt'; RETURN; ENDIF

 IF(ASSOCIATED(REGISFILES))    N=SIZE(REGISFILES)
 IF(ASSOCIATED(REGISFILES_TOP))N=SIZE(REGISFILES_TOP)

 !## read/process
 DO I=1,N
  
  IF(ASSOCIATED(REGISFILES))THEN
   WRITE(*,'(2(I10,1X),F10.2)') I,SIZE(REGISFILES),REAL(I*100)/REAL(SIZE(REGISFILES))
   J=INDEX(REGISFILES(I),'\',.TRUE.)+1; K=INDEX(REGISFILES(I),TRIM(WC(2:)))-1
   CTYPE=REGISFILES(I)(J:K)
  ELSE
   WRITE(*,'(2(I10,1X),F10.2)') I,SIZE(REGISFILES_TOP),REAL(I*100)/REAL(SIZE(REGISFILES_TOP))
   J=INDEX(REGISFILES_TOP(I),'\',.TRUE.)+1; K=INDEX(REGISFILES_TOP(I),'.',.TRUE.)-1
   CTYPE=REGISFILES_TOP(I)(J:K)
  ENDIF

  WRITE(* ,'(A)') 'Reading:'; WRITE(IU,'(A)') 'Reading:'

  IF(ASSOCIATED(REGISFILES))THEN

   !## try top
   WRITE(* ,'(A)') '-'//TRIM(REGISFILES(I)); WRITE(IU,'(A)') '-'//TRIM(REGISFILES(I))
   IF(.NOT.IDFREADSCALE(REGISFILES(I),TB(1),2,0,0.0D0,0))RETURN !## scale mean
   TB(1)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(REGISFILES(I)(INDEX(REGISFILES(I),'\',.TRUE.):))

   !## try bot
   FNAME=UTL_SUBST(REGISBOT,'*',TRIM(CTYPE))
   WRITE(* ,'(A)') '-'//TRIM(FNAME); WRITE(IU,'(A)') '-'//TRIM(FNAME)
   IF(.NOT.IDFREADSCALE(FNAME,TB(2),2,0,0.0D0,0))RETURN !## scale mean
   TB(2)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(FNAME(INDEX(FNAME,'\',.TRUE.):))

   !## try kh
   FNAME=UTL_SUBST(REGISKHV,'*',TRIM(CTYPE))
   WRITE(* ,'(A)') '-'//TRIM(FNAME); WRITE(IU,'(A)') '-'//TRIM(FNAME)
   IF(.NOT.IDFREADSCALE(FNAME,KH(1),3,0,0.0D0,1))RETURN !## scale geometric
   KH(1)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(FNAME(INDEX(FNAME,'\',.TRUE.):))

   !## try kv
   FNAME=UTL_SUBST(REGISKVV,'*',TRIM(CTYPE))
   WRITE(*,'(A)') '-'//TRIM(FNAME); WRITE(IU,'(A)') '-'//TRIM(FNAME)
   IF(.NOT.IDFREADSCALE(FNAME,KV(1),3,0,0.0D0,1))RETURN !## scale geometric
   KV(1)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(FNAME(INDEX(FNAME,'\',.TRUE.):))
  
  ELSE

   !## try top
   WRITE(* ,'(A)') '- TOP: '//TRIM(REGISFILES_TOP(I)); WRITE(IU,'(A)') '-'//TRIM(REGISFILES_TOP(I))
   IF(.NOT.IDFREADSCALE(REGISFILES_TOP(I),TB(1),2,0,0.0D0,0))RETURN !## scale mean
   TB(1)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(REGISFILES_TOP(I)(INDEX(REGISFILES_TOP(I),'\',.TRUE.):))

   !## try bot
   WRITE(* ,'(A)') '- BOT: '//TRIM(REGISFILES_BOT(I)); WRITE(IU,'(A)') '-'//TRIM(REGISFILES_BOT(I))
   IF(.NOT.IDFREADSCALE(REGISFILES_BOT(I),TB(2),2,0,0.0D0,0))RETURN !## scale mean
   TB(2)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(REGISFILES_BOT(I)(INDEX(REGISFILES_BOT(I),'\',.TRUE.):))

   !## try khv
   WRITE(* ,'(A)') '- KHV: '//TRIM(REGISFILES_KHV(I)); WRITE(IU,'(A)') '-'//TRIM(REGISFILES_KHV(I))
   IF(.NOT.IDFREADSCALE(REGISFILES_KHV(I),KH(1),3,0,0.0D0,1))RETURN !## scale geometric
   KH(1)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(REGISFILES_KHV(I)(INDEX(REGISFILES_KHV(I),'\',.TRUE.):))

   !## try kvv
   WRITE(* ,'(A)') '- KVV: '//TRIM(REGISFILES_KVV(I)); WRITE(IU,'(A)') '-'//TRIM(REGISFILES_KVV(I))
   IF(.NOT.IDFREADSCALE(REGISFILES_KVV(I),KV(1),3,0,0.0D0,1))RETURN !## scale geometric
   KV(1)%FNAME=TRIM(OUTPUTFOLDER)//'\REGIS\'//TRIM(REGISFILES_KVV(I)(INDEX(REGISFILES_KVV(I),'\',.TRUE.):))

  ENDIF

  DO ILAY=1,SIZE(FFRAC)
   FFRAC(ILAY)%FNAME=TRIM(OUTPUTFOLDER)//'\FFRAC\'//TRIM(CTYPE)//'_l'//TRIM(ITOS(ILAY))//'.idf'; FFRAC(ILAY)%X=0.0D0; FFRAC(ILAY)%NODATA=0.0D0
  ENDDO
  DO ILAY=1,SIZE(CFRAC)
   CFRAC(ILAY)%FNAME=TRIM(OUTPUTFOLDER)//'\CFRAC\'//TRIM(CTYPE)//'_l'//TRIM(ITOS(ILAY))//'.idf'; CFRAC(ILAY)%X=0.0D0; CFRAC(ILAY)%NODATA=0.0D0
  ENDDO
  
  !## process data
  WRITE(*,'(A)') 'Process data ...'
  DO IROW=1,TOPIDF(1)%NROW; DO ICOL=1,TOPIDF(1)%NCOL
   
   TR=TB(1)%X(ICOL,IROW); BR=TB(2)%X(ICOL,IROW)   
   IF(TR.EQ.TB(1)%NODATA.OR.BR.EQ.TB(2)%NODATA)CYCLE
      
   !## compute fractions for aquifers
   DO ILAY=1,NLAY
    XTOP=TOPIDF(ILAY)%X(ICOL,IROW); XBOT=BOTIDF(ILAY)%X(ICOL,IROW)
    Z1=MIN(TR,XTOP); Z2=MAX(BR,XBOT)
    !## fraction in aquifer
    IF(Z1.GT.Z2.AND.XTOP.GT.XBOT)THEN
     F=(Z1-Z2)/(XTOP-XBOT)
     !## sum up the total fractions
     KDFRACIDF(ILAY)%X(ICOL,IROW)=KDFRACIDF(ILAY)%X(ICOL,IROW)+F

     !## sum horizontal transmissivity for each model layer
     KVAL=MAX(0.0D0,KH(1)%X(ICOL,IROW))
     KDHIDF(ILAY)%X(ICOL,IROW)=KDHIDF(ILAY)%X(ICOL,IROW)+(Z1-Z2)*KVAL

     !## sum vertical resistance for each model layer
     KVAL=MAX(0.0D0,KV(1)%X(ICOL,IROW))
     IF(KVAL.GT.0.0D0)KDVIDF(ILAY)%X(ICOL,IROW)=KDVIDF(ILAY)%X(ICOL,IROW)+(Z1-Z2)/KVAL

     !## sum up the total fractions
     FFRAC(ILAY)%X(ICOL,IROW)=F 

    ENDIF
   ENDDO

   !## compute fractions for aquitards, resistance between mid aquifer1 and aquifer2
   DO ILAY=1,NLAY-1
    XTOP=BOTIDF(ILAY  )%X(ICOL,IROW)
    XBOT=TOPIDF(ILAY+1)%X(ICOL,IROW)
    Z1=MIN(TR,XTOP); Z2=MAX(BR,XBOT)
    !## fraction in aquitards
    IF(Z1.GT.Z2)THEN
     F=(Z1-Z2)/(XTOP-XBOT)
     CFRACIDF(ILAY)%X(ICOL,IROW)=CFRACIDF(ILAY)%X(ICOL,IROW)+F

     !## found vertical permeability
     KVAL=MAX(0.0D0,KV(1)%X(ICOL,IROW)) !KDVIDF(ILAY)%X(ICOL,IROW)) !KV(1)%X(ICOL,IROW))
     !## sum up the total resistance
     IF(KVAL.GT.0.0D0)CIDF(ILAY)%X(ICOL,IROW)=CIDF(ILAY)%X(ICOL,IROW)+(Z1-Z2)/KVAL
     
     !## sum up the total fractions
     CFRAC(ILAY)%X(ICOL,IROW)=F 
    
    ENDIF
   ENDDO
   
  ENDDO; ENDDO
    
!  WRITE(* ,'(A)') 'Writing:'; WRITE(IU,'(A)') 'Writing:'
!  DO ILAY=1,NLAY
!   IF(SUM(FFRAC(ILAY)%X).GT.0.0D0)THEN
!    WRITE(* ,'(A)') '- '//TRIM(FFRAC(ILAY)%FNAME)
!    WRITE(IU,'(A)') '- '//TRIM(FFRAC(ILAY)%FNAME)
!    IF(.NOT.IDFWRITE(FFRAC(ILAY),FFRAC(ILAY)%FNAME,1))RETURN
!   ENDIF
!  ENDDO
!  DO ILAY=1,NLAY-1
!   IF(SUM(CFRAC(ILAY)%X).GT.0.0D0)THEN
!    WRITE(* ,'(A)') '- '//TRIM(CFRAC(ILAY)%FNAME)
!    WRITE(IU,'(A)') '- '//TRIM(CFRAC(ILAY)%FNAME)
!    IF(.NOT.IDFWRITE(CFRAC(ILAY),CFRAC(ILAY)%FNAME,1))RETURN
!   ENDIF
!  ENDDO
  
 ENDDO
 CLOSE(IU)
 
 !## check and remove small thicknesses for aquifers/(aquitards)
 DO IROW=1,TOPIDF(1)%NROW; DO ICOL=1,TOPIDF(1)%NCOL; DO ILAY=1,NLAY
  TR=TOPIDF(ILAY)%X(ICOL,IROW); BR=BOTIDF(ILAY)%X(ICOL,IROW)   
  IF(TR.EQ.TOPIDF(ILAY)%NODATA.OR.BR.EQ.BOTIDF(ILAY)%NODATA)CYCLE
  IF(TR-BR.LT.MINTHICKNESS)THEN  !T-BR=0.1 MINTHICKNESS=2
   D=(TR-BR)   !2-0.1=1.9
   KDHIDF(ILAY)%X(ICOL,IROW)=0.0D0
   KDVIDF(ILAY)%X(ICOL,IROW)=0.0D0
   !## shift all below
   DO JLAY=ILAY,NLAY
    BOTIDF(JLAY)%X(ICOL,IROW)=BOTIDF(JLAY)%X(ICOL,IROW)+D
    IF(JLAY+1.LT.NLAY)TOPIDF(JLAY+1)%X(ICOL,IROW)=TOPIDF(JLAY+1)%X(ICOL,IROW)+D
   ENDDO
  ENDIF
 ENDDO; ENDDO; ENDDO
 
 !## write transmissivities/vertical resistances
 DO I=1,SIZE(KDHIDF)
  LINE=TRIM(OUTPUTFOLDER)//'\mdl_kd_l'//TRIM(ITOS(I))//'.idf'
  IF(.NOT.IDFWRITE(KDHIDF(I),TRIM(LINE),1))RETURN
 ENDDO
 DO I=1,SIZE(CIDF)
  LINE=TRIM(OUTPUTFOLDER)//'\mdl_vc_l'//TRIM(ITOS(I))//'.idf'
  IF(.NOT.IDFWRITE(CIDF(I)  ,TRIM(LINE),1))RETURN
 ENDDO

 DO IROW=1,TOPIDF(1)%NROW; DO ICOL=1,TOPIDF(1)%NCOL; DO ILAY=1,NLAY
  TR=TOPIDF(ILAY)%X(ICOL,IROW); BR=BOTIDF(ILAY)%X(ICOL,IROW)   
  IF(TR.EQ.TOPIDF(ILAY)%NODATA.OR.BR.EQ.BOTIDF(ILAY)%NODATA)CYCLE
  !## compute vertical anisotropy
  IF(TR-BR.LE.0.0D0)THEN
   KDHIDF(ILAY)%X(ICOL,IROW)=0.0D0
   KDVIDF(ILAY)%X(ICOL,IROW)=1.0D0
  ELSE
   !## compute horizontal permeability
   KDHIDF(ILAY)%X(ICOL,IROW)=KDHIDF(ILAY)%X(ICOL,IROW)/(TR-BR)
   !## compute vertical permeability
   IF(KDVIDF(ILAY)%X(ICOL,IROW).GT.0.0D0)KDVIDF(ILAY)%X(ICOL,IROW)=(TR-BR)/KDVIDF(ILAY)%X(ICOL,IROW)
   IF(KDHIDF(ILAY)%X(ICOL,IROW).EQ.0.0D0)THEN
    KDVIDF(ILAY)%X(ICOL,IROW)=1.0D0
   ELSE
    KDVIDF(ILAY)%X(ICOL,IROW)=KDVIDF(ILAY)%X(ICOL,IROW)/KDHIDF(ILAY)%X(ICOL,IROW)
   ENDIF
  ENDIF
  !## compute vertical permeability
  IF(ILAY.LT.NLAY)THEN
   TR=BOTIDF(ILAY)%X(ICOL,IROW); BR=TOPIDF(ILAY+1)%X(ICOL,IROW)   
   IF(TR.EQ.TOPIDF(ILAY)%NODATA.OR.BR.EQ.BOTIDF(ILAY)%NODATA)CYCLE
   IF(CIDF(ILAY)%X(ICOL,IROW).LE.0.0D0)THEN
    CIDF(ILAY)%X(ICOL,IROW)= 0.0D0
   ELSE
    CIDF(ILAY)%X(ICOL,IROW)=(TR-BR)/CIDF(ILAY)%X(ICOL,IROW)
   ENDIF
  ENDIF
 ENDDO; ENDDO; ENDDO
 
 DO I=1,SIZE(KDHIDF);    IF(.NOT.IDFWRITE(KDHIDF(I),KDHIDF(I)%FNAME,1))RETURN; ENDDO
 DO I=1,SIZE(KDVIDF);    IF(.NOT.IDFWRITE(KDVIDF(I),KDVIDF(I)%FNAME,1))RETURN; ENDDO
 DO I=1,SIZE(CIDF);      IF(.NOT.IDFWRITE(CIDF(I),CIDF(I)%FNAME,1))RETURN; ENDDO
 DO I=1,SIZE(KDFRACIDF); IF(.NOT.IDFWRITE(KDFRACIDF(I),KDFRACIDF(I)%FNAME,1))RETURN; ENDDO
 DO I=1,SIZE(CFRACIDF);  IF(.NOT.IDFWRITE(CFRACIDF(I),CFRACIDF(I)%FNAME,1))RETURN; ENDDO
 !## overwrite top/bot for small modifications
 DO I=1,SIZE(TOPIDF);    IF(.NOT.IDFWRITE(TOPIDF(I),TOPIDF(I)%FNAME,1))RETURN; ENDDO
 DO I=1,SIZE(BOTIDF);    IF(.NOT.IDFWRITE(BOTIDF(I),BOTIDF(I)%FNAME,1))RETURN; ENDDO

 !## make boundary
 DO ILAY=1,NLAY; DO IROW=1,TOPIDF(1)%NROW; DO ICOL=1,TOPIDF(1)%NCOL
  IF(TOPIDF(ILAY)%X(ICOL,IROW).EQ.TOPIDF(ILAY)%NODATA)THEN
   TOPIDF(ILAY)%X(ICOL,IROW)=0.0D0
  ELSE
   TOPIDF(ILAY)%X(ICOL,IROW)=1.0D0
  ENDIF
 ENDDO; ENDDO; TOPIDF(ILAY)%FNAME=TRIM(OUTPUTFOLDER)//'\mdl_ib_l'//TRIM(ITOS(ILAY))//'.idf'; ENDDO
 DO I=1,SIZE(TOPIDF); IF(.NOT.IDFWRITE(TOPIDF(I),TOPIDF(I)%FNAME,1))RETURN; ENDDO
 
 SOLID_CALC_KDC=.TRUE.
 
 END FUNCTION SOLID_CALC_KDC

 !###======================================================================
 LOGICAL FUNCTION SOLID_CALC_KDC_INIT(NLAY)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NLAY
 INTEGER :: I
 
 SOLID_CALC_KDC_INIT=.FALSE.
 
 ALLOCATE(CIDF(NLAY-1),CFRACIDF(NLAY-1),KDHIDF(NLAY),KDVIDF(NLAY),KDFRACIDF(NLAY),TOP(NLAY), &
          BOT(NLAY),TB(2),KH(1),KV(1),FFRAC(NLAY),CFRAC(NLAY-1))
 !## nullify idf-objects
 DO I=1,SIZE(CIDF);      CALL IDFNULLIFY(CIDF(I)) ; ENDDO
 DO I=1,SIZE(KDHIDF);    CALL IDFNULLIFY(KDHIDF(I)); ENDDO 
 DO I=1,SIZE(KDVIDF);    CALL IDFNULLIFY(KDVIDF(I)); ENDDO 
 DO I=1,SIZE(CFRACIDF);  CALL IDFNULLIFY(CFRACIDF(I)) ; ENDDO
 DO I=1,SIZE(KDFRACIDF); CALL IDFNULLIFY(KDFRACIDF(I)); ENDDO 
 DO I=1,SIZE(KH);        CALL IDFNULLIFY(KH(I));    ENDDO 
 DO I=1,SIZE(KV);        CALL IDFNULLIFY(KV(I));    ENDDO 
 DO I=1,SIZE(TB);        CALL IDFNULLIFY(TB(I));    ENDDO
 DO I=1,SIZE(FFRAC);     CALL IDFNULLIFY(FFRAC(I));    ENDDO
 DO I=1,SIZE(CFRAC);     CALL IDFNULLIFY(CFRAC(I));    ENDDO
 !## copy settings0
 DO I=1,SIZE(CIDF);  CALL IDFCOPY(TOPIDF(1),CIDF(I)) ; ENDDO 
 DO I=1,SIZE(KDHIDF); CALL IDFCOPY(TOPIDF(1),KDHIDF(I)); ENDDO 
 DO I=1,SIZE(KDVIDF); CALL IDFCOPY(TOPIDF(1),KDVIDF(I)); ENDDO 
 DO I=1,SIZE(CFRACIDF);  CALL IDFCOPY(TOPIDF(1),CFRACIDF(I)) ; ENDDO 
 DO I=1,SIZE(KDFRACIDF); CALL IDFCOPY(TOPIDF(1),KDFRACIDF(I)); ENDDO 
 DO I=1,SIZE(TB)   ; CALL IDFCOPY(TOPIDF(1),TB(I))   ; ENDDO 
 DO I=1,SIZE(FFRAC)   ; CALL IDFCOPY(TOPIDF(1),FFRAC(I))   ; ENDDO 
 DO I=1,SIZE(CFRAC)   ; CALL IDFCOPY(TOPIDF(1),CFRAC(I))   ; ENDDO 
 CALL IDFCOPY(TOPIDF(1),KH(1)); CALL IDFCOPY(TOPIDF(1),KV(1))
 
 !## write header information kdidf / cidf
 DO I=1,SIZE(KDHIDF)
  KDHIDF(I)%X=0.0D0; KDHIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\mdl_khv_l'//TRIM(ITOS(I))//'.idf'
 ENDDO
 DO I=1,SIZE(KDHIDF)
  KDVIDF(I)%X=0.0D0; KDVIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\mdl_kva_l'//TRIM(ITOS(I))//'.idf'
 ENDDO
 DO I=1,SIZE(CIDF)
  CIDF(I)%X=0.0D0;   CIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\mdl_kvv_l'//TRIM(ITOS(I))//'.idf'
 ENDDO
 DO I=1,SIZE(KDFRACIDF)
  KDFRACIDF(I)%X=0.0D0; KDFRACIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\mdl_kdfrac_l'//TRIM(ITOS(I))//'.idf'
 ENDDO
 DO I=1,SIZE(CFRACIDF)
  CFRACIDF(I)%X=0.0D0; CFRACIDF(I)%FNAME=TRIM(OUTPUTFOLDER)//'\mdl_cfrac_l'//TRIM(ITOS(I))//'.idf'
 ENDDO
  
 SOLID_CALC_KDC_INIT=.TRUE.
 
 END FUNCTION SOLID_CALC_KDC_INIT

 !###======================================================================
 SUBROUTINE SOLID_CALC_KDC_DEALLOCATE()
 !###======================================================================
 IMPLICIT NONE
 
 !## deallocate
 IF(ALLOCATED(TOPIDF))THEN;    CALL IDFDEALLOCATE(TOPIDF,SIZE(TOPIDF));       DEALLOCATE(TOPIDF); ENDIF
 IF(ALLOCATED(BOTIDF))THEN;    CALL IDFDEALLOCATE(BOTIDF,SIZE(BOTIDF));       DEALLOCATE(BOTIDF); ENDIF
 IF(ALLOCATED(CIDF))THEN;      CALL IDFDEALLOCATE(CIDF ,SIZE(CIDF));          DEALLOCATE(CIDF);   ENDIF
 IF(ALLOCATED(KDHIDF))THEN;    CALL IDFDEALLOCATE(KDHIDF,SIZE(KDHIDF));       DEALLOCATE(KDHIDF);  ENDIF
 IF(ALLOCATED(KDVIDF))THEN;    CALL IDFDEALLOCATE(KDVIDF,SIZE(KDVIDF));       DEALLOCATE(KDVIDF);  ENDIF
 IF(ALLOCATED(CFRACIDF))THEN;  CALL IDFDEALLOCATE(CFRACIDF ,SIZE(CFRACIDF));  DEALLOCATE(CFRACIDF); ENDIF
 IF(ALLOCATED(KDFRACIDF))THEN; CALL IDFDEALLOCATE(KDFRACIDF,SIZE(KDFRACIDF)); DEALLOCATE(KDFRACIDF); ENDIF
 IF(ASSOCIATED(REGISFILES))DEALLOCATE(REGISFILES)
 IF(ALLOCATED(TOP))DEALLOCATE(TOP) 
 IF(ALLOCATED(BOT))DEALLOCATE(BOT) 
 IF(ALLOCATED(KH))THEN; CALL IDFDEALLOCATE(KH,SIZE(KH)); DEALLOCATE(KH); ENDIF
 IF(ALLOCATED(KV))THEN; CALL IDFDEALLOCATE(KV,SIZE(KV)); DEALLOCATE(KV); ENDIF
 IF(ALLOCATED(TB))THEN; CALL IDFDEALLOCATE(TB,SIZE(TB)); DEALLOCATE(TB); ENDIF 
 IF(ALLOCATED(FFRAC))THEN; CALL IDFDEALLOCATE(FFRAC,SIZE(FFRAC)); DEALLOCATE(FFRAC); ENDIF 
 IF(ALLOCATED(CFRAC))THEN; CALL IDFDEALLOCATE(CFRAC,SIZE(CFRAC)); DEALLOCATE(CFRAC); ENDIF 
 
 END SUBROUTINE SOLID_CALC_KDC_DEALLOCATE

END MODULE MOD_SOLID