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


      subroutine isginit()
! description:
! ------------------------------------------------------------------------------
! init isg
!

! declaration section
! ------------------------------------------------------------------------------
      use mod_glbvar, only: iss,cdate,delt,simcsize,simbox,lqd
      use mod_basvar, only: delr,delc,nrow,ncol
      use global, only: iunit
      use m_mf2005_iu, only: iumet

      implicit none

! program section
! ------------------------------------------------------------------------------

! check if MET package is activated
      if (IUNIT(IUMET).le.0) then
         write(*,*) 'ISG cannot be used without met-package.'
         call ustop(' ')
      end if

      call isginit1(nrow,ncol)
      call isginit2(iss,cdate,delt,simcsize,simbox,delr,delc,nrow,ncol,lqd)

! end of program
      return
      end subroutine

      subroutine isginit1(im_nrow, im_ncol)
      use global, only: ncol, nrow
      use mod_basvar, only: delr, delc
      use pksmpi_mod,only: nrproc, gncol, gnrow     
                           
      implicit none
      integer, intent(out) :: im_nrow, im_ncol
      
      if(nrproc.gt.1) then
         im_ncol = gncol
         im_nrow = gnrow
      else
         im_ncol = ncol
         im_nrow = nrow
      end if   

      if(.not.allocated(delr)) allocate(delr(0:im_ncol))
      if(.not.allocated(delc)) allocate(delc(0:im_nrow))

! end of program
      return
      end subroutine

      subroutine isginit2(im_iss,im_cdate,im_delt,im_simcsize,im_simbox,im_delr,im_delc,im_nrow,im_ncol,lqd)

      use global, only: ncol, nrow, delr, delc, iunit, issflg, nper
      use gwfmetmodule, only: coord_xll, coord_yll, coord_xur, coord_yur,&
                              gcoord_xll, gcoord_yll, gcoord_xur, gcoord_yur,&
                              time_cstring, iss, ieq, cdelr, cdelc
      use m_mf2005_main, only: timesteptime, kper
      use gwfbasmodule, only: delt
      use pksmpi_mod,only: nrproc, gncol, gnrow 

      implicit none

! arguments
      real, intent(inout) :: im_delt
      real(kind=8), intent(inout) :: im_simcsize
      integer, intent(inout) :: im_iss
      character(len=20), intent(inout) :: im_cdate
      real(KIND=8), dimension(4), intent(inout) :: im_simbox
      real(KIND=8), dimension(0:im_ncol), intent(inout) :: im_delr
      real(KIND=8), dimension(0:im_nrow), intent(inout) :: im_delc
      integer, intent(inout) :: im_ncol, im_nrow
      logical, intent(inout) :: lqd

! local variables
      double precision :: t
      integer :: i, icol, irow, date, hour, minute, second

! program section
! ------------------------------------------------------------------------------

! find out or model is steady state or not
      im_iss = 1
      if (.not.associated(iss)) then
         write(*,*) 'Error, initialization ISG'
         call ustop(' ')
      else
         im_iss = iss
      end if
      if (im_iss.eq.0) im_iss = 2

      im_delt = delt

! check if grid is uniform
      lqd = .true.
      if (.not.associated(ieq)) then
         write(*,*) 'Error, initialization ISG'
         call ustop(' ')
      else
         if (ieq.eq.1) lqd = .false.
      end if

      if (lqd) then
         im_simcsize = delr(1)
      else
       if(nrproc.gt.1)then
         write(*,*) 'Error: non-uniform grids not yet supported icm PKS.'
         call ustop(' ')
       endif
      end if

      if (.not.associated(cdelr).or..not.associated(cdelc)) then
         write(*,*) 'Error, initialization ISG'
         call ustop(' ')
      else
         if (nrproc.eq.1) then 
            im_delr(0:im_ncol) = cdelr(0:ncol)
            im_delc(0:im_nrow) = cdelc(0:nrow)
         else    
            im_delr(0) = gcoord_xll
            do icol = 1, im_ncol
               im_delr(icol) = gcoord_xll + real(icol)*im_simcsize
            end do
            im_delc(0) = gcoord_yur
            do irow = 1, im_nrow
               im_delc(irow) = gcoord_yur - real(irow)*im_simcsize
            end do
         end if      
      end if

      if (nrproc.eq.1) then 
         im_simbox(1) = coord_xll
         im_simbox(2) = coord_yll
         im_simbox(3) = coord_xur
         im_simbox(4) = coord_yur
      else
         im_simbox(1) = gcoord_xll
         im_simbox(2) = gcoord_yll
         im_simbox(3) = gcoord_xur
         im_simbox(4) = gcoord_yur
      end if          
      
      t = timesteptime
      call cfn_mjd2datehms(t,date,hour,minute,second)
      write(im_cdate,'(i8)') date
      im_cdate = adjustl(im_cdate)

! end of program
      return
      end subroutine

      subroutine isgfinalize()
      use mod_basvar, only: delr, delc
      implicit none

      if(allocated(delr)) deallocate(delr)
      if(allocated(delc)) deallocate(delc)

! end of program
      return
      end subroutine

!###====================================================================
SUBROUTINE PCK1RPISG(NISG,FNAME,ILAY)
!###====================================================================
USE IMOD_UTL, ONLY : IMOD_UTL_IDATETOJDATE,IMOD_UTL_FILENAME,IMOD_UTL_PRINTTEXT,IMOD_UTL_ITOS, &
    IMOD_UTL_OPENIDF,IMOD_UTL_OPENASC,IMOD_UTL_INTERSECT_EQUI,IMOD_UTL_INTERSECT_NONEQUI,IMOD_UTL_GETRCL, &
   XA,YA,ZA,FA,LN,CA,RA,INTERSECT_DEALLOCATE,UTL_GETRECORDLENGTH
USE MOD_GLBVAR, ONLY : LINE,ISS,CDATE,DELT,SIMCSIZE,SIMBOX,LQD
USE IMOD_UTL, ONLY : ICF,UTL_IDFIROWICOL
USE MOD_BASVAR, ONLY : DELR,DELC,NROW,NCOL
USE MOD_ISGVAR
USE IMOD_IDF_PAR
USE IMOD_IDF, ONLY : IDFNULLIFY,IDFALLOCATEX,IDFDEALLOCATEX,IDFDEALLOCATE
IMPLICIT NONE
REAL(KIND=8),PARAMETER :: MAXWIDTH=1000.0D0
INTEGER, INTENT(OUT) :: NISG
CHARACTER(LEN=*),INTENT(IN) :: FNAME
INTEGER,INTENT(IN) :: ILAY 
REAL(KIND=8),POINTER,DIMENSION(:,:) :: ISGLIST_BU
CHARACTER(LEN=32) :: CNAME
INTEGER,DIMENSION(8) :: ISGIU
INTEGER :: I,II,III,J,JJ,K,SDATE,EDATE,TTIME,MX,N,IROW,ICOL,ICRS,NCRS,IREC,ICROS,N2DIM,IC,IR,ISGTYPE
INTEGER :: MAXCROS,JCRS,NCLC,IREF,JREC,ICLC,NR,ISTW,NSTW,NS1,NS2,NN,NODE,ICSTYPE
INTEGER(KIND=2) :: ZM
INTEGER(KIND=1) :: ZC,ZP
REAL(KIND=8) ::  DXY,C,WP,BH,RWIDTH,WETPER,ISGLEN,MIND,D,FCT,H1,H2,Z,WL,INFF,ZCHK
REAL(KIND=8),ALLOCATABLE,DIMENSION(:,:) :: QSORT,RVAL
REAL(KIND=8),ALLOCATABLE,DIMENSION(:) :: DIST
REAL(KIND=8),ALLOCATABLE,DIMENSION(:) :: XNR
INTEGER,ALLOCATABLE,DIMENSION(:) :: IPOS
REAL(KIND=8) ::  NODATA=-9999.0D0
REAL(KIND=4) :: XS,YS,DS,BS,KS
LOGICAL :: LNODAT
REAL(KIND=8) ::  WATERDEPTH=0.10D0 !## minimal waterdepth
TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: IDF
TYPE(IDFOBJ) :: ICROSS,PCROSS
LOGICAL :: ISGATTRIBUTES_2DCROSS_READ
INTEGER :: IOS,NA,MA
REAL(KIND=8),PARAMETER :: NODATAIDF=0.0D0 !## do not change !!!

IF(ALLOCATED(QSORT))DEALLOCATE(QSORT)
call INTERSECT_DEALLOCATE()
MX=1000; ALLOCATE(XA(MX),YA(MX),FA(MX),LN(MX))

!## translate cdate in to julian date - for transient simulations only!
IF(ISS.EQ.2)THEN
 IF(LEN_TRIM(CDATE).EQ.8)THEN
  READ(CDATE,*,IOSTAT=IOS) SDATE
 ELSE
  !## make sure sdate is only 8 integer (yyyymmdd)
  READ(CDATE,'(I8)',IOSTAT=IOS) SDATE
 ENDIF
 IF(IOS.EQ.0)THEN
  SDATE=IMOD_UTL_IDATETOJDATE(SDATE)
  EDATE=SDATE+MAX(1,INT(DELT))
  TTIME=EDATE-SDATE 
 ELSE
  CALL IMOD_UTL_PRINTTEXT('iMODFLOW can not convert entered date ('//TRIM(CDATE)//') into an integer value.',2)
 ENDIF
ELSEIF(ISS.EQ.1)THEN
 TTIME=1
ENDIF
ALLOCATE(QSORT(TTIME,MAXITEMS),XNR(MAXITEMS))

!## initialize unitnumbers
ISGIU=0
I=INDEX(FNAME,'.',.TRUE.)-1
!## open *.ISG file
LINE=FNAME(:I)//'.ISG'
CALL IMOD_UTL_FILENAME(LINE)
CALL IMOD_UTL_OPENASC(ISGIU(1),LINE,'R')
!## open *.ISP file
LINE=FNAME(:I)//'.ISP'
CALL IMOD_UTL_FILENAME(LINE)
J=UTL_GETRECORDLENGTH(LINE) 
ISGTYPE=4; IF(J.EQ.16)ISGTYPE=8
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(2),LINE,'R',8)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(2),LINE,'R',16)
!## open *.ISD1 file
LINE=FNAME(:I)//'.ISD1'
CALL IMOD_UTL_FILENAME(LINE)
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(3),LINE,'R',44)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(3),LINE,'R',48)
!## open *.ISD2 file
LINE=FNAME(:I)//'.ISD2'
CALL IMOD_UTL_FILENAME(LINE)
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(4),LINE,'R',20)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(4),LINE,'R',44)
!## open *.ISC1 file
LINE=FNAME(:I)//'.ISC1'
CALL IMOD_UTL_FILENAME(LINE)
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(5),LINE,'R',44)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(5),LINE,'R',48)
!## open *.ISC2 file
LINE=FNAME(:I)//'.ISC2'
CALL IMOD_UTL_FILENAME(LINE)
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(6),LINE,'R',12)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(6),LINE,'R',24)
!## open *.IST1 file
LINE=FNAME(:I)//'.IST1'
CALL IMOD_UTL_FILENAME(LINE)
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(7),LINE,'R',44)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(7),LINE,'R',48)
!## open *.IST2 file
LINE=FNAME(:I)//'.IST2'
CALL IMOD_UTL_FILENAME(LINE)
IF(ISGTYPE.EQ.4)CALL IMOD_UTL_OPENIDF(ISGIU(8),LINE,'R',12)
IF(ISGTYPE.EQ.8)CALL IMOD_UTL_OPENIDF(ISGIU(8),LINE,'R',28)

READ(ISGIU(1),*) NISGH

NXY=0
NISG=0

ALLOCATE(DWP(1),CROS(1,1))

!## start with this
ALLOCATE(ISGDATA(NROW*NCOL,11))

!## number of two-dimensional cross-sections
N2DIM=0

DO I=1,NISGH
 READ(ISGIU(1),*) CISGH,ISEG,NSEG,ICLC,NCLC,ICRS,NCRS,ISTW,NSTW
 !## skip elements with le 1 segment
 IF(NSEG.LT.2)CYCLE; IF(NCLC.LT.2)CYCLE
 !## cannot process weirs the way iMOD does it
 IF(NSTW.GT.0)THEN
  CALL IMOD_UTL_PRINTTEXT('Segment '//TRIM(IMOD_UTL_ITOS(I))//' ['//TRIM(CISGH)//'] has weirs, effect of this is computed to nearest calc. node.',0)
 ENDIF

 !## increase memory for segments (including points from calculation-nodes!)
 IF(NSEG+NCLC+2*NSTW.GT.NXY)THEN
  NXY=NSEG+NCLC+2*NSTW
  IF(ALLOCATED(X))   DEALLOCATE(X);    IF(ALLOCATED(Y))   DEALLOCATE(Y)
  IF(ALLOCATED(RVAL))DEALLOCATE(RVAL); IF(ALLOCATED(DIST))DEALLOCATE(DIST)
  IF(ALLOCATED(IPOS))DEALLOCATE(IPOS)
  ALLOCATE(X(NXY),Y(NXY),DIST(NXY),IPOS(NXY),RVAL(MAXITEMS,0:NXY))
 ENDIF

 !## read segments
 IREC=ISEG-1
 DO J=1,NSEG
  IF(ISGTYPE.EQ.4)THEN
   READ(ISGIU(2),REC=IREC+J+ICF) XS,YS
   X(J)=XS; Y(J)=YS
  ELSE
   READ(ISGIU(2),REC=IREC+J+ICF) X(J),Y(J)
  ENDIF
 END DO

 !## evaluate whether line is within model areas, otherwise skip rest of procedure ...
 IF(MAXVAL(X(1:NSEG)).GT.SIMBOX(1).AND.MAXVAL(Y(1:NSEG)).GT.SIMBOX(2).AND. &
    MINVAL(X(1:NSEG)).LT.SIMBOX(3).AND.MINVAL(Y(1:NSEG)).LT.SIMBOX(4))THEN

  !## read cross-section information from *.isc
  IF(NCRS.GT.SIZE(DWP))THEN 
   IF(ALLOCATED(DWP))DEALLOCATE(DWP); ALLOCATE(DWP(NCRS))
  ENDIF

  !## read cross-sectional data
  MAXCROS=0
  IREC=ICRS-1
  DO J=1,NCRS
   READ(ISGIU(5),REC=IREC+J+ICF) DWP(J)%NCROS
   !## do not read in two-dimensional cross-sections
   MAXCROS=MAXCROS+MAX(0,DWP(J)%NCROS)
  END DO

  IF(MAXCROS.GT.SIZE(CROS,2).OR.NCRS.GT.SIZE(CROS,1))THEN !MAXCROS.GT.DIMCROS.OR.NCRS.GT.DIMDWP)THEN
   IF(MAXCROS.GT.0.AND.NCRS.GT.0)THEN
    IF(ALLOCATED(CROS))DEALLOCATE(CROS); ALLOCATE(CROS(NCRS,MAXCROS))
   ENDIF
  ENDIF

  IPOS=0
  !## include calculation nodes in between segments
  CALL PCK7RPISG(ICLC,ISGIU(3),NCLC,NXY,NSEG,X,Y,DIST,IPOS,+1,ISGTYPE)
  !## include structure nodes in between segments
  CALL PCK7RPISG(ISTW,ISGIU(7),NSTW,NXY,NSEG,X,Y,DIST,IPOS,-1,ISGTYPE)

  !## read/prepare all information for current date/segment where information is connected to
  J=0
  DO
   J=J+1
   IF(IPOS(J).GT.0)THEN
    READ(ISGIU(3),REC=IPOS(J)+ICF) NR,IREF
    CALL PCK2RPISG(SDATE,EDATE,TTIME,QSORT,XNR,ISGIU(4),MAXITEMS,RVAL(1,J),NR,IREF,I,1,NODATA,ISGTYPE)
   ENDIF
   IF(IPOS(J).LT.0)THEN
    READ(ISGIU(7),REC=ABS(IPOS(J))+ICF) NR,IREF
    CALL PCK2RPISG(SDATE,EDATE,TTIME,QSORT,XNR,ISGIU(8),MAXITEMS,RVAL(1,J),NR,IREF,I,-1,NODATA,ISGTYPE)
    J=J+1
    !## put waterlevel down on waterlevel
    RVAL(1,J)=RVAL(2,J-1)
   ENDIF
   IF(J.EQ.NSEG)EXIT
  END DO

  !## interpolate all segments in between and give proper values!
  CALL PCK5RPISG(X,Y,DIST,NXY,NSEG)

  !## determine flow-direction
  LNODAT=.TRUE.
  DO JJ=1,4; IF(RVAL(JJ,1)   .EQ.NODATA)THEN; LNODAT=.FALSE.; EXIT; ENDIF; ENDDO
  DO JJ=1,4; IF(RVAL(JJ,NSEG).EQ.NODATA)THEN; LNODAT=.FALSE.; EXIT; ENDIF; ENDDO
  IF(LNODAT)THEN
   H1=RVAL(1,1)
   H2=RVAL(1,NSEG)
   !## flow direction other-way-around
   J=0
   DO
    J=J+1
    IF(J.GT.NSEG)EXIT
    !## structure found
    IF(IPOS(J).LT.0)THEN
     IF(H2.GT.H1)THEN
      Z          =RVAL(1,J)
      RVAL(1,J)  =RVAL(1,J+1)
      RVAL(1,J+1)=Z
     ELSE
      Z          =RVAL(1,J)
     ENDIF
     J          =J+1
    ENDIF
   ENDDO

   !## interpolate waterlevel,waterbottom,inf.factor,c-value - do not interupt it by structures!
   DO JJ=1,4
    DO J=1,NSEG-1
     IF(IPOS(J).GT.0.AND.RVAL(JJ,J).NE.NODATA)THEN
      DO II=J+1,NSEG
       IF(IPOS(II).GT.0.AND.RVAL(JJ,II).NE.NODATA)EXIT
      ENDDO
      D=MAX(0.0,DIST(II)-DIST(J))
      DO K=J+1,II-1
       FCT=0.0; IF(D.GT.0.0)FCT=(DIST(K)-DIST(J))/D
       RVAL(JJ,K)=RVAL(JJ,J)+((RVAL(JJ,II)-RVAL(JJ,J))*FCT)
      END DO
     ENDIF
    END DO
   ENDDO

   !## interpolate with structures
   IF(NSTW.GT.0)THEN
    DO J=1,NSEG-1
     IF(IPOS(J).NE.0.AND.RVAL(1,J).NE.NODATA)THEN
      D=0.0
      DO II=J+1,NSEG
       D=D+(DIST(II)-DIST(J))
       IF(IPOS(II).NE.0.AND.RVAL(1,II).NE.NODATA)EXIT
      ENDDO
      DO K=J+1,II-1
       FCT=0.0; IF(D.GT.0.0)FCT=(DIST(K)-DIST(J))/D
       RVAL(1,K)=RVAL(1,J)+((RVAL(1,II)-RVAL(1,J))*FCT)
      END DO
     ENDIF
    END DO
   ENDIF

   !## read cross-sectional data
   IREC=ICRS-1
   DO J=1,NCRS
    IF(ISGTYPE.EQ.4)THEN
     READ(ISGIU(5),REC=IREC+J+ICF) DWP(J)%NCROS,IREF,DS
     DWP(J)%DIST=REAL(DS,8)
    ELSE
     READ(ISGIU(5),REC=IREC+J+ICF) DWP(J)%NCROS,IREF,DWP(J)%DIST
    ENDIF
    !## skip two-dimensional cross-sections
    IF(DWP(J)%NCROS.LE.0)THEN; N2DIM=N2DIM+1; CYCLE; ENDIF
    JREC=IREF-1
    DO ICROS=1,DWP(J)%NCROS
     JREC=JREC+1
     IF(ISGTYPE.EQ.4)THEN
      READ(ISGIU(6),REC=JREC+ICF) DS,BS,KS
      CROS(J,ICROS)%DIST=  REAL(DS,8)
      CROS(J,ICROS)%BOTTOM=REAL(BS,8)
      CROS(J,ICROS)%KM=    REAL(KS,8)
     ELSE
      READ(ISGIU(6),REC=JREC+ICF) CROS(J,ICROS)%DIST,CROS(J,ICROS)%BOTTOM,CROS(J,ICROS)%KM
     ENDIF
    END DO
    !## make sure lowest bottom = 0.0
    MIND=MINVAL(CROS(J,1:DWP(J)%NCROS)%BOTTOM)
    DO JJ=1,DWP(J)%NCROS
     CROS(J,JJ)%BOTTOM=CROS(J,JJ)%BOTTOM-MIND
    ENDDO
   END DO

   !## start to intersect all segment/segmentpoints to the model-grid
   ISGLEN=0.0
   DO ISEG=2,NSEG

    X1 =X(ISEG-1); Y1 =Y(ISEG-1); X2 =X(ISEG); Y2 =Y(ISEG)
    !## distance between two points with information
    DXY=((X(ISEG)-X(ISEG-1))**2.0)+((Y(ISEG)-Y(ISEG-1))**2.0)
    IF(DXY.GT.0.0)THEN
     DXY=SQRT(DXY)

     DO JJ=1,4
      RVAL(JJ,0)=(RVAL(JJ,ISEG)-RVAL(JJ,ISEG-1))/DXY
     ENDDO

     N=0
     IF(LQD)THEN
      !## intersect line with rectangular-regular-equidistantial-grid
      CALL IMOD_UTL_INTERSECT_EQUI(DELR(0),DELR(NCOL),DELC(NROW),DELC(0),SIMCSIZE,SIMCSIZE,X1,X2,Y1,Y2,N,.FALSE.)
     ELSE
      CALL IMOD_UTL_INTERSECT_NONEQUI(DELR,DELC,NROW,NCOL,X1,X2,Y1,Y2,N,.FALSE.)
     ENDIF

     !## fill result array
     DXY=0.0
     DO J=1,N

      ISGLEN=ISGLEN+LN(J)
      !## which cross-section is active within current segment
      CALL PCK6RPISG(JCRS,NCRS,ISGLEN)

      IF(LN(J).GT.0.0)THEN
       ICOL=cA(J)
       IROW=rA(J)
       !## within model-domain
       IF(ICOL.GE.1.AND.IROW.GE.1.AND. &
          ICOL.LE.NCOL.AND.IROW.LE.NROW)THEN
        !## interpolate to centre of line ...
        DXY=DXY+(LN(J)/2.0)

        WP=RVAL(1,ISEG-1)+DXY*RVAL(1,0) !waterlevel
        BH=RVAL(2,ISEG-1)+DXY*RVAL(2,0) !waterbottom

        !## correct bottomlevel in case bottomlevel is above waterlevel
        IF(WP-BH.LT.WATERDEPTH)BH=WP-WATERDEPTH

        !## translate from local to global coordinates and get proper wetted perimeter and width of channel!
        CALL PCK4RPISG(JCRS,BH,WP,RWIDTH,WETPER)

        IF(WETPER.GT.0.0)THEN

         NISG = NISG+1

         IF(NISG.GT.SIZE(isgdata,1))THEN
          !## increase memory
          NA=SIZE(ISGDATA,1); MA=SIZE(ISGDATA,2)
          ALLOCATE(ISGDATA_BU(NA+1000,MA)); ISGDATA_BU=0.0
          DO II=1,NA; DO JJ=1,MA; ISGDATA_BU(II,JJ)=ISGDATA(II,JJ); ENDDO; ENDDO
          DEALLOCATE(ISGDATA); ISGDATA=>ISGDATA_BU
         ENDIF
         
          ISGDATA(NISG,1) = ILAY
          ISGDATA(NISG,2) = IROW
          ISGDATA(NISG,3) = ICOL
          ISGDATA(NISG,4) = WP
          C               = RVAL(3,ISEG-1)+DXY*RVAL(3,0)  !C-VALUE
          ISGDATA(NISG,5) = LN(J)*WETPER/C                !CONDUCTANCES
          ISGDATA(NISG,6) = BH
          ISGDATA(NISG,7) = RVAL(4,ISEG-1)+DXY*RVAL(4,0)  !INF.FACTORS

          ! NECESSARY FOR IFVDL=1
          ISGDATA(NISG,8)  = LN(J)
          ISGDATA(NISG,9)  = WETPER
          ISGDATA(NISG,10) = C
          ISGDATA(NISG,11) = RWIDTH

        ENDIF

        DXY=DXY+(LN(J)/2.0)

       ENDIF
      ENDIF
     ENDDO
    ENDIF
   ENDDO
  ENDIF
 ENDIF
ENDDO

 ALLOCATE(IDF(10)); DO I=1,SIZE(IDF); CALL IDFNULLIFY(IDF(I)); IDF(I)%IU=0; ENDDO
 IDF%XMIN=SIMBOX(1); IDF%YMIN=SIMBOX(2); IDF%XMAX=SIMBOX(3); IDF%YMAX=SIMBOX(4); IDF%NODATA=NODATAIDF
 IDF%IXV=0; IDF%ITB=0; IDF%NCOL=NCOL; IDF%NROW=NROW
 IF(LQD)THEN
  IDF%DX=SIMCSIZE; IDF%IEQ=0; IDF%DY=IDF%DX
 ELSE
  IDF%DX=HUGE(1.0); DO I=1,NCOL; IDF%DX=MIN(DELR(I)-DELR(I-1),IDF%DX); ENDDO
  IDF%DY=HUGE(1.0); DO I=1,NROW; IDF%DY=MIN(DELC(I-1)-DELC(I),IDF%DY); ENDDO
  IDF%IEQ=1
  DO I=1,SIZE(IDF)
   ALLOCATE(IDF(I)%SX(0:NCOL),IDF(I)%SY(0:NROW))
   DO J=0,NCOL; IDF(I)%SX(J)=DELR(J); ENDDO
   DO J=0,NROW; IDF(I)%SY(J)=DELC(J); ENDDO
  ENDDO
 ENDIF
 
 DO I=1,SIZE(IDF);
  IF(.NOT.IDFALLOCATEX(IDF(I)))THEN
   CALL IMOD_UTL_PRINTTEXT('Can not allocate memory for idf(i)',1)
  ENDIF
  IDF(I)%X=0.0
 ENDDO

 !## create grid with waterlevels,bottomlevels and resistances
 DO I=1,NISG
  IROW=isgdata(I,2); ICOL=isgdata(I,3)
  !## skip nodata
  IF(isgdata(I,5).EQ.NODATA)CYCLE
  IDF(1) %X(ICOL,IROW)=IDF(1) %X(ICOL,IROW)+ isgdata(I,5)                !## conductance
  IDF(2) %X(ICOL,IROW)=IDF(2) %X(ICOL,IROW)+(isgdata(I,4)*isgdata(I,5))  !## wl
  IDF(3) %X(ICOL,IROW)=IDF(3) %X(ICOL,IROW)+(isgdata(I,6)*isgdata(I,5))  !## bh
  IDF(4) %X(ICOL,IROW)=IDF(4) %X(ICOL,IROW)+(isgdata(I,7)*isgdata(I,5))  !## inf
  IDF(5) %X(ICOL,IROW)=IDF(5) %X(ICOL,IROW)+(isgdata(I,10)*isgdata(I,5)) !## c
  IDF(7) %X(ICOL,IROW)=IDF(7) %X(ICOL,IROW)+(isgdata(I,11)*isgdata(I,5)) !## rwidth
  IDF(9) %X(ICOL,IROW)=IDF(9) %X(ICOL,IROW)+ isgdata(I,8)                !## total length
  IDF(10)%X(ICOL,IROW)=IDF(10)%X(ICOL,IROW)+(isgdata(I,9)*isgdata(I,5))  !## wetper
 ENDDO
 DO IROW=1,NROW; DO ICOL=1,NCOL
  IF(IDF(1)%X(ICOL,IROW).GT.0.0)THEN
   IDF(2) %X(ICOL,IROW)=IDF(2) %X(ICOL,IROW)/IDF(1)%X(ICOL,IROW)
   IDF(3) %X(ICOL,IROW)=IDF(3) %X(ICOL,IROW)/IDF(1)%X(ICOL,IROW)
   IDF(4) %X(ICOL,IROW)=IDF(4) %X(ICOL,IROW)/IDF(1)%X(ICOL,IROW)
   IDF(5) %X(ICOL,IROW)=IDF(5) %X(ICOL,IROW)/IDF(1)%X(ICOL,IROW)
   IDF(7) %X(ICOL,IROW)=IDF(7) %X(ICOL,IROW)/IDF(1)%X(ICOL,IROW)
   IDF(9) %X(ICOL,IROW)=IDF(9) %X(ICOL,IROW) 
   IDF(10)%X(ICOL,IROW)=IDF(10)%X(ICOL,IROW)/IDF(1)%X(ICOL,IROW)
  ENDIF
 ENDDO; ENDDO
  
 !## existence of two-dimensional cross-sections
 IF(N2DIM.GT.0)THEN 

  REWIND(ISGIU(1)); READ(ISGIU(1),*) NISGH
  DO I=1,NISGH
   READ(ISGIU(1),*) CISGH,ISEG,NSEG,ICLC,NCLC,ICRS,NCRS,ISTW,NSTW
   IF(NSEG.LT.2)CYCLE; IF(NCLC.LT.2)CYCLE

   !## read segments
   IREC=ISEG-1
   DO J=1,NSEG
    IF(ISGTYPE.EQ.4)THEN
     READ(ISGIU(2),REC=IREC+J+ICF) XS,YS
     X(J)=XS
     Y(J)=YS
    ELSE
     READ(ISGIU(2),REC=IREC+J+ICF) X(J),Y(J)
    ENDIF
   END DO

   !## evaluate whether line is within model areas, otherwise skip rest of procedure ...
   IF(MAXVAL(X(1:NSEG)).LE.SIMBOX(1).OR.MAXVAL(Y(1:NSEG)).LE.SIMBOX(2).OR. &
      MINVAL(X(1:NSEG)).GE.SIMBOX(3).OR.MINVAL(Y(1:NSEG)).GE.SIMBOX(4))CYCLE

   !## read cross-section information from *.isc
   IF(NCRS.GT.SIZE(DWP))THEN; IF(ALLOCATED(DWP))DEALLOCATE(DWP);ALLOCATE(DWP(NCRS)); ENDIF

   !## read cross-sectional data
   MAXCROS=0; IREC=ICRS-1
   DO J=1,NCRS
    READ(ISGIU(5),REC=IREC+J+ICF) DWP(J)%NCROS
    !## take two-dimensional cross-sections only
    MAXCROS=MAXCROS+MIN(0,DWP(J)%NCROS)
   END DO
   MAXCROS=ABS(MAXCROS)
   !## skip this segment since no two-dimensional cross-sections available
   IF(MAXCROS.EQ.0)CYCLE

   IF(MAXCROS.GT.SIZE(CROS,2).OR.NCRS.GT.SIZE(CROS,1))THEN
    IF(ALLOCATED(CROS))DEALLOCATE(CROS); ALLOCATE(CROS(NCRS,MAXCROS))
   ENDIF

   !## read cross-sectional data
   IREC=ICRS-1
   DO J=1,NCRS
    IF(ISGTYPE.EQ.4)THEN
     READ(ISGIU(5),REC=IREC+J+ICF) DWP(J)%NCROS,IREF,DS,CNAME
     DWP(J)%DIST=REAL(DS,8)
    ELSE
     READ(ISGIU(5),REC=IREC+J+ICF) DWP(J)%NCROS,IREF,DWP(J)%DIST,CNAME
    ENDIF
    IF(DWP(J)%NCROS.GE.0)CYCLE

    IF(DWP(J)%NCROS.GT.0)THEN
     IF(DWP(J)%NCROS.LE.2)THEN
      CALL IMOD_UTL_PRINTTEXT('iMODFLOW can not apply a 1D cross-section with less or equal 2 points',1)
      CALL IMOD_UTL_PRINTTEXT('Cross-section name '//TRIM(CNAME),2)
     ENDIF
    ELSE
     IF(ABS(DWP(J)%NCROS).LE.1)THEN
      CALL IMOD_UTL_PRINTTEXT('iMODFLOW can not apply a 2D cross-section with no locations',1)
      CALL IMOD_UTL_PRINTTEXT('Cross-section name '//TRIM(CNAME),2)
     ENDIF
    ENDIF

    JREC=IREF-1; ICSTYPE=0
    DO ICROS=1,ABS(DWP(J)%NCROS)
     JREC=JREC+1
     IF(ICROS.EQ.1.OR.ICSTYPE.EQ.0)THEN
      IF(ISGTYPE.EQ.4)THEN
       READ(ISGIU(6),REC=JREC+ICF) DS,BS,KS
       CROS(J,ICROS)%DIST=  REAL(DS,8)
       CROS(J,ICROS)%BOTTOM=REAL(BS,8)
       CROS(J,ICROS)%KM=    REAL(KS,8)
      ELSE
       READ(ISGIU(6),REC=JREC+ICF) CROS(J,ICROS)%DIST,CROS(J,ICROS)%BOTTOM,CROS(J,ICROS)%KM
      ENDIF
      !## read bathemetry with thressholds
      IF(CROS(J,ICROS)%DIST.LT.0.AND.CROS(J,ICROS)%BOTTOM.LT.0.0)ICSTYPE=1
      CROS(J,ICROS)%ZP=0.0
     ELSE
      IF(ISGTYPE.EQ.4)THEN
       READ(ISGIU(6),REC=JREC+ICF) DS,BS,ZM,ZC,ZP
       CROS(J,ICROS)%DIST=  REAL(DS,8)
       CROS(J,ICROS)%BOTTOM=REAL(BS,8)
      ELSE
       READ(ISGIU(6),REC=JREC+ICF) CROS(J,ICROS)%DIST,CROS(J,ICROS)%BOTTOM,ZM,ZC,ZP
      ENDIF
      CROS(J,ICROS)%KM=REAL(ZM)+REAL(ZC)/100.0
      CROS(J,ICROS)%ZP=REAL(ZP) 
     ENDIF
    END DO

    CALL ISGADJUSTCOMPUTEXY(J)         !## compute correct x/y coordinate of current cross-section
    CALL UTL_IDFIROWICOL(IDF(1),IROW,ICOL,ISGX,ISGY)               !## get location in raster
    !## skip if outside current model network
    IF(IROW.EQ.0.OR.ICOL.EQ.0)CYCLE
    IF(ISGATTRIBUTES_2DCROSS_READ(J,ICROSS,PCROSS,ZCHK))THEN               !## read bathymetry current cross-section
     WL  =IDF(2)%X(ICOL,IROW)                                  !## waterlevel at cross-section
     INFF=IDF(4)%X(ICOL,IROW)                                  !## infiltration factor at location of cross-section
     C   =IDF(5)%X(ICOL,IROW)                                  !## resistance at location of cross-section
     !## intersection migth miss the cell
     IF(C.LE.0.0)THEN
      !## look around
IRLOOP: DO IR=MAX(1,IROW-1),MIN(NROW,IROW+1)
       DO IC=MAX(1,ICOL-1),MIN(NCOL,ICOL+1)
        !## infiltration factor at location of cross-section
        INFF=IDF(4)%X(IC,IR)                                  
        !## waterlevel at cross-section
        WL=IDF(2)%X(IC,IR)   
        !## resistance
        C=IDF(5)%X(IC,IR)
        IF(C.NE.0.0)EXIT IRLOOP
       ENDDO
      ENDDO IRLOOP
     ENDIF
     CALL ISG2GRID_BATHEMETRY(IDF,SIZE(IDF),ICROSS,PCROSS,ZCHK,WL,C,INFF)  !## adjust stage grid for bathemetrie
     CALL IDFDEALLOCATEX(ICROSS); CALL IDFDEALLOCATEX(PCROSS)
    ENDIF

   END DO
  ENDDO
 
 ENDIF
 
 !## extent grids based upon their width
 CALL ISG2GRID_EXTENT_WITH_WIDTH(SIZE(IDF),IDF,MAXWIDTH)

 !## clean list
 DO J=1,SIZE(isgdata,1); DO I=1,SIZE(isgdata,2); isgdata(J,I)=0.0; ENDDO; ENDDO

 !## reuse isg from the grid, that is the consequence of using 2d cross-sections
 NISG=0
 DO IROW=1,NROW; DO ICOL=1,NCOL
  IF(IDF(1)%X(ICOL,IROW).GT.0.0)THEN
   NISG=NISG+1
   isgdata(NISG,1) = ILAY
   isgdata(NISG,2) = IROW
   isgdata(NISG,3) = ICOL
   isgdata(NISG,5)= IDF(1)%X(ICOL,IROW)  ! cond
   isgdata(NISG,4)= IDF(2)%X(ICOL,IROW)  ! stage
   isgdata(NISG,6)= IDF(3)%X(ICOL,IROW)  ! rbot
   isgdata(NISG,7)= IDF(4)%X(ICOL,IROW)  ! inf.
   isgdata(NISG,8)= IDF(9)%X(ICOL,IROW)  ! length
   isgdata(NISG,9)= IDF(10)%X(ICOL,IROW) ! wettedperimeter
   isgdata(NISG,10)=IDF(5)%X(ICOL,IROW)  ! resis

  ENDIF
 ENDDO; ENDDO

 CALL IDFDEALLOCATE(IDF,SIZE(IDF))
 CALL IMOD_UTL_PRINTTEXT('Finished gridding 2d cross-sections, no. of isg elements '//TRIM(IMOD_UTL_ITOS(NISG)),3)

DO I=1,8; CLOSE(ISGIU(I)); END DO

IF(ALLOCATED(QSORT))DEALLOCATE(QSORT); IF(ALLOCATED(DIST))DEALLOCATE(DIST)
IF(ALLOCATED(CROS))DEALLOCATE(CROS);   IF(ALLOCATED(RVAL))DEALLOCATE(RVAL)
IF(ALLOCATED(IPOS))DEALLOCATE(IPOS);   IF(ALLOCATED(DWP))DEALLOCATE(DWP)
call INTERSECT_DEALLOCATE()

RETURN
END SUBROUTINE

!###====================================================================
SUBROUTINE PCK2RPISG(SDATE,EDATE,TTIME,QSORT,XNR,LU,MAXITEMS,RVAL, &
                     NR,IREF,IISGH,ITYPE,NODATA,ISGTYPE)
!###====================================================================
USE MOD_GLBVAR, ONLY : LINE,ISS
USE IMOD_UTL, ONLY: IMOD_UTL_IDATETOJDATE,IMOD_UTL_GETMED2,IMOD_UTL_ITOS,IMOD_UTL_PRINTTEXT,ICF
IMPLICIT NONE
REAL(KIND=8),INTENT(IN) :: NODATA
INTEGER,INTENT(IN) :: SDATE,EDATE,TTIME,MAXITEMS,LU,IISGH,NR,IREF,ITYPE,ISGTYPE
REAL(KIND=8),INTENT(INOUT),DIMENSION(TTIME,MAXITEMS) :: QSORT
REAL(KIND=8),DIMENSION(MAXITEMS),INTENT(INOUT) :: RVAL,XNR
REAL(KIND=4),DIMENSION(MAXITEMS) :: SVAL
INTEGER :: IR,I,J,I1,I2,IDATE,NDATE,NAJ,IREC,N
INTEGER :: IAVERAGE=2 !## median values
REAL(KIND=8) ::  X
CHARACTER(LEN=8) :: CTIME !## time in hhmmss

IREC=IREF-1

IF(NR.GT.0)THEN

 IF(ISS.EQ.1)QSORT= 0.0
 IF(ISS.EQ.2)QSORT=NODATA 

 I1=1; IR=0; XNR=0.0
 
 DO 
  
  IR=IR+1; IF(IR.GT.NR)EXIT
    
  IREC=IREC+1

  RVAL=0.0
  IF(ITYPE.EQ. 1)THEN
   IF(ISGTYPE.EQ.4)THEN
    READ(LU,REC=IREC+ICF) IDATE,(SVAL(I),I=1,MAXITEMS)
    DO I=1,MAXITEMS; RVAL(I)=REAL(SVAL(I),8); ENDDO
   ELSE
    READ(LU,REC=IREC+ICF) IDATE,CTIME,(RVAL(I),I=1,MAXITEMS)
   ENDIF
  ENDIF
  IF(ITYPE.EQ.-1)THEN
   IF(ISGTYPE.EQ.4)THEN
    READ(LU,REC=IREC+ICF) IDATE,(SVAL(I),I=1,2)
    DO I=1,2; RVAL(I)=REAL(SVAL(I),8); ENDDO
   ELSE
    READ(LU,REC=IREC+ICF) IDATE,CTIME,(RVAL(I),I=1,2)
   ENDIF
  ENDIF
!  IF(ITYPE.EQ. 1)READ(LU,REC=IREC+ICF) IDATE,(RVAL(I),I=1,MAXITEMS)
!  IF(ITYPE.EQ.-1)READ(LU,REC=IREC+ICF) IDATE,(RVAL(I),I=1,2)

  !## don't bother for steady-state, take the mean!
  IF(ISS.EQ.1)THEN
   DO I=1,MAXITEMS; IF(RVAL(I).NE.NODATA)THEN; QSORT(1,I)=QSORT(1,I)+RVAL(I); XNR(I)=XNR(I)+1.0; ENDIF; ENDDO
  !## transient simulation
  ELSEIF(ISS.EQ.2)THEN
   
   IDATE=IMOD_UTL_IDATETOJDATE(IDATE)
   
   !## try to speed up the search for transient searches, only for more than 10 entries
   IF(IR.EQ.1.AND.NR.GT.10)THEN
    READ(LU,REC=IREC+1+ICF) NDATE
    !## number of days in between
    N=IMOD_UTL_IDATETOJDATE(NDATE)-IDATE
    N=(SDATE-IDATE)/N; N=MAX(N-1,0)
    IF(IR+N.LT.NR)THEN
     READ(LU,REC=IREC+1+N+ICF) NDATE
     !## okay, take this jump in time
     IF(IMOD_UTL_IDATETOJDATE(NDATE).LE.SDATE)THEN
      IR  =IR  +N
      IREC=IREC+N
     ENDIF
    ENDIF
   ENDIF

   NDATE=EDATE
   IF(IR.LT.NR)THEN
    READ(LU,REC=IREC+1+ICF) NDATE
    NDATE=IMOD_UTL_IDATETOJDATE(NDATE)
   ENDIF
   NDATE=MIN(NDATE,EDATE)

   !## stop searchin for data, outside modeling window!
   IF(IDATE.GT.EDATE)EXIT

   !## within modeling window
   IF(NDATE.GE.SDATE)THEN

    N=NDATE-SDATE
    IF(IDATE.GT.SDATE)N=N-(IDATE-SDATE)
    I2=I1+N-1

    DO I=1,MAXITEMS
     QSORT(I1:I2,I)=RVAL(I)
    END DO

    I1=I2+1

   ENDIF
  ENDIF

 END DO

 IF(ISS.EQ.1)THEN
  !## determine for each period appropriate attribute term
  IF(NR.EQ.0)CALL IMOD_UTL_PRINTTEXT('No data found for steady-state in ISG file!',2)
  DO I=1,MAXITEMS
   IF(XNR(I).GT.0.0)THEN
    RVAL(I)=QSORT(1,I)/XNR(I)
   ELSE
    RVAL(I)=NODATA
   ENDIF
  ENDDO
 ELSEIF(ISS.EQ.2)THEN
  IF(IAVERAGE.EQ.1)THEN
   !## arithmetic mean
   XNR=0.0; RVAL=0.0
   DO I=1,MAXITEMS
    DO J=1,TTIME
     IF(QSORT(J,I).NE.NODATA)THEN
      RVAL(I)=RVAL(I)+QSORT(J,I)
      XNR(I) =XNR(I)+1.0
     ENDIF
    ENDDO
    IF(XNR(I).GT.0.0)THEN
     RVAL(I)=RVAL(I)/XNR(I)
    ELSE
     RVAL(I)=NODATA
    ENDIF
   ENDDO 
  ELSEIF(IAVERAGE.EQ.2)THEN
   !## median
   DO I=1,MAXITEMS
    RVAL(I)=IMOD_UTL_GETMED2(QSORT(:,I),TTIME,NODATA,NAJ) 
    IF(NAJ.EQ.0)CALL IMOD_UTL_PRINTTEXT('No data found for current stress period in ISG file!',2)
   ENDDO
  ENDIF
 ENDIF

ELSE
 IF(ITYPE.EQ.1)CALL IMOD_UTL_PRINTTEXT('No waterlevels found for SEGMENT '//TRIM(IMOD_UTL_ITOS(IISGH)),2)
ENDIF

RETURN
END SUBROUTINE

!###====================================================================
SUBROUTINE PCK4RPISG(JCRS,BH,WP,RWIDTH,WETPER)
!###====================================================================
!RE.: In case waterlevel exceeds cross-section, rwidth and wetper will
!     yield max. values based upon cross-section information only!
!###====================================================================
USE MOD_ISGVAR, ONLY : DWP,CROS
IMPLICIT NONE
INTEGER,INTENT(IN) :: JCRS
REAL(KIND=8),INTENT(IN) :: BH,WP
REAL(KIND=8),INTENT(OUT) :: RWIDTH,WETPER
INTEGER :: ICROS
REAL(KIND=8) ::  DISTB,DISTW,FACTOR,X1,X2,WPCOR

WPCOR =WP-BH
WETPER=0.0

!## use default cross-section whenever no cross-section has been fouond
IF(JCRS.EQ.0)THEN
 RWIDTH=5.0
 WETPER=WPCOR+RWIDTH
 RETURN
ENDIF

!## find left-x-coordinate
X1=CROS(JCRS,1)%DIST
DO ICROS=1,DWP(JCRS)%NCROS-1
 IF(CROS(JCRS,ICROS)%BOTTOM  .GE.WPCOR.AND. &
    CROS(JCRS,ICROS+1)%BOTTOM.LE.WPCOR)THEN
  DISTB  = CROS(JCRS,ICROS)%BOTTOM-CROS(JCRS,ICROS+1)%BOTTOM
  DISTW  = CROS(JCRS,ICROS)%BOTTOM-WPCOR
  
  IF(DISTB.NE.0.0)THEN
  FACTOR = DISTW/DISTB
  X1     = FACTOR*(CROS(JCRS,ICROS+1)%DIST-CROS(JCRS,ICROS)%DIST)
  X1     = CROS(JCRS,ICROS)%DIST+X1

  !## wetted perimeter
  FACTOR  =(X1-CROS(JCRS,ICROS+1)%DIST)**2.0+(WPCOR-CROS(JCRS,ICROS+1)%BOTTOM)**2.0
  IF(FACTOR.NE.0.0)FACTOR=SQRT(FACTOR)
  WETPER=WETPER+FACTOR
  ENDIF
  
 ENDIF
END DO

!## find right-x-coordinate
X2=CROS(JCRS,DWP(JCRS)%NCROS)%DIST
DO ICROS=DWP(JCRS)%NCROS-1,1,-1
 IF(CROS(JCRS,ICROS+1)%BOTTOM.GE.WPCOR.AND. &
    CROS(JCRS,ICROS)%BOTTOM  .LE.WPCOR)THEN
  DISTB  = CROS(JCRS,ICROS+1)%BOTTOM-CROS(JCRS,ICROS)%BOTTOM
  DISTW  = CROS(JCRS,ICROS+1)%BOTTOM-WPCOR

  IF(DISTB.NE.0.0)THEN
  FACTOR = DISTW/DISTB
  X2     = FACTOR*(CROS(JCRS,ICROS+1)%DIST-CROS(JCRS,ICROS)%DIST)
  X2     = CROS(JCRS,ICROS+1)%DIST-X2

  !## wetted perimeter
  FACTOR =(X2-CROS(JCRS,ICROS)%DIST)**2.0+(WPCOR-CROS(JCRS,ICROS)%BOTTOM)**2.0
  IF(FACTOR.NE.0.0)FACTOR=SQRT(FACTOR)
  WETPER=WETPER+FACTOR
  ENDIF
  
 ENDIF
END DO

!## proces wetted perimeter for 'in-between' sections
DO ICROS=1,DWP(JCRS)%NCROS-1
 IF(CROS(JCRS,ICROS)%BOTTOM  .LT.WPCOR.AND. &
    CROS(JCRS,ICROS+1)%BOTTOM.LT.WPCOR)THEN
  !## wetted perimeter
  FACTOR =(CROS(JCRS,ICROS)%DIST-CROS(JCRS,ICROS+1)%DIST)**2.0+ &
          (CROS(JCRS,ICROS)%BOTTOM-CROS(JCRS,ICROS+1)%BOTTOM)**2.0
  IF(FACTOR.NE.0.0)FACTOR=SQRT(FACTOR)
  WETPER=WETPER+FACTOR
 ENDIF
END DO

RWIDTH=X2-X1

RETURN
END SUBROUTINE

!###====================================================================
SUBROUTINE PCK5RPISG(X,Y,DIST,NXY,NSEG)
!###====================================================================
IMPLICIT NONE
INTEGER,INTENT(IN) :: NXY,NSEG
REAL(KIND=8),INTENT(IN),DIMENSION(NXY) :: X,Y
REAL(KIND=8),INTENT(OUT),DIMENSION(NXY) :: DIST
INTEGER :: ISEG
REAL(KIND=8) ::  DXY

DIST=0.0
DO ISEG=2,NSEG
 DXY=((X(ISEG)-X(ISEG-1))**2.0)+((Y(ISEG)-Y(ISEG-1))**2.0)
 IF(DXY.GT.0.0)DXY=SQRT(DXY)
 DIST(ISEG)=DIST(ISEG-1)+DXY
END DO

RETURN
END SUBROUTINE

!###====================================================================
SUBROUTINE PCK6RPISG(JCRS,NCRS,ISGLEN)
!###====================================================================
USE MOD_ISGVAR, ONLY : DWP
IMPLICIT NONE
REAL(KIND=8),INTENT(IN) :: ISGLEN
INTEGER,INTENT(IN) :: NCRS
INTEGER,INTENT(OUT) :: JCRS
INTEGER :: IC1,IC2,ICRS
REAL(KIND=8) ::  DC,D,DCL,DCR

JCRS=0

IF(NCRS.EQ.1.AND.DWP(1)%NCROS.GT.0)THEN; JCRS=1; RETURN; ENDIF

IC1 =0  !## cross-section left
IC2 =0  !## cross-section right
DCL =10.0E10  !left
DCR =10.0E10  !right
DC  =10.0E10
DO ICRS=1,NCRS
 !## take only one-dimensional cross-sections
 IF(DWP(ICRS)%NCROS.LE.0)CYCLE
 D=DWP(ICRS)%DIST-ISGLEN
 IF(D.GT.0)THEN !## lies right
  !## take whenever closer than next cross-section in the front
  IF(D.LE.DCL)THEN
   DCL=D
   IC2=ICRS
  ENDIF
 ELSE           !## lies behind
  !## take whenever closer than previous cross-section in the back
  IF(ABS(D).LE.DCR)THEN
   DCR=ABS(D)
   IC1=ICRS
  ENDIF
 ENDIF
ENDDO
IF(IC1.EQ.0)IC1=IC2
JCRS=IC1

RETURN
END SUBROUTINE

!###====================================================================
SUBROUTINE PCK7RPISG(ITYP,LU,NTYP,NXY,NSEG,X,Y,DIST,IPOS,ITYPE,ISGTYPE)
!###====================================================================
USE IMOD_UTL, ONLY: ICF
IMPLICIT NONE
INTEGER,INTENT(IN) :: NXY,ITYP,LU,NTYP,ITYPE,ISGTYPE
INTEGER,INTENT(INOUT) :: NSEG
REAL(KIND=8),INTENT(INOUT),DIMENSION(NXY) :: X,Y
REAL(KIND=8),INTENT(INOUT),DIMENSION(NXY) :: DIST
INTEGER,INTENT(OUT),DIMENSION(NXY) :: IPOS
INTEGER :: IREC,ISEG,IREF,I,J,II
REAL(KIND=8) ::  DXY,D1,D2,F,XC,YC,SDIST
REAL(KIND=4) :: DS

!## include calculation nodes as segments!
IREC=ITYP-1
!## determine which nodes possess heads/etc.
DO I=1,NTYP
 !## get tot distances
 CALL PCK5RPISG(X,Y,DIST,NXY,NSEG)
 IREC=IREC+1
 IF(ISGTYPE.EQ.4)THEN
  READ(LU,REC=IREC+ICF) J,IREF,DS
  SDIST=REAL(DS,8)
 ELSE
  READ(LU,REC=IREC+ICF) J,IREF,SDIST
 ENDIF
 DXY=SDIST  !k=number,dxy=distance from beginning

 !## find position in between segments
 DO ISEG=2,NSEG
  IF(DXY.GE.DIST(ISEG-1).AND.DXY.LE.DIST(ISEG))EXIT
 END DO

 !## caused by inaccuracy of comparison of dxy and dist()
 ISEG=MIN(ISEG,NSEG)
 !## distance current segment
 D1= DIST(ISEG)-DIST(ISEG-1)
 D2= DXY-DIST(ISEG-1)
 !## segment itself is zero
 F =0.0
 IF(D1.NE.0.0)F=D2/D1

 !## put in extra coordinate
 IF(F.LE.0.0.AND.ITYPE.EQ.1)THEN
  IPOS(ISEG-1)=IREC*ITYPE  !##put data to current node
 ELSEIF(F.GE.1.0.AND.ITYPE.EQ.1)THEN
  IPOS(ISEG)=IREC*ITYPE    !##put data to current node
 ELSE
  XC=X(ISEG-1)+((X(ISEG)-X(ISEG-1))*F)
  YC=Y(ISEG-1)+((Y(ISEG)-Y(ISEG-1))*F)
  !## position coordinates in between
  DO II=NSEG+1,ISEG+1,-1
   X(II)   =X(II-1)
   Y(II)   =Y(II-1)
   IPOS(II)=IPOS(II-1)
  ENDDO
  X(ISEG)            =XC
  Y(ISEG)            =YC
  !## increase number of segments
  NSEG               =NSEG+1     
  !##put data to current node
  IPOS(ISEG)         =IREC*ITYPE 
 ENDIF

!## duplicate point in case of structure
 IF(ITYPE.EQ.-1)THEN
  !## position coordinates in between
  DO II=NSEG+1,ISEG+1,-1
   X(II)   =X(II-1)
   Y(II)   =Y(II-1)
   IPOS(II)=IPOS(II-1)
  ENDDO 
  X(ISEG)            =XC
  Y(ISEG)            =YC
  !## increase number of segments
  NSEG               =NSEG+1     
  !## put data to current node
  IPOS(ISEG)         =I*ITYPE   
 ENDIF

END DO

!## put last ipos() for secerity reasons ... could happen in case of accuray of coordinates
IF(IPOS(NSEG).EQ.0)THEN
 DO I=NSEG,1,-1
  IF(IPOS(I).GT.0)THEN
   IPOS(NSEG)=IPOS(I)
   EXIT
  ENDIF
 END DO
ENDIF

RETURN
END SUBROUTINE

!###====================================================================
SUBROUTINE PCK3RPISG(ISGLIST,NISG)
!###====================================================================
USE IMOD_UTL, ONLY : IMOD_UTL_PRINTTEXT,IMOD_UTL_RTOS
USE MOD_BASVAR, ONLY : DELR,DELC
USE GLOBAL, ONLY: CC, CV, BOTM, LBOTM, KDSV
USE GWFRIVMODULE, ONLY: IFVDL,ISFT,SFT
IMPLICIT NONE
INTEGER, INTENT(IN) :: NISG
REAL(KIND=4),DIMENSION(NISG,11),INTENT(INOUT) :: ISGLIST
REAL(KIND=4) ::  LEKCONDUCTANCE_Y
INTEGER :: I, ILAY, IROW, ICOL
REAL(KIND=4) ::  DXY,LI,BIN,C0,C1,KH,KV,COND_IN,FCT,D
REAL(KIND=4),PARAMETER :: TINY=1.0D-20
LOGICAL :: LUSED ! PKS

IF (IFVDL.EQ.0) RETURN

DO I = 1, NISG
 ILAY = ISGLIST(I,1)
 IROW = ISGLIST(I,2)
 ICOL = ISGLIST(I,3)
 DXY=(DELR(ICOL)-DELR(ICOL-1))*(DELC(IROW-1)-DELC(IROW))
 CALL PKS7MPITRN(ICOL,IROW,ILAY,LUSED) ! PKS
 IF (LUSED) THEN ! PKS
 FCT=ISGLIST(I,7)   !## infiltration factor
 LI =ISGLIST(I,8)   !## length of river segment
 BIN=ISGLIST(I,9)   !## wetted perimeter
 C0 = ISGLIST(I,10) !## resistance
 C1=CV(ICOL,IROW,ILAY)
 IF(C1.GT.0.0)THEN
  C1 =DXY/C1                                    !## resistance of aquitard (1/day -> day)
  IF(ISFT.EQ.1)THEN
   D  =SFT(ICOL,IROW,1)                         !## thickness of aquifer
   KH =SFT(ICOL,IROW,2)                         !## permeability
  ELSE
   D = BOTM(ICOL,IROW,LBOTM(ILAY)-1)-BOTM(ICOL,IROW,LBOTM(ILAY)) !## thickness of aquifer
   KH =KDSV(ICOL,IROW,ILAY)/(D+TINY)              !## permeability
  ENDIF
  KV =KH/10.0                                   !## vertical permeability

  IF(KH.LE.0.0D0)CALL IMOD_UTL_PRINTTEXT('Error KH='//TRIM(IMOD_UTL_RTOS(KH,'F',7)),2)
  IF(D .LE.0.0D0)CALL IMOD_UTL_PRINTTEXT('Error D ='//TRIM(IMOD_UTL_RTOS(D ,'F',7)),2)

  !## process formulae van de Lange - infiltratie weerstand
  COND_IN=0.0; IF(FCT.GT.0.0.AND.FCT.LE.1.0)COND_IN=LEKCONDUCTANCE_Y(DXY,D,KV,KH,C1,LI,BIN,C0/FCT)
  !## process formulae van de Lange - drainage weerstand
  ISGLIST(I,5)=LEKCONDUCTANCE_Y(DXY,D,KV,KH,C1,LI,BIN,C0)
  !## infiltration factor
  ISGLIST(I,7)=FCT; IF(ISGLIST(I,5).GT.0.0)ISGLIST(I,7)=COND_IN/ISGLIST(I,5)

  IF(ISGLIST(I,7).GT.1.0.OR.ISGLIST(I,7).LT.0.0)THEN
   CALL IMOD_UTL_PRINTTEXT('Error C0/FCT    ='//TRIM(IMOD_UTL_RTOS(C0/FCT,'F',7)),0)
   CALL IMOD_UTL_PRINTTEXT('Error COND_INF  ='//TRIM(IMOD_UTL_RTOS(COND_IN,'F',7)),0)
   CALL IMOD_UTL_PRINTTEXT('Error C0        ='//TRIM(IMOD_UTL_RTOS(C0,'F',7)),0)
   CALL IMOD_UTL_PRINTTEXT('Error COND_DRN  ='//TRIM(IMOD_UTL_RTOS(ISGLIST(NISG,5),'F',7)),0)
   CALL IMOD_UTL_PRINTTEXT('Error Inf.Factor='//TRIM(IMOD_UTL_RTOS(ISGLIST(NISG,7),'F',7)),0)
  ENDIF
 ELSE
  ISGLIST(I,5)=0.0
 ENDIF
 ENDIF !PKS
END DO

RETURN
END SUBROUTINE

!###===============================================================================
 SUBROUTINE ISGADJUSTCOMPUTEXY(J)
 !###===============================================================================
 USE MOD_ISGVAR
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: J
 INTEGER :: I,K
 REAL(KIND=8) ::  DXY,F,TD,DIST

 DIST=DWP(J)%DIST

 TD=0.0
 DO I=2,NSEG
  DXY=((X(I)-X(I-1))**2.0)+((Y(I)-Y(I-1))**2.0)
  IF(DXY.GT.0.0)DXY=SQRT(DXY)
  TD=TD+DXY; IF(TD.GE.DIST)EXIT
 END DO

 !## distance current segment
 F   =(DIST-(TD-DXY))/DXY
 ISGX= X(I-1)+(X(I)-X(I-1))*F
 ISGY= Y(I-1)+(Y(I)-Y(I-1))*F

 END SUBROUTINE ISGADJUSTCOMPUTEXY

 !###======================================================================
 LOGICAL FUNCTION ISGATTRIBUTES_2DCROSS_READ(J,IDF,PIDF,ZCHK)
 !###======================================================================
 USE IMOD_UTL, ONLY : IMOD_UTL_PRINTTEXT
 USE IMOD_IDF_PAR
 USE IMOD_IDF
 USE MOD_ISGVAR
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: J
 REAL(KIND=8),INTENT(OUT) :: ZCHK
 INTEGER :: I,IROW,ICOL,ICHK
 TYPE(IDFOBJ),INTENT(OUT) :: IDF,PIDF
 LOGICAL :: ISGATTRIBUTES_2DCROSS_IDFDIM

 ISGATTRIBUTES_2DCROSS_READ=.FALSE.

 CALL IDFNULLIFY(IDF); CALL IDFNULLIFY(PIDF)
 IF(.NOT.ISGATTRIBUTES_2DCROSS_IDFDIM(IDF,J))THEN
  CALL IMOD_UTL_PRINTTEXT('iMOD can not present results in table. No cellsizes can be determined',1)
 ENDIF

 CALL IDFCOPY(IDF,PIDF)
 IF(.NOT.IDFALLOCATEX(IDF))RETURN
 IF(.NOT.IDFALLOCATEX(PIDF))RETURN
 IDF%NODATA=-9999; IDF%X=IDF%NODATA
 PIDF%NODATA=-9999; PIDF%X=0.0  !## on default, inundated whenever stage higher than riverbed

 ICHK=0; IF(CROS(J,1)%DIST.LT.0.0.AND.CROS(J,1)%BOTTOM.LT.0.0)THEN
  ICHK=1; ZCHK=CROS(J,1)%KM
 ENDIF

 !## skip first (dx,dy)
 DO I=2,ABS(DWP(J)%NCROS)
  CALL UTL_IDFIROWICOL(IDF,IROW,ICOL,CROS(J,I)%DIST,CROS(J,I)%BOTTOM) 
  IDF%X(ICOL,IROW)=CROS(J,I)%KM 
  !## store inundation pointer
  PIDF%X(ICOL,IROW)=1.0; IF(ICHK.EQ.1)PIDF%X(ICOL,IROW)=CROS(J,I)%ZP
 ENDDO 

 ISGATTRIBUTES_2DCROSS_READ=.TRUE.

 END FUNCTION ISGATTRIBUTES_2DCROSS_READ

 !###======================================================================
 LOGICAL FUNCTION ISGATTRIBUTES_2DCROSS_IDFDIM(IDF,J)
 !###======================================================================
 USE IMOD_IDF_PAR
 USE MOD_ISGVAR
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: J
 TYPE(IDFOBJ),INTENT(OUT) :: IDF
 INTEGER :: I,II
 REAL(KIND=8) ::  D

 ISGATTRIBUTES_2DCROSS_IDFDIM=.FALSE.

 IDF%DX=ABS(CROS(J,1)%DIST)
 IDF%DY=ABS(CROS(J,1)%BOTTOM)

 IF(IDF%DX.LE.0.0.OR.IDF%DY.LE.0.0)RETURN

 IDF%XMIN=MINVAL(CROS(J,2:ABS(DWP(J)%NCROS))%DIST)  -0.5*IDF%DX
 IDF%XMAX=MAXVAL(CROS(J,2:ABS(DWP(J)%NCROS))%DIST)  +0.5*IDF%DX
 IDF%YMIN=MINVAL(CROS(J,2:ABS(DWP(J)%NCROS))%BOTTOM)-0.5*IDF%DY
 IDF%YMAX=MAXVAL(CROS(J,2:ABS(DWP(J)%NCROS))%BOTTOM)+0.5*IDF%DY

 IDF%IEQ=0; IDF%ITB=0; IDF%IXV=0
 IDF%NCOL=(IDF%XMAX-IDF%XMIN)/IDF%DX
 IDF%NROW=(IDF%YMAX-IDF%YMIN)/IDF%DY

 ISGATTRIBUTES_2DCROSS_IDFDIM=.TRUE.

 END FUNCTION ISGATTRIBUTES_2DCROSS_IDFDIM

 !###====================================================================
 SUBROUTINE ISG2GRID_BATHEMETRY(IDF,NIDF,ICROSS,PCROSS,ZCHK,WL,C,INFF)
 !###====================================================================
 USE IMOD_IDF_PAR
 USE IMOD_IDF, ONLY : IDFGETAREA
 USE IMOD_UTL, ONLY : UTL_IDFIROWICOL,UTL_IDFGETLOC
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NIDF
 TYPE(IDFOBJ),DIMENSION(NIDF),INTENT(INOUT) :: IDF
 TYPE(IDFOBJ),INTENT(IN) :: ICROSS,PCROSS
 REAL(KIND=8),INTENT(IN) :: WL,C,INFF,ZCHK
 INTEGER :: IR1,IR2,IC1,IC2,IROW,ICOL,JROW,JCOL
 REAL(KIND=8) ::  XC,YC,CR

 !## fill in 3d-cross-section into model-network
 DO IROW=1,ICROSS%NROW; DO ICOL=1,ICROSS%NCOL

  !## nodata found for waterlevel
  IF(ICROSS%X(ICOL,IROW).EQ.ICROSS%NODATA)CYCLE
  !## waterlevel equal or less than riverbot
  IF(WL.LE.ICROSS%X(ICOL,IROW))CYCLE
  !## if inundation-criterion applied, only inundate if zchk criterion is met
  IF(PCROSS%X(ICOL,IROW).LT.0.0.AND.WL.LE.ZCHK)CYCLE

  !## all okay - continue
  CALL UTL_IDFGETLOC(ICROSS,IROW,ICOL,XC,YC)
  CALL UTL_IDFIROWICOL(IDF(1),JROW,JCOL,XC,YC)
  !## outside current model
  IF(JROW.EQ.0.OR.JCOL.EQ.0)CYCLE

  !## manipulate resistance for ...
  CR=C; IF(PCROSS%X(ICOL,IROW).NE.0.0)CR=CR*ABS(PCROSS%X(ICOL,IROW))
  IF(IDF(9)%X(JCOL,JROW).EQ.0.0)THEN
   IDF(1)%X(JCOL,JROW)=IDFGETAREA(ICROSS,ICOL,IROW)/CR
  ELSE
   IDF(1)%X(JCOL,JROW)=IDF(1)%X(JCOL,JROW)+IDFGETAREA(ICROSS,ICOL,IROW)/CR
  ENDIF
  IDF(9)%X(JCOL,JROW)=IDF(9)%X(JCOL,JROW)+1.0 !## counter how many times it passes
  IDF(2)%X(JCOL,JROW)=WL                      !## waterlevel
  IDF(3)%X(JCOL,JROW)=ICROSS%X(ICOL,IROW)     !## bottomlevel
  IDF(4)%X(JCOL,JROW)=INFF                    !## infiltration factor
 ENDDO; ENDDO

 END SUBROUTINE ISG2GRID_BATHEMETRY
 
 !###====================================================================
 SUBROUTINE ISG2GRID_EXTENT_WITH_WIDTH(NIDF,IDF,MAXWIDTH)
 !###====================================================================
 USE IMOD_IDF_PAR
 USE IMOD_IDF, ONLY : IDFGETAREA
 IMPLICIT NONE
 REAL(KIND=8),PARAMETER :: FDX=1.5D0 !## apply extent whenever width is 1.5*dx cell
 INTEGER,INTENT(IN) :: NIDF
 REAL(KIND=8),INTENT(IN) :: MAXWIDTH
 TYPE(IDFOBJ),INTENT(INOUT),DIMENSION(NIDF) :: IDF
 INTEGER :: IROW,ICOL,N,NN,IR,IC,IRR,ICC,I,IRAT,IRAT1,NRE,NCE
 REAL(KIND=8) ::  W,L,X,Y,X1,Y1,F,A
 logical :: pks7mpimasterwrite ! PKS
 
 DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL

  !## already visited by bathymetry routine - so skip it
  IF(IDF(6)%X(ICOL,IROW).NE.0.0)THEN
   IDF(6)%X(ICOL,IROW)=-1.0*IDF(6)%X(ICOL,IROW); CYCLE
  ENDIF

  !## get aggregation (erosion-matrix) pointer to one for all initially
  IF(IDF(1)%X(ICOL,IROW).NE.IDF(1)%NODATA)IDF(6)%X(ICOL,IROW)=1.0
 
  !## conductance assume to be cell-filled; f=dx/length-segment*dx/width
  A=IDFGETAREA(IDF(1),ICOL,IROW)  
  W=IDF(7)%X(ICOL,IROW)
  L=IDF(9)%X(ICOL,IROW)
       
  !## take into account the fact that the line could be diagonal
  IF(W.GT.FDX*SQRT(A))THEN
   !## conductance as 1/d
   IDF(1)%X(ICOL,IROW)=IDF(1)%X(ICOL,IROW)/(L*W)
   !## fill this one
   IDF(1)%X(ICOL,IROW)=A*IDF(1)%X(ICOL,IROW)
  ENDIF
 ENDDO; ENDDO
 
 if (pks7mpimasterwrite()) WRITE(*,'(A)') 'Start computing erosion matrix' 
 DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL

  !## already processed by bathymetry
  IF(IDF(6)%X(ICOL,IROW).LT.0.0)CYCLE
  !## river available
  IF(IDF(1)%X(ICOL,IROW).NE.IDF(1)%NODATA)THEN

   !## get mean width up to maximum width
   W=MIN(IDF(7)%X(ICOL,IROW),MAXWIDTH)
   !## current cellsize is less than current width of river
   A=IDFGETAREA(IDF(1),ICOL,IROW)  

   IF(W.GT.FDX*SQRT(A))THEN 

    !## create antialiased erosion/fattening matrix
    CALL ISG2GRID_EROSION_MATRIX(A,W,NIDF,IDF,ICOL,IROW) 

   ENDIF
  ENDIF
 ENDDO; ENDDO
 if (pks7mpimasterwrite()) WRITE(*,'(A)') 'Finished computing erosion matrix' 
 
 !## clean for conductances le zero
 DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL
  IF(IDF(1)%X(ICOL,IROW).LE.0.0)THEN
   DO I=1,4; IDF(I)%X(ICOL,IROW)=IDF(I)%NODATA; ENDDO
  ENDIF
 ENDDO; ENDDO

 END SUBROUTINE ISG2GRID_EXTENT_WITH_WIDTH

 !###====================================================================
 SUBROUTINE ISG2GRID_EROSION_MATRIX(A,W,NIDF,IDF,ICOL,IROW)
 !###====================================================================
 USE IMOD_IDF_PAR
 USE IMOD_IDF, ONLY : IDFGETAREA,IDFGETEDGE
 USE IMOD_UTL, ONLY : UTL_IDFIROWICOL,UTL_IDFGETLOC
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: ICOL,IROW,NIDF
 TYPE(IDFOBJ),INTENT(INOUT),DIMENSION(NIDF) :: IDF
 REAL(KIND=8),INTENT(IN) :: W,A
 INTEGER :: I,J,K,IC,IR,IC1,IC2,IR1,IR2
 REAL(KIND=8) ::  X,Y,X1,X2,Y1,Y2,TRAD,C,F,DX,DY
 REAL(KIND=8),DIMENSION(9) :: XP,YP
  
 !## half of the river width
 TRAD=W/2.0

 !## get current location
 CALL UTL_IDFGETLOC(IDF(1),IROW,ICOL,X,Y)
 !## get window for brush-function
 CALL UTL_IDFIROWICOL(IDF(1),IR2,IC1,X-TRAD,Y-TRAD)
 CALL UTL_IDFIROWICOL(IDF(1),IR1,IC2,X+TRAD,Y+TRAD)

 !## store conductance,stage,bottom,entry-restistance
 C=IDF(1)%X(ICOL,IROW)/A

 DO IR=IR1,IR2
  DO IC=IC1,IC2
    
   !## skip already filled in bathymetry/or previous wide profile
   IF(IDF(6)%X(IC,IR).LT.0.0.OR.IDF(6)%X(IC,IR).GE.1.0)CYCLE

   CALL IDFGETEDGE(IDF(1),IR,IC,X1,Y1,X2,Y2)
   XP(1)= X1        ; YP(1)=(Y1+Y2)/2.0
   XP(2)=(X1+X2)/2.0; YP(2)= Y1
   XP(3)= X2        ; YP(3)= YP(1)
   XP(4)= XP(2)     ; YP(4)= Y2
   XP(5)= X1        ; YP(5)= Y1  
   XP(6)= X1        ; YP(6)= Y2  
   XP(7)= X2        ; YP(7)= Y2  
   XP(8)= X2        ; YP(8)= Y1  
   XP(9)= XP(2)     ; YP(9)= YP(1)  

   K=0
   DO I=1,9
    DX=(XP(I)-X)**2.0
    DY=(YP(I)-Y)**2.0
    IF(SQRT(DX+DY).LT.TRAD)K=K+1
   ENDDO 

   IF(K.GT.0)THEN
    F=REAL(K)/9.0

    !## only overrule if fully covered or f larger than existing value
    IF(F.EQ.1.0.OR.F.GT.IDF(6)%X(IC,IR))THEN
     IDF(6)%X(IC,IR)=F

     !## multiply conductance for area
     IDF(1)%X(IC,IR)=C*F*IDFGETAREA(IDF(1),IC,IR)
     !## copy then all
     DO I=2,4; IDF(I)%X(IC,IR)=IDF(I)%X(ICOL,IROW); ENDDO

    ENDIF   

   ENDIF

  ENDDO
 ENDDO

 END SUBROUTINE ISG2GRID_EROSION_MATRIX