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

USE WINTERACTER
USE RESOURCE
USE MOD_SOBEK_PAR  
USE MOD_UTL, ONLY : UTL_WSELECTFILE,UTL_GETHELP,UTL_PUTRECORDLENGTH,ITIMETOGTIME
USE MOD_PREF_PAR, ONLY : PREFVAL
USE MOD_OSD, ONLY : OSD_OPEN
USE MOD_UTL, ONLY : ITOS,RTOS,UTL_GETUNIT,UTL_CAP,UTL_JDATETOIDATE,UTL_CREATEDIR,UTL_DIALOGSHOW,UTL_DIST
USE IMODVAR, ONLY : DP_KIND,SP_KIND,PI
USE MOD_QKSORT, ONLY : QKSORT
USE MOD_ISG_UTL, ONLY : UTL_GETUNITSISG
USE MOD_ISG_PAR, ONLY : ISGDOUBLE
USE IO_NETCDF
USE NETCDF
USE IO_UGRID
USE M_READCROSSSECTIONS
USE M_CROSSSECTIONS

INTEGER,PARAMETER :: ISOBEK=1
TYPE(T_UG_MESHGEOM) :: MESHGEOM
CHARACTER(LEN=UG_IDSLEN), ALLOCATABLE :: NBRANCHIDS(:), NNODEIDS(:), NODEIDS(:)       
CHARACTER(LEN=UG_IDSLONGNAMESLEN), ALLOCATABLE :: NBRANCHLONGNAMES(:), NNODELONGNAMES(:), NODELONGNAMES(:) 
TYPE(T_NETWORK) :: NETWORK
REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: S1
INTEGER,ALLOCATABLE,DIMENSION(:) :: IYMD,IHMS

CONTAINS

 !###======================================================================
 LOGICAL FUNCTION DFLOWFM1CALC()
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=13) :: DUMMY
 INTEGER :: IERR,IONCID,NMESH,IM,NETWORKINDEX,NCID,ID_S1,NT,ID_TIMEDIM,I,ID_TIME,TITLELENGTH,IYR,IMH,IDY,IHR,IMT,ISC,NDAY
! TYPE(T_UG_NETWORK) :: NETID 
 LOGICAL :: INCLUDEARRAYS !< (OPTIONAL) WHETHER OR NOT TO INCLUDE COORDINATE ARRAYS AND CONNECTIVITY TABLES. DEFAULT: .FALSE., I.E., DIMENSION COUNTS ONLY.
 INTEGER :: START_INDEX   !< (OPTIONAL) THE START INDEX   
 CHARACTER(LEN=256) :: NETWORK1DNAME, MESH1DNAME,REFDAT_MAP
! REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: MAPTIMES
 INTEGER,ALLOCATABLE,DIMENSION(:) :: IMAPTIMES
! INTEGER :: ID_NETNODEZ
 
 DFLOWFM1CALC=.FALSE.
 
 IERR=IONC_GETFULLVERSIONSTRING_IO_NETCDF(DUMMY); IF(IERR.NE.0)RETURN
 
 IERR=IONC_OPEN(TRIM(DFLOWFMDIR)//'\FLOWFM_NET.NC',NF90_NOWRITE,IONCID); IF(IERR.NE.0)RETURN
 
 IERR=IONC_GET_MESH_COUNT(IONCID,NMESH); IF(IERR.NE.0)RETURN
 
 IERR = IONC_GET_NCID(IONCID, NCID)
 
 INCLUDEARRAYS=.TRUE.
 START_INDEX=1
 DO IM=1,NMESH
  !## gathers mesh-info (dims, nodes, edges)
  IERR=IONC_GET_MESHGEOM(IONCID,IM,NETWORKINDEX,MESHGEOM); IF(IERR.NE.0)RETURN
  IF(MESHGEOM%DIM.EQ.1)THEN
         !Save meshgeom for later writing of the 1d network names
   !## retrieve the 1d geometry twice one time in meshgeom to process in this subroutine, one time in meshgeom1d, that can be used later during initialisation.
   IERR=IONC_GET_MESHGEOM(IONCID,IM,NETWORKINDEX,MESHGEOM,START_INDEX,INCLUDEARRAYS,NBRANCHIDS,NBRANCHLONGNAMES, &
             NNODEIDS,NNODELONGNAMES,NODEIDS,NODELONGNAMES,NETWORK1DNAME,MESH1DNAME); IF(IERR.NE.0)RETURN
   IERR=IONC_INQ_VARID_BY_STANDARD_NAME(IONCID,IM,UG_LOC_NODE,'SEA_SURFACE_HEIGHT',ID_S1)
   
   IERR = NF90_INQ_DIMID        (NCID, 'time'    , ID_TIMEDIM )
   IF(IERR.EQ.NF90_NOERR)THEN
    IERR = NF90_INQUIRE_DIMENSION(NCID, ID_TIMEDIM, LEN=NT) 
    IF(IERR.EQ.NF90_NOERR)THEN
     ALLOCATE(S1(MESHGEOM%NUMNODE,NT))
     IF(IERR.EQ.NF90_NOERR)THEN
      DO I=1,NT
       IERR=NF90_GET_VAR(NCID,ID_S1,S1(:,I),START=(/1,I/))
      ENDDO
     ENDIF
    
     IERR = NF90_INQ_VARID(NCID,'time', ID_TIME)
     IERR = NF90_INQUIRE_DIMENSION(NCID, ID_TIMEDIM, LEN=NT)
     IERR = NF90_INQUIRE_VARIABLE(NCID, ID_TIMEDIM, XTYPE=I)
     SELECT CASE(I)
      CASE(NF90_DOUBLE)
       WRITE(*,*) 'DOUBLE'
      CASE (NF90_FLOAT)
       WRITE(*,*) 'FLOAT'
      CASE (NF90_INT)
       WRITE(*,*) 'INT'
       ALLOCATE(iMAPTIMES(NT))
       IERR = NF90_GET_VAR(NCID, ID_TIME, iMAPTIMES)
     END SELECT
     IERR = NF90_INQUIRE_ATTRIBUTE(NCID, ID_TIME, 'units', LEN = TITLELENGTH)
     IERR = NF90_GET_ATT(NCID, ID_TIME, "units", REFDAT_MAP)
     READ(REFDAT_MAP(15:18),*) IYR
     READ(REFDAT_MAP(20:21),*) IMH
     READ(REFDAT_MAP(23:24),*) IDY
     READ(REFDAT_MAP(26:27),*) IHR
     READ(REFDAT_MAP(29:30),*) IMT
     READ(REFDAT_MAP(32:33),*) ISC
   
     ALLOCATE(IYMD(NT),IHMS(NT))
     DO I=1,NT
      ISC=ISC+IMAPTIMES(I)

      DO 
       IF(ISC.LT.60)EXIT
       IMT=IMT+1; ISC=ISC-60
      ENDDO
      DO 
       IF(IMT.LT.60)EXIT
       IHR=IHR+1; IMT=IMT-60
      ENDDO
      DO 
       IF(IHR.LT.24)EXIT
       IDY=IDY+1; IHR=IHR-24
      ENDDO
      DO 
       NDAY=WDATEDAYSINMONTH(IYR,IMH)
       IF(IDY.LE.NDAY)EXIT
       IMH=IMH+1; IDY=IDY-NDAY
      ENDDO
      DO 
       IF(IMH.LE.12)EXIT
       IMH=IMH-12; IYR=IYR+1
      ENDDO 
    
      IYMD(I)=IYR*10000+IMH*100+IDY
      IHMS(I)=IHR*3600+IMT*60 +ISC
     ENDDO
   
    ENDIF
   ELSE
    ISFINAL=1
   ENDIF

   CALL READCROSSSECTIONDEFINITIONS(NETWORK,TRIM(DFLOWFMDIR)//'\CRSDEF.INI')
   CALL READCROSSSECTIONLOCATIONFILE(NETWORK,TRIM(DFLOWFMDIR)//'\CRSLOC.INI')
    
   CALL OPEN_ISGFILES(1)
   CALL DFLOWFM_IMPORT1_MAIN() 
   CALL SOBEK1CLOSE()

  ENDIF
 ENDDO
 
!Deze file is zeer vergelijkbaar met de netwerk-file. Sterker nog, het is hetzelfde, maar dan uitgebreid met resultaatwaarden op het
!D- en het 2D-grid. Het grid staat er dus helemaal in. Dit zogenaamde output grid hoeft overigens niet altijd helemaal gelijk te zijn
!aan het input grid uit de _net.nd file, maar doorgaans wijkt het maar weinig af.
!In feite kun je dus beter de ISG’s genereren uit deze _map.nc file i.p.v. uit de _net.nc file, omdat er niet alleen het rooster in
 !staat, maar ook de waterstand.
!De variable mesh1d_s1 bevat de waterstanden op de 1D-gridpunten, in dezelfde volgorde als mesh1d_node_x / mesh1d_node_y /
 !mesh1d_node_id, die je waarschijnlijk al gebruikt in je huidige reader. De mesh1d_s1 variable heeft een extra dimensie (de
 !eerste dimensie), n.l. de tijd. Uit deze 2D-array kun je dus met één call de tijdreeks op een punt opvragen.
 
 IF(ALLOCATED(IYMD))DEALLOCATE(IYMD)
 IF(ALLOCATED(IHMS))DEALLOCATE(IHMS)
 IF(ALLOCATED(S1))  DEALLOCATE(S1)

 DFLOWFM1CALC=.TRUE.
 
 END FUNCTION DFLOWFM1CALC

 !###======================================================================
 SUBROUTINE DFLOWFM_IMPORT1_MAIN()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: ISEG,ICRS,IB,NNP,ICLC,NSEG,NCLC,NCRS,IREFSD,IREFSC,IREFST, &
            NSTW,ISTW,NQHR,IQHR,IREFSQ
 REAL(KIND=DP_KIND) :: DIST

 NNODE=MESHGEOM%NNODES
 NBRCH=MESHGEOM%NBRANCHES

 WRITE(IU(ISG),*) NBRCH,0

 IF(IBATCH.EQ.0)THEN; CALL WINDOWOUTSTATUSBAR(4,''); ELSE; WRITE(*,'(/A/)'); ENDIF
 
 ISEG  =0  !## segment-points
 ICRS  =0  !## cross-sections
 ICLC  =0  !## data
 ISTW  =0
 IQHR  =0
 IREFSD=0
 IREFSC=0
 IREFST=0
 IREFSQ=0
 NNP   =0

 DO IB=1,NBRCH
  CALL DFLOWFM_CREATEIPS(IB,ISEG,NSEG,DIST)
  CALL DFLOWFM_CREATEISD(IB,ICLC,NCLC,IREFSD,DIST)
  CALL DFLOWFM_CREATEISC(IB,ICRS,NCRS,IREFSC,DIST)
  NSTW=0
  NQHR=0

  LINE='"'//TRIM(NBRANCHIDS(IB))//'",'// &
       TRIM(ITOS(ISEG-NSEG+1))//','//TRIM(ITOS(NSEG))//','// &
       TRIM(ITOS(ICLC-NCLC+1))//','//TRIM(ITOS(NCLC))//','// &
       TRIM(ITOS(ICRS-NCRS+1))//','//TRIM(ITOS(NCRS))//','// &
       TRIM(ITOS(ISTW-NSTW+1))//','//TRIM(ITOS(NSTW))//','// &
       TRIM(ITOS(IQHR-NQHR+1))//','//TRIM(ITOS(NQHR))
  WRITE(IU(ISG),*) TRIM(LINE)

  IF(IBATCH.EQ.0)THEN
   CALL WINDOWOUTSTATUSBAR(4,'Progress '//TRIM(RTOS(REAL(IB*100,8)/REAL(NBRCH,8),'F',2))//' %')
  ELSE
   WRITE(6,'(A,F10.4,A)') '+Progress ',REAL(IB*100)/REAL(NBRCH),'%         '
  ENDIF
  
 END DO

 WRITE(IU(ISP) ,REC=1) UTL_PUTRECORDLENGTH(16) !8*256+247
 WRITE(IU(ISD1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(ISD2),REC=1) UTL_PUTRECORDLENGTH(44) !20*256+247
 WRITE(IU(ISC1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(ISC2),REC=1) UTL_PUTRECORDLENGTH(24) !12*256+247
 WRITE(IU(IST1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(IST2),REC=1) UTL_PUTRECORDLENGTH(28) !12*256+247
 WRITE(IU(ISQ1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(ISQ2),REC=1) UTL_PUTRECORDLENGTH(32) !16*256+247

 IF(IBATCH.NE.0)THEN
  WRITE(*,'(A)') 'Number of records in:'
  WRITE(*,'(A,I10)')  TRIM(FNAME(ISG))  //' ',NBRCH
  WRITE(*,'(A,I10)')  TRIM(FNAME(ISP))  //' ',NNP
  WRITE(*,'(A,2I10)') TRIM(FNAME(ISD1))//' ',ICLC,IREFSD
  WRITE(*,'(A,2I10)') TRIM(FNAME(ISC1))//' ',ICRS,IREFSC
  WRITE(*,'(A,2I10)') TRIM(FNAME(IST1))//' ',ISTW,IREFST
  WRITE(*,'(A,2I10)') TRIM(FNAME(ISQ1))//' ',IQHR,IREFSQ
 ENDIF
 
 END SUBROUTINE DFLOWFM_IMPORT1_MAIN
 
 !###======================================================================
 SUBROUTINE DFLOWFM_CREATEIPS(IB,ISEG,NSEG,DIST)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ISEG
 INTEGER,INTENT(OUT) :: NSEG
 REAL(KIND=DP_KIND),INTENT(OUT) :: DIST
 INTEGER :: I,N,IS
 
 IS=0; IF(IB.GT.1)IS=SUM(MESHGEOM%NBRANCHGEOMETRYNODES(1:IB-1))
 IS=IS+1
 N=MESHGEOM%NBRANCHGEOMETRYNODES(IB)

! OFFSET VAN DE NODES
! MESHGEOM%NODEOFFSETS(1) = 0.000000000000000D+000 _ CALC. POINTS

 NSEG=0; DIST=0.0D0
 DO I=IS,IS+N-1
  IF(I.GT.IS)THEN
   IF(MESHGEOM%NGEOPOINTX(I).NE.MESHGEOM%NGEOPOINTX(I-1).OR. &
      MESHGEOM%NGEOPOINTY(I).NE.MESHGEOM%NGEOPOINTX(I-1))THEN
    DIST=DIST+UTL_DIST(MESHGEOM%NGEOPOINTX(I),MESHGEOM%NGEOPOINTY(I),MESHGEOM%NGEOPOINTX(I-1),MESHGEOM%NGEOPOINTY(I-1))
    ISEG=ISEG+1; NSEG=NSEG+1
    WRITE(IU(ISP),REC=ISEG+1) MESHGEOM%NGEOPOINTX(I),MESHGEOM%NGEOPOINTY(I)
   ENDIF
  ELSE
   ISEG=ISEG+1; NSEG=NSEG+1
   WRITE(IU(ISP),REC=ISEG+1) MESHGEOM%NGEOPOINTX(I),MESHGEOM%NGEOPOINTY(I)
  ENDIF
 END DO

 IF(ABS(MESHGEOM%NBRANCHLENGTHS(IB)-DIST).GT.1.0D0)THEN
  WRITE(IU(IOUT),'(2(A,F15.7))') 'Actual length ',MESHGEOM%NBRANCHLENGTHS(IB),' not equal to computed length',DIST
 ENDIF
 
 END SUBROUTINE DFLOWFM_CREATEIPS

 !###======================================================================
 SUBROUTINE DFLOWFM_CREATEISD(IB,ICLC,NCLC,IREFSD,TDIST)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ICLC,IREFSD
 INTEGER,INTENT(OUT) :: NCLC
 REAL(KIND=DP_KIND),INTENT(IN) :: TDIST
 INTEGER :: I,J,JJ,II,III,N,ID1,ID2,IH,IM,IS
 REAL(KIND=DP_KIND) :: DIST,WP,BH,MIND1,MIND2,D,F1,F2
 CHARACTER(LEN=MAXLEN) :: CROSNAME
 CHARACTER(LEN=8) :: CTIME
 INTEGER,DIMENSION(2) :: IFT
 INTEGER,ALLOCATABLE,DIMENSION(:) :: IDS
 
 WP=0.0D0
 BH=0.0D0
 
 ALLOCATE(IDS(MESHGEOM%NUMNODE)); IDS=0
 
 NCLC=0; DO JJ=1,MESHGEOM%NUMEDGE
  IFT(1)=MESHGEOM%EDGE_NODES(1,JJ)!+1
  IFT(2)=MESHGEOM%EDGE_NODES(2,JJ)!+1

  !## is on current current branch
  IF(MESHGEOM%EDGEBRANCHIDX(JJ).NE.IB)CYCLE

  DO III=1,2
   I=IFT(III)
   !## already processed
   IF(IDS(I).EQ.1)CYCLE; IDS(I)=1
! !## add calculation nodes for current branch
! NCLC=0; DO I=1,MESHGEOM%NUMNODE
   !## not on current branch - set distance to zero
   IF(MESHGEOM%NODEBRANCHIDX(I).NE.IB)THEN
    DIST=0.0D0 !CYCLE
   ELSE
    DIST=MESHGEOM%NODEOFFSETS(I)
   ENDIF 
  
   DIST=MIN(DIST,TDIST)
  
   NCLC    =NCLC+1
   ICLC    =ICLC+1
   IF(ISFINAL.EQ.0)THEN
    N      =SIZE(S1,2)
   ELSE
    N      =1
   ENDIF
   CROSNAME=NODEIDS(I)

   WRITE(IU(ISD1),REC=ICLC+1) N,IREFSD+1,DIST,CROSNAME

   !## find bedlevel
   ID1=0; ID2=0; MIND1=HUGE(1.0); MIND2=MIND1; DO II=1,NETWORK%CRS%COUNT
    !## not on current branch
    IF(NETWORK%CRS%CROSS(II)%BRANCHID.NE.NBRANCHIDS(IB))CYCLE
    D=NETWORK%CRS%CROSS(II)%CHAINAGE
    D=D-DIST
    !## check closest before
    IF(D.LE.0.0D0.AND.ABS(D).LT.MIND1)THEN
     MIND1=ABS(D); ID1=II
    ENDIF
    !## check closest after
    IF(D.GE.0.0D0.AND.D.LT.MIND2)THEN
     MIND2=D; ID2=II
    ENDIF
   ENDDO
   
   IF(ID1.EQ.0)THEN
    BH=NETWORK%CRS%CROSS(ID2)%BEDLEVEL
   ELSEIF(ID2.EQ.0)THEN
    BH=NETWORK%CRS%CROSS(ID1)%BEDLEVEL
   ELSE
    F1=MIND2/(MIND1+MIND2)
    F2=1.0D0-F1
    BH=F1*NETWORK%CRS%CROSS(ID1)%BEDLEVEL+ &
       F2*NETWORK%CRS%CROSS(ID2)%BEDLEVEL
   ENDIF
 
   IF(ALLOCATED(S1))THEN
    IF(ISFINAL.EQ.0)THEN
     DO J=1,SIZE(S1,2)
      CALL ITIMETOGTIME(IHMS(J),IH,IM,IS)
      WP=S1(I,J)
      WP=MAX(WP,BH)
      IREFSD=IREFSD+1
      WRITE(CTIME,'(3(I2.2,A1))') IH,':',IM,':',IS
      WRITE(IU(ISD2),REC=IREFSD+1) IYMD(J),CTIME,WP,BH,1.0D0,0.3D0
     ENDDO
    ELSE
     J=SIZE(S1,2)
     CALL ITIMETOGTIME(IHMS(J),IH,IM,IS)
     WP=S1(I,J)
     WP=MAX(WP,BH)
     IREFSD=IREFSD+1
     WRITE(CTIME,'(3(I2.2,A1))') IH,':',IM,':',IS
     WRITE(IU(ISD2),REC=IREFSD+1) IYMD(J),CTIME,WP,BH,1.0D0,0.3D0
    ENDIF
   ELSE
    IREFSD=IREFSD+1
    WRITE(IU(ISD2),REC=IREFSD+1) 20200101,'00:00:00',BH+1.0D0,BH,1.0D0,0.3D0
   ENDIF
  ENDDO
 ENDDO
 
 END SUBROUTINE DFLOWFM_CREATEISD
  
 !###======================================================================
 SUBROUTINE DFLOWFM_CREATEISC(IB,ICRS,NCRS,IREFSC,TDIST)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ICRS,IREFSC
 INTEGER,INTENT(OUT) :: NCRS
 REAL(KIND=DP_KIND),INTENT(IN) :: TDIST
 INTEGER :: I,J,II,N
 REAL(KIND=DP_KIND) :: DIST,MRC,OR,DOR,RD,MY,XMID
 CHARACTER(LEN=MAXLEN) :: CROSNAME
 TYPE(T_CSTYPE),POINTER :: PCSDEF
! TYPE(T_CROSSSECTION),POINTER :: CRSDEF
 REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: XPROF,YPROF,XP,YP
 
 MRC=0.03D0
 
 NCRS=0; DO I=1,NETWORK%CRS%COUNT
  !## not on current branch
  IF(NETWORK%CRS%CROSS(I)%BRANCHID.NE.NBRANCHIDS(IB))CYCLE
  DIST=NETWORK%CRS%CROSS(I)%CHAINAGE

  !## get correct cross-section definition
  DO II=1,NETWORK%CSDEFINITIONS%COUNT
   PCSDEF=>NETWORK%CSDEFINITIONS%CS(II)
   IF(NETWORK%CRS%CROSS(I)%TABDEF%ID.EQ.PCSDEF%ID)THEN
    CROSNAME=NETWORK%CRS%CROSS(I)%TABDEF%ID
    EXIT
   ENDIF
  ENDDO  
  
  NCRS    =NCRS+1
  ICRS    =ICRS+1
!  CROSNAME=NETWORK%CRS%CROSS(I)%CSID 
!  CRSDEF=>NETWORK%CRS%CROSS(I)
  SELECT CASE(PCSDEF%CROSSTYPE)
   !## yz table, z is rising but y can decrease which cannot in iMOD
   CASE (CS_TABULATED)

    N=PCSDEF%LEVELSCOUNT
    ALLOCATE(XPROF(N*2),YPROF(N*2),XP(N),YP(N))
    XPROF(1:N)=PCSDEF%FLOWWIDTH(1:N)
    YPROF(1:N)=PCSDEF%HEIGHT(1:N)

    XP(1)=XPROF(1)
    YP(1)=YPROF(1)
    J=1; DO II=2,N
     !## skip as cross-section goes the opposite direction
     IF(XPROF(II).GE.XPROF(J))THEN
      J=J+1
      XP(J)=XPROF(II)
      YP(J)=YPROF(II)
!     ELSE
!      WRITE(*,*) 
     ENDIF
    END DO
    
    N=J
    
    !## divide by 2
    XP=XP/2.0D0

    II=0; DO J=1,N
     II=II+1; XPROF(II)=-XP(N-J+1); YPROF(II)=YP(N-J+1)
    ENDDO
    DO J=1,N
     II=II+1; XPROF(II)= XP(J); YPROF(II)=YP(J)
    ENDDO
    
    N=II
    DEALLOCATE(XP,YP)
    
   !## closed circle
   CASE (CS_CIRCLE)
    
    N=11; ALLOCATE(XPROF(N),YPROF(N))
    OR =0.5D0*(2.0D0*PI)
    DOR=0.05D0*(2.0D0*PI)
    OR=OR+DOR
    DO II=1,N
     OR=OR-DOR
     XPROF(II)=0.0D0+RD*COS(OR)
     YPROF(II)=      RD-RD*SIN(OR)
    ENDDO
     
   CASE (CS_YZ_PROF)
    N=PCSDEF%LEVELSCOUNT
    ALLOCATE(XPROF(N),YPROF(N))
    DO J=1,N
     XPROF(J)=PCSDEF%Y(J)
     YPROF(J)=PCSDEF%Z(J)
    END DO 
    CALL QKSORT(N,XPROF,V2=YPROF)
  END SELECT

  !## set on zero as base
  MY=MINVAL(YPROF(1:N)); DO II=1,N; YPROF(II)=YPROF(II)-MY; ENDDO
  !## shift such that centre is mid
  XMID=(XPROF(1)+XPROF(N))/2.0D0
  !## set new x-values
  XPROF=XPROF-XMID
  
  WRITE(IU(ISC1),REC=ICRS+1) N,IREFSC+1,DIST,CROSNAME
  DO J=1,N
   IREFSC=IREFSC+1
   WRITE(IU(ISC2),REC=IREFSC+1) XPROF(J),YPROF(J),MRC
  END DO 

  DEALLOCATE(XPROF,YPROF)

 ENDDO
 
 END SUBROUTINE DFLOWFM_CREATEISC
  
 !###======================================================================
 SUBROUTINE SOBEK1MAIN()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: ITYPE
 TYPE(WIN_MESSAGE) :: MESSAGE
 CHARACTER(LEN=256) :: FNAME
 
 CALL WDIALOGLOAD(ID_DIMPORTSOBEK,ID_DIMPORTSOBEK)

 CALL WDIALOGPUTIMAGE(ID_OPEN1,ID_ICONSAVEAS)
 CALL WDIALOGPUTIMAGE(ID_OPEN2,ID_ICONOPEN)
 CALL WDIALOGPUTIMAGE(ID_OPEN3,ID_ICONOPEN)
 CALL WDIALOGPUTIMAGE(ID_OPEN4,ID_ICONOPEN)  

 CALL WDIALOGSELECT(ID_DIMPORTSOBEK)
 CALL UTL_DIALOGSHOW(-1,-1,0,3)
 IBATCH=0
 
 DO
  CALL WMESSAGE(ITYPE,MESSAGE)
  SELECT CASE (ITYPE)

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

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

   CASE(PUSHBUTTON)
    SELECT CASE (MESSAGE%VALUE1)
     CASE (ID_OPEN1)
      FNAME='' 
      IF(UTL_WSELECTFILE('iMOD ISG file|*.isg|', &
                   SAVEDIALOG+PROMPTON+DIRCHANGE+APPENDEXT,FNAME,&
                   'Save iMOD ISG File (*.isg)'))THEN
       CALL WDIALOGPUTSTRING(IDF_STRING1,FNAME)
      ENDIF
     CASE (ID_OPEN2)
      FNAME='' 
      IF(UTL_WSELECTFILE('Sobek Network File|*.tp|', &
                   LOADDIALOG+MUSTEXIST+PROMPTON+DIRCHANGE+APPENDEXT,FNAME,&
                   'Load Sobek Network File (network.tp)'))THEN
       CALL WDIALOGPUTSTRING(IDF_STRING2,FNAME)
      ENDIF
     CASE (ID_OPEN3)
      FNAME='' 
      IF(UTL_WSELECTFILE('SOBEK Calcpnt His-file|*.his|', &
                   LOADDIALOG+MUSTEXIST+PROMPTON+DIRCHANGE+APPENDEXT,FNAME,&
                   'Load SOBEK His file (*.his)'))THEN
       CALL WDIALOGPUTSTRING(IDF_STRING3,FNAME)
      ENDIF
     CASE (ID_OPEN4)
      FNAME='' 
      IF(UTL_WSELECTFILE('SOBEK Struc His-file|*.his|', &
                   LOADDIALOG+MUSTEXIST+PROMPTON+DIRCHANGE+APPENDEXT,FNAME,&
                   'Load SOBEK His file (*.his)'))THEN
       CALL WDIALOGPUTSTRING(IDF_STRING4,FNAME)
      ENDIF
     CASE (IDOK)
      CALL SOBEK1FIELDS()
      IF(SOBEK1CALC())THEN; ENDIF !EXIT
     CASE (IDHELP)
      CALL UTL_GETHELP('5.5.1','TMO.IT.SOBEK')
     !## cancel modeling
     CASE (IDCANCEL)
      EXIT
    END SELECT

  END SELECT
 ENDDO
 
 CALL WDIALOGSELECT(ID_DIMPORTSOBEK)
 CALL WDIALOGUNLOAD()

 END SUBROUTINE SOBEK1MAIN

 !###======================================================================
 SUBROUTINE SOBEK1FIELDS
 !###======================================================================
 IMPLICIT NONE
 
 CALL WDIALOGGETSTRING(IDF_STRING1,ISGNAME)
 CALL WDIALOGGETSTRING(IDF_STRING2,SOBEKDIR)
 CALL WDIALOGGETSTRING(IDF_STRING3,CALCPNTHISNAME)
 CALL WDIALOGGETSTRING(IDF_STRING4,STRUCHISNAME)

 END SUBROUTINE SOBEK1FIELDS

 !###======================================================================
 LOGICAL FUNCTION SOBEK1CALC()
 !###======================================================================
 IMPLICIT NONE
 
 SOBEK1CALC=.FALSE.
 
 !## input
 FNAME(ITP)= TRIM(SOBEKDIR)//'\NETWORK.TP'
 FNAME(ICR)= TRIM(SOBEKDIR)//'\NETWORK.CR'
 FNAME(ICP)= TRIM(SOBEKDIR)//'\NETWORK.CP'
 FNAME(IGR)= TRIM(SOBEKDIR)//'\NETWORK.GR'
 FNAME(IST)= TRIM(SOBEKDIR)//'\NETWORK.ST'
 FNAME(IAT)= TRIM(SOBEKDIR)//'\PROFILE.DAT'
 FNAME(IEF)= TRIM(SOBEKDIR)//'\PROFILE.DEF'
 FNAME(IFR)= TRIM(SOBEKDIR)//'\FRICTION.DAT'
 FNAME(IHIS)=CALCPNTHISNAME 
 FNAME(SHIS)=STRUCHISNAME   
 
 CALL OPEN_ISGFILES(0)

 !## get his in data set
 IOS(IHIS)=0
 IF(IBATCH.EQ.0)THEN
  CALL WINDOWOUTSTATUSBAR(4,'Reading '//TRIM(FNAME(IHIS))//'...')
 ELSE
  WRITE(*,'(/1X,A)') 'Reading '//TRIM(FNAME(IHIS))//'...'
 ENDIF
 HISDATASET=DIOPLTGETDATASET(FNAME(IHIS))
 IF(.NOT.DIOPLTOPENEDOK(HISDATASET))THEN
  CALL DIOPLTDESTROY(HISDATASET)
  IOS(IHIS)=-1
 ENDIF
 IF(FNAME(SHIS).NE.'')THEN
  IF(IBATCH.EQ.0)THEN
   CALL WINDOWOUTSTATUSBAR(4,'Reading '//TRIM(FNAME(SHIS))//'...')
  ELSE
   WRITE(*,'(/1X,A)') 'Reading '//TRIM(FNAME(SHIS))//'...'
  ENDIF
  SHISDATASET=DIOPLTGETDATASET(FNAME(SHIS))
  IF(.NOT.DIOPLTOPENEDOK(SHISDATASET))THEN
   CALL DIOPLTDESTROY(SHISDATASET)
   IOS(SHIS)=-1
  ENDIF
 ENDIF
 
 DIMTIMES=0

 IF(SUM(IOS).EQ.0)THEN
  
  WRITE(IU(IOUT),'(A)') ' Reads:   '//TRIM(FNAME(ITP)) !NETWORK.TP'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ICR)) !NETWORK.CR'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ICP)) !NETWORK.CP'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IGR)) !NETWORK.GR'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IST)) !NETWORK.ST'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IAT)) !PROFILE.DAT'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IEF)) !PROFILE.DEF'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IFR)) !FRICTION.DAT'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IHIS))!CLCPNT.HIS'
  WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(SHIS))!STRUC.HIS'

  CALL SOBEK1IMPORT_MAIN()
 
 ENDIF
 
 CALL SOBEK1CLOSE()
  
 SOBEK1CALC=.TRUE.

 END FUNCTION SOBEK1CALC

 !###======================================================================
 SUBROUTINE SOBEK1CLOSE()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I
 LOGICAL :: LEX
 
 IF(IBATCH.NE.0)WRITE(*,*)
 DO I=1,NIU
  IF(IOS(I).EQ.0)THEN
   INQUIRE(UNIT=IU(I),OPENED=LEX)
   IF(LEX)CLOSE(IU(I))
  ELSE
   IF(IBATCH.EQ.0)THEN
    CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot open/find: '//TRIM(FNAME(I)),'Error')
   ELSE
    WRITE(*,'(1X,A)') 'Cannot open/find: '//TRIM(FNAME(I))
   ENDIF
  ENDIF
 END DO
 !## destroy dataset
 IF(IOS(IHIS).EQ.0)CALL DIOPLTDESTROY(HISDATASET)
 
 CLOSE(IU(IOUT))
 
 IF(SUM(IOS).EQ.0)THEN
  IF(IBATCH.EQ.0)THEN
   CALL WMESSAGEBOX(OKONLY,INFORMATIONICON,COMMONOK,'Sobek Network and results from selected *.HIS file(s):'//CHAR(13)// &
    TRIM(FNAME(IHIS))//CHAR(13)//TRIM(FNAME(SHIS))//CHAR(13)//'Successfully imported to iMOD-ISG file:'//CHAR(13)// &
    TRIM(FNAME(ISG)),'Information')
  ELSE
   WRITE(*,'(/1X,A/)') ' Program Successfully completed '
  ENDIF
 ENDIF
 
 END SUBROUTINE SOBEK1CLOSE

 !###======================================================================
 SUBROUTINE OPEN_ISGFILES(IOPTION)
 !###======================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IOPTION
 INTEGER :: I,J
 INTEGER,ALLOCATABLE,DIMENSION(:) :: JU
 LOGICAL :: LEX
 
 !## create folder
 I=INDEX(ISGNAME,'\',.TRUE.)-1
 CALL UTL_CREATEDIR(ISGNAME(:I))

 FNAME(IOUT)=ISGNAME(:I)//'\imod2sobek.log'  
 IU(IOUT)=UTL_GETUNIT(); CALL OSD_OPEN(IU(IOUT),FILE=FNAME(IOUT),STATUS='UNKNOWN',FORM='FORMATTED',IOSTAT=IOS(IOUT))

 I=INDEX(ISGNAME,'.',.TRUE.)-1
 IF(I.EQ.-1)I=LEN_TRIM(ISGNAME) !## no extention found
  
 FNAME(ISG) =ISGNAME(:I)//'.ISG'
 FNAME(ISP) =ISGNAME(:I)//'.ISP'
 FNAME(ISD1)=ISGNAME(:I)//'.ISD1'
 FNAME(ISD2)=ISGNAME(:I)//'.ISD2'
 FNAME(ISC1)=ISGNAME(:I)//'.ISC1'
 FNAME(ISC2)=ISGNAME(:I)//'.ISC2'
 FNAME(IST1)=ISGNAME(:I)//'.IST1'
 FNAME(IST2)=ISGNAME(:I)//'.IST2'
 FNAME(ISQ1)=ISGNAME(:I)//'.ISQ1'
 FNAME(ISQ2)=ISGNAME(:I)//'.ISQ2'
  
 IOS=0

 !## read sobek-files
 IF(IOPTION.EQ.0)THEN
  DO J=1,8
   I=IL(J)
   IU(I)=UTL_GETUNIT()
   CALL OSD_OPEN(IU(I),FILE=FNAME(I),STATUS='OLD',ACTION='READ',IOSTAT=IOS(I))
   !## error opening file
   IF(IOS(I).NE.0)THEN
    IF(IBATCH.EQ.0)THEN
     CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot find:'//CHAR(13)//TRIM(FNAME(I)),'Error')
    ELSE
     WRITE(*,'(A)') 'Cannot find:'//CHAR(13)//TRIM(FNAME(I))
    ENDIF
    DO I=1,NIU
     IF(IU(I).GT.0)THEN; INQUIRE(UNIT=IU(I),OPENED=LEX); IF(LEX)CLOSE(IU(I)); ENDIF
    ENDDO
    RETURN
   ENDIF
  END DO
 ENDIF

 !## always create a double precision ISG
 ISGDOUBLE=8
 IF(ALLOCATED(JU))DEALLOCATE(JU); ALLOCATE(JU(10)); JU=0
 CALL UTL_GETUNITSISG(JU,FNAME(ISG),'REPLACE')
 !## something went wrong
 IF(SUM(JU).EQ.0)RETURN
 IU(ISG) =JU(1)
 IU(ISP) =JU(2)
 IU(ISD1)=JU(3)
 IU(ISD2)=JU(4)
 IU(ISC1)=JU(5)
 IU(ISC2)=JU(6)
 IU(IST1)=JU(7)
 IU(IST2)=JU(8)
 IU(ISQ1)=JU(9)
 IU(ISQ2)=JU(10)
 DEALLOCATE(JU)

 WRITE(IU(IOUT),'(A)') ' Creates: '//TRIM(FNAME(ISG)) !IMOD.ISG'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISP)) !IMOD.ISP'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISD1))!IMOD.ISD1'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISD2))!IMOD.ISD2'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISC1))!IMOD.ISC1'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISC2))!IMOD.ISC2'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IST1))!IMOD.IST1'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(IST2))!IMOD.IST2'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISQ1))!IMOD.ISQ1'
 WRITE(IU(IOUT),'(A)') '          '//TRIM(FNAME(ISQ2))!IMOD.ISQ2'
 WRITE(IU(IOUT),'(A/)') '          '//TRIM(FNAME(IOUT))!LOG-FILE'

 END SUBROUTINE OPEN_ISGFILES

 !###======================================================================
 SUBROUTINE SOBEK1IMPORT_MAIN()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: I
 
 !## allocate memory for SOBEK-system
 IF(IBATCH.EQ.0)THEN
  CALL WINDOWOUTSTATUSBAR(4,'Allocated ...')
 ELSE
  WRITE(*,'(1X,A)') 'Allocated ...'
 ENDIF

 CALL MAIN1ALNETWORKTP()
 CALL MAIN1ALNETWORKCP()
 CALL MAIN1ALNETWORKCR()
 CALL MAIN1ALNETWORKGR()
 CALL MAIN1ALNETWORKST()
 CALL MAIN1ALPROFILEDAT()
 CALL MAIN1ALPROFILEDEF()
 CALL MAIN1ALFRICTIONDAT()

 ALLOCATE(TP1(NNODE),TP2(NBRCH),TP3(NNDLK),CP(NCPNT),CR(NCROS),GR(NCALCP), &
          PDAT(NDAT),PDEF(NDEF),BDFR(NBDFR),ST(NSTUW),QH(NCALCP))
 DO I=1,SIZE(PDEF); NULLIFY(PDEF(I)%XPROF,PDEF(I)%YPROF); ENDDO
 
 IF(IBATCH.EQ.0)THEN
  CALL WINDOWOUTSTATUSBAR(4,'Allocated ...')
 ELSE
  WRITE(*,'(1X,A)') 'Read ...'
 ENDIF
 
 CALL MAIN1RDNETWORKTP()
 CALL MAIN1RDNETWORKCP()
 CALL MAIN1RDNETWORKCR()
 CALL MAIN1RDNETWORKGR()
 CALL MAIN1RDNETWORKST()
 CALL MAIN1RDPROFILEDAT()
 CALL MAIN1RDPROFILEDEF()
 CALL MAIN1RDFRICTIONDAT()

 CALL SOBEK_CREATEISG()

 DO I=1,SIZE(PDEF)
  IF(ASSOCIATED(PDEF(I)%XPROF))DEALLOCATE(PDEF(I)%XPROF)
  IF(ASSOCIATED(PDEF(I)%YPROF))DEALLOCATE(PDEF(I)%YPROF)
 ENDDO
 
 DEALLOCATE(TP1,TP2,CP,CR,GR,PDAT,PDEF,BDFR,ST,QH)

 END SUBROUTINE SOBEK1IMPORT_MAIN

 !##=====================================================================
 SUBROUTINE SOBEK_CREATEISG()
 !##=====================================================================
 IMPLICIT NONE
 INTEGER :: ISEG,ICRS,IB,NNP,ICLC,NSEG,NCLC,NCRS,IREFSD,IREFSC,IREFST, &
            NSTW,ISTW,NQHR,IQHR,IREFSQ

 ALLOCATE(X(MAXCPNT),Y(MAXCPNT))

 WRITE(IU(ISG),*) NBRCH,0

 IF(IBATCH.EQ.0)THEN
  CALL WINDOWOUTSTATUSBAR(4,'')
 ELSE
  WRITE(*,'(/A/)')
 ENDIF
 
 ISEG  =0  !## segment-points
 ICRS  =0  !## cross-sections
 ICLC  =0  !## data
 ISTW  =0
 IQHR  =0
 IREFSD=0
 IREFSC=0
 IREFST=0
 IREFSQ=0
 NNP   =0
 DO IB=1,NBRCH

  IF(.NOT.SOBEK_CREATEISP(IB,ISEG,NNP,NSEG))EXIT
  CALL SOBEK_CREATEISD1(ICLC,NCLC,IREFSD,IB)
  CALL SOBEK_CREATEISC1(ICRS,IB,NCRS,IREFSC)
  CALL SOBEK_IST1CREATE(ISTW,IB,NSTW,IREFST)
  CALL SOBEK_ISQ1CREATE(IQHR,IB,NQHR,IREFSQ)

  LINE='"'//TRIM(TP2(IB)%ID)//'",'//TRIM(ITOS(ISEG-NSEG+1))//','//TRIM(ITOS(NSEG))//','// &
                                    TRIM(ITOS(ICLC-NCLC+1))//','//TRIM(ITOS(NCLC))//','// &
                                    TRIM(ITOS(ICRS-NCRS+1))//','//TRIM(ITOS(NCRS))//','// &
                                    TRIM(ITOS(ISTW-NSTW+1))//','//TRIM(ITOS(NSTW))//','// &
                                    TRIM(ITOS(IQHR-NQHR+1))//','//TRIM(ITOS(NQHR))
  WRITE(IU(ISG),*) TRIM(LINE)

  IF(IBATCH.EQ.0)THEN
   CALL WINDOWOUTSTATUSBAR(4,'Progress '//TRIM(RTOS(REAL(IB*100,8)/REAL(NBRCH,8),'F',2))//' %')
  ELSE
   WRITE(6,'(A,F10.4,A)') '+Progress ',REAL(IB*100)/REAL(NBRCH),'%         '
  ENDIF
  
 END DO

 !0,16,48,44,48,24,48,28,48,32
 
 WRITE(IU(ISP) ,REC=1) UTL_PUTRECORDLENGTH(16) !8*256+247
 WRITE(IU(ISD1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(ISD2),REC=1) UTL_PUTRECORDLENGTH(44) !20*256+247
 WRITE(IU(ISC1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(ISC2),REC=1) UTL_PUTRECORDLENGTH(24) !12*256+247
 WRITE(IU(IST1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(IST2),REC=1) UTL_PUTRECORDLENGTH(28) !12*256+247
 WRITE(IU(ISQ1),REC=1) UTL_PUTRECORDLENGTH(48) !44*256+247 
 WRITE(IU(ISQ2),REC=1) UTL_PUTRECORDLENGTH(32) !16*256+247

 IF(IBATCH.NE.0)THEN
  WRITE(*,'(A)') 'Number of records in:'
  WRITE(*,'(A,I10)') TRIM(FNAME(ISG))  //' ',NBRCH
  WRITE(*,'(A,I10)') TRIM(FNAME(ISP))  //' ',NNP
  WRITE(*,'(A,2I10)') TRIM(FNAME(ISD1))//' ',ICLC,IREFSD
  WRITE(*,'(A,2I10)') TRIM(FNAME(ISC1))//' ',ICRS,IREFSC
  WRITE(*,'(A,2I10)') TRIM(FNAME(IST1))//' ',ISTW,IREFST
  WRITE(*,'(A,2I10)') TRIM(FNAME(ISQ1))//' ',IQHR,IREFSQ
 ENDIF
 
 IF(ALLOCATED(X))DEALLOCATE(X)
 IF(ALLOCATED(Y))DEALLOCATE(Y)
 IF(ALLOCATED(TIMINDICES))DEALLOCATE(TIMINDICES)
 IF(ALLOCATED(JULIANTIMES))DEALLOCATE(JULIANTIMES)
 IF(ALLOCATED(SELECTEDVALUES))DEALLOCATE(SELECTEDVALUES)

 END SUBROUTINE SOBEK_CREATEISG

 !##=====================================================================
 LOGICAL FUNCTION SOBEK_CREATEISP(IB,ISEG,NNP,NSEG)
 !##=====================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),PARAMETER :: RAD=360.0D0/(2.0D0*3.1415926)
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ISEG,NNP
 INTEGER,INTENT(OUT) :: NSEG
 REAL(KIND=DP_KIND) :: DIST,D
 REAL(KIND=DP_KIND) :: XC,YC
 INTEGER :: I,J,K

 SOBEK_CREATEISP=.FALSE.
 
 !## begin node
 J=1
 
 DO I=1,NNODE
  IF(TRIM(TP1(I)%ID).EQ.TRIM(TP2(IB)%CBN))EXIT 
 END DO
 IF(I.LE.NNODE)THEN
  XC=TP1(I)%PX; YC=TP1(I)%PY
 ELSE 
  !## try from flow-linkage nodes, if available
  DO I=1,NNDLK
   IF(TRIM(TP3(I)%ID).EQ.TRIM(TP2(IB)%CBN))EXIT 
  END DO
  IF(I.LE.NNDLK)THEN
   XC=TP3(I)%PX; YC=TP3(I)%PY
  ELSE
   IF(IBATCH.EQ.0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot match identification ('//TRIM(TP2(IB)%CBN)// &
    ') for branch: '//TRIM(ITOS(IB)),'Error')
   IF(IBATCH.EQ.1)WRITE(*,'(A,I10)') 'Cannot match identification ('//TRIM(TP2(IB)%CBN)//') for branch: ',IB
   RETURN
  ENDIF
 ENDIF
 X(J)=XC; Y(J)=YC

 !## put in between curving points
 D=0.0D0
 DO I=1,NCPNT
 !## curving points id equal to current branch-id
  IF(TRIM(CP(I)%ID).EQ.TRIM(TP2(IB)%ID))THEN

   !## angle->radians
   DIST=CP(I)%SLEN-D
   D   =D+DIST
   XC  =XC+SIN(CP(I)%ANGLE/RAD)*DIST
   YC  =YC+COS(CP(I)%ANGLE/RAD)*DIST
   J   =J+1
   IF(J.GT.MAXCPNT)THEN
    STOP 'increase MAXCPNT'
   ENDIF
   X(J)=XC
   Y(J)=YC

  ENDIF
 END DO

 !## skip double points
 K=1
 DO I=2,J
  IF(X(I).NE.X(I-1).OR.Y(I).NE.Y(I-1))THEN
   K   =K+1
   X(K)=X(I)
   Y(K)=Y(I)
  ENDIF
 END DO

 J=K

 DIST=0.0D0
 DO I=1,J
  IF(I.GT.1)THEN
   IF(X(I).NE.X(I-1).OR.Y(I).NE.Y(I-1))THEN
    DIST=DIST+SQRT((X(I)-X(I-1))**2.0D0+(Y(I)-Y(I-1))**2.0D0)
    ISEG=ISEG+1
    WRITE(IU(ISP),REC=ISEG+1) X(I),Y(I)
   ENDIF
  ELSE
   ISEG=ISEG+1
   WRITE(IU(ISP),REC=ISEG+1) X(I),Y(I)
  ENDIF
 END DO

 IF(ABS(TP2(IB)%AL-DIST).GT.1.0D0)THEN
  WRITE(IU(IOUT),'(2(A,F15.7))') 'Actual length ',TP2(IB)%AL,' not equal to computed length',DIST
 ENDIF
 
 !## adjust correct segment length
 TP2(IB)%AL=DIST

 NSEG=J
 NNP =NNP+J

 SOBEK_CREATEISP=.TRUE.

 END FUNCTION SOBEK_CREATEISP

 !##=====================================================================
 SUBROUTINE SOBEK_CREATEISD1(ICLC,NCLC,IREFSD,IB)
 !##=====================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),PARAMETER :: NODATA=-999.99D0
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ICLC,IREFSD
 INTEGER,INTENT(OUT) :: NCLC
 INTEGER :: ID,I,J,N
 CHARACTER(LEN=MAXLEN) :: CROSNAME

 !##from zomer/winter
 !CALL ISD1GETDATA('calcpnt.HIS')

 !ID=20040101

 !##from
 !ICLC  =ICLC+1
 !IREFSD=IREFSD+1
 !WRITE(IU(ISD1),REC=ICLC+1) 1,IREFSD,0.0D0,'CalcPnt_FR'
 !WRITE(IU(ISD2),REC=IREFSD+1) ID,WL(1),WB(1),RS(1),FC(1)

 !##to
 !ICLC  =ICLC+1
 !IREFSD=IREFSD+1
 !WRITE(IU(ISD1),REC=ICLC+1) 1,IREFSD,REAL(TP2(IB)%AL),'CalcPnt_TO'
 !WRITE(IU(ISD2),REC=IREFSD+1) ID,WL(3),WB(3),RS(3),FC(3)

 !NCLC  =2
 !RETURN

 !10 CONTINUE

 NCLC=0
 DO I=1,NCALCP
  IF(TRIM(TP2(IB)%ID).EQ.TRIM(GR(I)%CI))THEN

   !## from his-file
   CALL HIS1GETDATA(GR(I)%ID,1)

   IF(NTIMES.GT.0)THEN

    NCLC    =NCLC+1
    ICLC    =ICLC+1

    N       =NTIMES
    CROSNAME=GR(I)%ID

    WRITE(IU(ISD1),REC=ICLC+1) N,IREFSD+1,GR(I)%LC,CROSNAME
    DO J=1,N
     IREFSD=IREFSD+1
     ID    =UTL_JDATETOIDATE(INT(JULIANTIMES(J)))
     WRITE(IU(ISD2),REC=IREFSD+1) ID,'00:00:00',SELECTEDVALUES(1,1,J),SELECTEDVALUES(2,1,J),1.0D0,0.3D0
    END DO

   ENDIF
  ENDIF
 ENDDO

 !IF(NCLC.EQ.0)WRITE(*,*) TRIM(TP2(IB)%ID)//' not found'

 !#check whether last calculationpoint lenth is equal to actual length of segment
 !READ(IU(ISD1),REC=ICLC+1)  N,IREFSD,LENGTH,CROSNAME
 !WRITE(IU(ISD1),REC=ICLC+1) N,IREFSD,REAL(TP2(IB)%AL),CROSNAME

 END SUBROUTINE SOBEK_CREATEISD1

 !##=====================================================================
 SUBROUTINE HIS1GETDATA(CLCID,ITYPE)
 !##=====================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: CLCID
 INTEGER,INTENT(IN) :: ITYPE
 INTEGER :: NPAR,NLOC
 CHARACTER(LEN=DIOMAXPARLEN),DIMENSION(:),POINTER :: PARS
 CHARACTER(LEN=DIOMAXLOCLEN),DIMENSION(:),POINTER :: LOCS
 INTEGER,DIMENSION(2) :: WATERLEVELINDICES
 INTEGER,DIMENSION(1) :: SELECTEDLOCATIONINDICES
 INTEGER :: I
 LOGICAL :: LWP,LWB
  
 !QUANTITYID(1,1)='Waterlevel  (m AD)'
 !QUANTITYID(2,1)='Waterdepth  (m)'
 !QUANTITYID(1,2)='Discharge (m|/s)'
 !QUANTITYID(2,2)='Waterlevel up (m AD)'
 !QUANTITYID(3,2)='Waterlevel down (m AD)'
 !QUANTITYID(4,2)='Crest level (m AD)'
 !QUANTITYID(5,2)='Gate lower edge level (m AD)'

 ! find the index of the requested quantity
 IF(ITYPE.EQ.1)THEN
  NPAR= DIOPLTGETNPAR(HISDATASET)
  PARS=>DIOPLTGETPARS(HISDATASET)
  LWP=.FALSE.; LWB=.FALSE.
  IF(TRIM(PARS(1)).EQ.'Waterlevel  (m AD)'.OR. &
     TRIM(PARS(1)).EQ.'Waterlevel mean (m AD)'.OR. &
     TRIM(PARS(1)).EQ.'Waterlevel max. (m AD)'.OR. &
     TRIM(PARS(1)).EQ.'Waterlevel mean (m A')LWP=.TRUE.
  IF(NPAR.GT.1)THEN
   IF(PARS(2).EQ.'Waterdepth  (m)'.OR. &
      PARS(2).EQ.'Waterdepth max. (m)'.OR. &
      PARS(2).EQ.'Waterdepth mean (m)')LWB=.TRUE.
  ENDIF
  IF(.NOT.LWP)THEN
   WRITE(IU(IOUT),'(A)') 'HIS-file does not contain the appropriate quantities:'
   WRITE(IU(IOUT),'(A)') '- Waterlevel  (m AD)'
   WRITE(IU(IOUT),'(A)') '  or'
   WRITE(IU(IOUT),'(A)') '- Waterlevel mean (m A)'
   IF(.NOT.LWB)WRITE(IU(IOUT),'(A)') '- Waterdepth  (m) <--- optinal'
   WRITE(IU(IOUT),'(A)') 'Available parameters are:'
   DO I=1,SIZE(PARS); WRITE(IU(IOUT),'(I5,A)') I,TRIM(PARS(I)); ENDDO
   IF(IBATCH.EQ.1)STOP 'error occured check tmp\imod2sobek.log'
   RETURN
  ENDIF
  WATERLEVELINDICES=1
  IF(LWB)WATERLEVELINDICES(2)=2
 ELSEIF(ITYPE.EQ.2)THEN
  NPAR= DIOPLTGETNPAR(SHISDATASET)
  PARS=>DIOPLTGETPARS(SHISDATASET)
  WATERLEVELINDICES(1)=2
  WATERLEVELINDICES(2)=3
 ENDIF

 IF(ITYPE.EQ.1)THEN
  NLOC= DIOPLTGETNLOC(HISDATASET)
  LOCS=>DIOPLTGETLOCS(HISDATASET)
 ELSEIF(ITYPE.EQ.2)THEN
  NLOC= DIOPLTGETNLOC(SHISDATASET)
  LOCS=>DIOPLTGETLOCS(SHISDATASET)
 ENDIF

 SELECTEDLOCATIONINDICES=0
 DO I=1,NLOC
  IF(TRIM(UTL_CAP(LOCS(I),'U')).EQ.TRIM(UTL_CAP(CLCID,'U')))SELECTEDLOCATIONINDICES(1)=I
 ENDDO
 IF(SELECTEDLOCATIONINDICES(1).EQ.0)THEN
  IF(ITYPE.EQ.1)WRITE(IU(IOUT),'(A)') 'CALCULATION POINT '//TRIM(CLCID)//' NOT FOUND IN '//TRIM(FNAME(IHIS))
  IF(ITYPE.EQ.2)WRITE(IU(IOUT),'(A)') 'WEIR '//TRIM(CLCID)//' NOT FOUND IN '//TRIM(FNAME(SHIS))
  NTIMES=0
  RETURN
 ENDIF

 !## fill time indices with 1-NTimes
 IF(ITYPE.EQ.1)NTIMES=DIOPLTGETNTIMES(HISDATASET)
 IF(ITYPE.EQ.2)NTIMES=DIOPLTGETNTIMES(SHISDATASET)

 IF(NTIMES.GT.DIMTIMES)THEN
  IF(ALLOCATED(SELECTEDVALUES))DEALLOCATE(SELECTEDVALUES)
  IF(ALLOCATED(JULIANTIMES))DEALLOCATE(JULIANTIMES)
  IF(ALLOCATED(TIMINDICES))DEALLOCATE(TIMINDICES)
  ALLOCATE(SELECTEDVALUES(2,1,NTIMES))
  ALLOCATE(TIMINDICES(NTIMES))
  ALLOCATE(JULIANTIMES(NTIMES))
  DIMTIMES=NTIMES
 ENDIF

 DO I=1,NTIMES
  TIMINDICES(I)=I
 ENDDO

 IF(ITYPE.EQ.1)JULIANTIMES=DIOPLTGETTIMES(HISDATASET)
 IF(ITYPE.EQ.2)JULIANTIMES=DIOPLTGETTIMES(SHISDATASET)

 !## Allocate result array (#parameters * 2 locations * #timestep), and get selection
 IF(ITYPE.EQ.1)THEN
  IF(DIOPLTGETSELECTION(HISDATASET,2,WATERLEVELINDICES,  &
                        1,SELECTEDLOCATIONINDICES,      &
                        NTIMES,TIMINDICES,SELECTEDVALUES))THEN
   IF(WATERLEVELINDICES(1).EQ.WATERLEVELINDICES(2))THEN
    !## create bottom-level from waterlevel
    SELECTEDVALUES(2,1,:)=SELECTEDVALUES(1,1,:) 
   ELSE
    !## create bottom-level from waterlevel and waterdepth
    SELECTEDVALUES(2,1,:)=SELECTEDVALUES(1,1,:)-SELECTEDVALUES(2,1,:)
   ENDIF
  ELSE
  ENDIF
 ELSEIF(ITYPE.EQ.2)THEN
  IF(DIOPLTGETSELECTION(SHISDATASET,2,WATERLEVELINDICES,  &
                        1,SELECTEDLOCATIONINDICES,      &
                        NTIMES,TIMINDICES,SELECTEDVALUES))THEN
  ELSE
  ENDIF

 ENDIF

 CALL ISDPROCESSDATA()

 END SUBROUTINE HIS1GETDATA

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

 !## compute mean waterlevel/bottomlevel
 J =INT(JULIANTIMES(1))
 K =1
 II=1
 DO I=2,NTIMES
  IF(INT(JULIANTIMES(I)).EQ.J)THEN
   K=K+1; SELECTEDVALUES(:,1,II)=SELECTEDVALUES(:,1,II)+SELECTEDVALUES(:,1,I)
  ELSE
   SELECTEDVALUES(:,1,II)=SELECTEDVALUES(:,1,II)/REAL(K,8)
   JULIANTIMES(II)       =J+1
   J                     =INT(JULIANTIMES(I))
   II                    =II+1
   SELECTEDVALUES(:,1,II)=SELECTEDVALUES(:,1,I)
   K                     =1
  ENDIF
 END DO

 SELECTEDVALUES(:,1,II)=SELECTEDVALUES(:,1,II)/REAL(K,8)
 JULIANTIMES(II)=J+1

 NTIMES=II

 END SUBROUTINE ISDPROCESSDATA

 !##=====================================================================
 SUBROUTINE SOBEK_CREATEISC1(ICRS,IB,NCRS,IREFSC)
 !##=====================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ICRS,IREFSC
 INTEGER,INTENT(OUT) :: NCRS
 INTEGER :: I,J,K,IC,II
 CHARACTER(LEN=MAXLEN) :: CROSNAME
 REAL(KIND=DP_KIND) :: MRC,MIND

 !ICRS  =ICRS+1
 !IREFSC=IREFSC+1
 !WRITE(IU(ISC1),REC=ICRS+1) 3,IREFSC,REAL(TP2(IB)%AL)/2.0,'Profile'
 !WRITE(IU(ISC2),REC=IREFSC+1) -5.0,5.0,25.0
 !IREFSC=IREFSC+1
 !WRITE(IU(ISC2),REC=IREFSC+1)  0.0D0,0.0D0,25.0
 !IREFSC=IREFSC+1
 !WRITE(IU(ISC2),REC=IREFSC+1)  5.0,5.0,25.0
 !NCRS=1
 !RETURN

 !## initialize kM
 MRC=25.0D0
 
 DO I=1,NBDFR
  IF(TRIM(TP2(IB)%ID).EQ.TRIM(BDFR(I)%CI))THEN
   MRC=BDFR(I)%MRC
   EXIT
  ENDIF
 END DO

 NCRS=0
 DO I=1,NCROS
  IF(TRIM(TP2(IB)%ID).EQ.TRIM(CR(I)%CI))THEN
   DO J=1,NDAT
    IF(TRIM(CR(I)%ID).EQ.TRIM(PDAT(J)%ID))THEN
     DO K=1,NDEF
      IF(TRIM(PDAT(J)%DI).EQ.TRIM(PDEF(K)%ID))THEN

       NCRS    =NCRS+1
       CROSNAME=PDEF(K)%ID

       ICRS    =ICRS+1
       WRITE(IU(ISC1),REC=ICRS+1) PDEF(K)%NNXY,IREFSC+1,CR(I)%LC,CROSNAME

       !## minimal value = zero!
       MIND         =MINVAL(PDEF(K)%YPROF(1:PDEF(K)%NNXY))
       PDEF(K)%YPROF=PDEF(K)%YPROF-MIND
       MIND         =0.0D0
       II           =0
       DO IC=1,PDEF(K)%NNXY
        IF(PDEF(K)%YPROF(IC).EQ.0.0D0)THEN
         MIND=MIND+PDEF(K)%XPROF(IC)
         II  =II+1
        ENDIF
       ENDDO
       MIND=MIND/REAL(II,8)
       PDEF(K)%XPROF=PDEF(K)%XPROF-MIND

       DO IC=1,PDEF(K)%NNXY
        IREFSC=IREFSC+1
        WRITE(IU(ISC2),REC=IREFSC+1) PDEF(K)%XPROF(IC),PDEF(K)%YPROF(IC),MRC
       END DO

      ENDIF
     ENDDO
    ENDIF
   ENDDO
  ENDIF
 ENDDO

 IF(NCRS.EQ.0)WRITE(IU(IOUT),'(A)') '>>> No profile found for '//TRIM(TP2(IB)%ID)//' <<<'

 END SUBROUTINE SOBEK_CREATEISC1

 !##=====================================================================
 SUBROUTINE SOBEK_IST1CREATE(ISTW,IB,NSTW,IREFST)
 !##=====================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: ISTW,IREFST
 INTEGER,INTENT(OUT) :: NSTW
 CHARACTER(LEN=MAXLEN) :: STUWNAME
 INTEGER :: I,J,N,ID

 !ISTW  =ISTW+1
 !IREFST=IREFST+1
 !WRITE(IU(IST1),REC=ISTW+1) 3,IREFST,REAL(TP2(IB)%AL)/2.0,'Stuw'
 !WRITE(IU(IST2),REC=IREFST+1) 20040101,20.0D0,15.0
 !IREFST=IREFST+1
 !WRITE(IU(IST2),REC=IREFST+1) 20040110,19.0,13.0
 !IREFST=IREFST+1
 !WRITE(IU(IST2),REC=IREFST+1) 20040131,17.0,11.0D0
 !NSTW=1

 NSTW=0
 DO I=1,NSTUW
  IF(TRIM(TP2(IB)%ID).EQ.TRIM(ST(I)%CI))THEN

 !#from his-file
   CALL HIS1GETDATA(ST(I)%ID,2)

   IF(NTIMES.GT.0)THEN

    NSTW    =NSTW+1
    ISTW    =ISTW+1
    N       =NTIMES
    STUWNAME=ST(I)%ID

    WRITE(IU(IST1),REC=ISTW+1) N,IREFST+1,ST(I)%LC,STUWNAME

    DO J=1,N
     IREFST=IREFST+1
     ID    =UTL_JDATETOIDATE(INT(JULIANTIMES(J)))
     WRITE(IU(IST2),REC=IREFST+1) ID,'00:00:00',SELECTEDVALUES(1,1,J),SELECTEDVALUES(2,1,J)
    END DO
   ENDIF

  ENDIF
 ENDDO

 END SUBROUTINE SOBEK_IST1CREATE

 !##=====================================================================
 SUBROUTINE SOBEK_ISQ1CREATE(IQHR,IB,NQHR,IREFSQ)
 !##=====================================================================
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: IB
 INTEGER,INTENT(INOUT) :: IQHR,IREFSQ
 INTEGER,INTENT(OUT) :: NQHR
 CHARACTER(LEN=MAXLEN) :: CNAME
 INTEGER :: I,N

 !IQHR  =IQHR+1
 !IREFSQ=IREFSQ+1                                          !123456789012
 !CNAME ='Qh '//TRIM(ITOS(IB))
 !WRITE(IU(ISQ1),REC=IQHR+1) 3,IREFSQ,REAL(TP2(IB)%AL)/2.0,CNAME
 !WRITE(IU(ISQ2),REC=IREFSQ+1) 20.0D0,15.0,10,6.0
 !IREFSQ=IREFSQ+1
 !WRITE(IU(ISQ2),REC=IREFSQ+1) 19.0,13.0,8.5,4.5
 !IREFSQ=IREFSQ+1
 !WRITE(IU(ISQ2),REC=IREFSQ+1) 17.0,11.0D0,6.0,2.3
 !NQHR=1

 NQHR=0
 DO I=1,NQHREL
  IF(TRIM(TP2(IB)%ID).EQ.TRIM(QH(I)%CI))THEN

   NQHR =NQHR+1
   IQHR =IQHR+1
   N    =2
   CNAME=QH(I)%ID

   IREFSQ=IREFSQ+1
   WRITE(IU(ISQ1),REC=IQHR+1) N,IREFSQ,QH(I)%LC,CNAME

   WRITE(IU(ISQ2),REC=IREFSQ+1) 20.0D0,11.0D0,10.0D0,2.3D0
   IREFSQ=IREFSQ+1
   WRITE(IU(ISQ2),REC=IREFSQ+1) 50.0D0,15.0,50.0D0,6.0D0

  ENDIF
 ENDDO

 END SUBROUTINE SOBEK_ISQ1CREATE

 !###======================================================================
 SUBROUTINE MAIN1ALNETWORKTP()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS,I

 NNODE=0; NNDLK=0; NBRCH=0
 READ(IU(ITP),*)
 DO
  READ(IU(ITP),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'NODE')THEN
   NNODE=NNODE+1
   I =INDEX(LINE,' id ')
   IF(I.LE.0)CALL MAINRDERROR('id',LINE)
  ELSEIF(INDEX(LINE,'NDLK').GT.0)THEN
   NNDLK=NNDLK+1
   I =INDEX(LINE,' id ')
   IF(I.LE.0)CALL MAINRDERROR('id',LINE)
  ELSEIF(INDEX(LINE,'BRCH').GT.0)THEN
   NBRCH=NBRCH+1
  ENDIF
 ENDDO

 REWIND(IU(ITP))

 WRITE(IU(IOUT),'(8X,I10,A)') NNODE,' nodes from '//TRIM(FNAME(ITP))
 WRITE(IU(IOUT),'(8X,I10,A)') NBRCH,' branches from '//TRIM(FNAME(ITP))
 WRITE(IU(IOUT),'(8X,I10,A)') NNDLK,' flow linkage nodes from '//TRIM(FNAME(ITP))

 END SUBROUTINE MAIN1ALNETWORKTP

 !###======================================================================
 SUBROUTINE MAIN1ALNETWORKCP()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !read number of curve-point from NETWORK.CP
 NCPNT=0
 READ(IU(ICP),*)
 DO
  READ(IU(ICP),*,IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)  
  IF(LINE(1:4).EQ.'BRCH')THEN
   READ(IU(ICP),*,IOSTAT=IOS) LINE
   IF(IOS.NE.0)EXIT
   IF(INDEX(LINE,'TBLE').GT.0)THEN
    DO
     READ(IU(ICP),*,IOSTAT=IOS) LINE
     IF(IOS.NE.0)EXIT
     IF(INDEX(LINE,'tble').GT.0)EXIT
     NCPNT=NCPNT+1
    END DO
   ENDIF
  ENDIF
 ENDDO

 REWIND(IU(ICP))

 WRITE(IU(IOUT),'(8X,I10,A)') NCPNT,' curve-points from '//TRIM(FNAME(ICP))

 END SUBROUTINE MAIN1ALNETWORKCP

 !###======================================================================
 SUBROUTINE MAIN1ALNETWORKCR()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !read number of crosssection from NETWORK.CR
 NCROS=0
 READ(IU(ICR),*)
 DO
  READ(IU(ICR),*,IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'CRSN')NCROS=NCROS+1
 ENDDO

 REWIND(IU(ICR))

 WRITE(IU(IOUT),'(8X,I10,A)') NCROS,' cross-sections from '//TRIM(FNAME(ICR))

 END SUBROUTINE MAIN1ALNETWORKCR

 !###======================================================================
 SUBROUTINE MAIN1ALNETWORKGR()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !read number of crosssection from NETWORK.CR
 NCALCP=0
 READ(IU(IGR),*)
 DO
  READ(IU(IGR),*,IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'GRID')THEN
   READ(IU(IGR),*,IOSTAT=IOS) LINE
   IF(IOS.NE.0)EXIT
   IF(INDEX(LINE,'GridPoint Table').GT.0)THEN
    READ(IU(IGR),*,IOSTAT=IOS) LINE
    IF(INDEX(LINE,'TBLE').GT.0)THEN
     DO
      READ(IU(IGR),*,IOSTAT=IOS) LINE
      IF(IOS.NE.0)EXIT
      IF(INDEX(LINE,'tble').GT.0)EXIT
      NCALCP=NCALCP+1
     END DO
    ENDIF
   ENDIF
  ENDIF
 ENDDO

 REWIND(IU(IGR))

 WRITE(IU(IOUT),'(8X,I10,A)') NCALCP,' calculation points from '//TRIM(FNAME(IGR))

 END SUBROUTINE MAIN1ALNETWORKGR

 !###======================================================================
 SUBROUTINE MAIN1ALNETWORKST()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !## read number of crosssection from NETWORK.CR
 NSTUW=0
 READ(IU(IST),*)
 DO
  READ(IU(IST),*,IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'STRU')NSTUW=NSTUW+1
 ENDDO

 REWIND(IU(IST))

 WRITE(IU(IOUT),'(8X,I10,A)') NSTUW,' weirs points from '//TRIM(FNAME(IST))

 END SUBROUTINE MAIN1ALNETWORKST

 !###======================================================================
 SUBROUTINE MAIN1ALPROFILEDAT()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !get profile id's from PROFILE.DAT
 DO
  READ(IU(IAT),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'CRSN')NDAT=NDAT+1
 ENDDO

 REWIND(IU(IAT))

 WRITE(IU(IOUT),'(8X,I10,A)') NDAT,' profile dat from '//TRIM(FNAME(IAT))

 END SUBROUTINE MAIN1ALPROFILEDAT

 !###======================================================================
 SUBROUTINE MAIN1ALPROFILEDEF()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !## get profile id's from PROFILE.DEF
 NDEF=0
 DO
  READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'CRDS')NDEF=NDEF+1
 ENDDO

 REWIND(IU(IEF))

 WRITE(IU(IOUT),'(8X,I10,A)') NDEF,' profile definitions from '//TRIM(FNAME(IEF))

 END SUBROUTINE MAIN1ALPROFILEDEF

 !###======================================================================
 SUBROUTINE MAIN1ALFRICTIONDAT()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS

 !## get friction id's from FRICTION.DAT
 NBDFR=0
 DO
  READ(IU(IFR),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'BDFR')NBDFR=NBDFR+1
 ENDDO

 REWIND(IU(IFR))

 WRITE(IU(IOUT),'(8X,I10,A)') NBDFR,' bedfriction definitions from '//TRIM(FNAME(IFR))

 END SUBROUTINE MAIN1ALFRICTIONDAT

 !###======================================================================
 SUBROUTINE MAIN1RDNETWORKTP()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: NN,NB,NL,I,IOS

 NN=0
 NB=0
 NL=0
 READ(IU(ITP),*)
 DO
  READ(IU(ITP),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'NODE')THEN

   NN=NN+1
   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) TP1(NN)%ID
   I=INDEX(LINE,' px '); IF(I.LE.0)CALL MAINRDERROR('px',LINE); READ(LINE(I+4:),*) TP1(NN)%PX
   I=INDEX(LINE,' py '); IF(I.LE.0)CALL MAINRDERROR('py',LINE); READ(LINE(I+4:),*) TP1(NN)%PY

  ELSEIF(LINE(1:4).EQ.'NDLK')THEN

   NL=NL+1
   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) TP3(NL)%ID
   I=INDEX(LINE,' px '); IF(I.LE.0)CALL MAINRDERROR('px',LINE); READ(LINE(I+4:),*) TP3(NL)%PX
   I=INDEX(LINE,' py '); IF(I.LE.0)CALL MAINRDERROR('py',LINE); READ(LINE(I+4:),*) TP3(NL)%PY
   I=INDEX(LINE,' ci '); IF(I.LE.0)CALL MAINRDERROR('ci',LINE); READ(LINE(I+4:),*) TP3(NL)%CI
   I=INDEX(LINE,' lc '); IF(I.LE.0)CALL MAINRDERROR('lc',LINE); READ(LINE(I+4:),*) TP3(NL)%LC

  ELSEIF(LINE(1:4).EQ.'BRCH')THEN

   NB=NB+1

   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) TP2(NB)%ID
   I=INDEX(LINE,' bn '); IF(I.LE.0)CALL MAINRDERROR('bn',LINE); READ(LINE(I+4:),*) TP2(NB)%CBN
   I=INDEX(LINE,' en '); IF(I.LE.0)CALL MAINRDERROR('en',LINE); READ(LINE(I+4:),*) TP2(NB)%CEN
   I=INDEX(LINE,' al '); IF(I.LE.0)CALL MAINRDERROR('al',LINE); READ(LINE(I+4:),*) TP2(NB)%AL

  ENDIF
 ENDDO

 WRITE(IU(IOUT),'(8X,I10,A)') NN,' nodes from '//TRIM(FNAME(ITP))
 WRITE(IU(IOUT),'(8X,I10,A)') NB,' branches from '//TRIM(FNAME(ITP))
 WRITE(IU(IOUT),'(8X,I10,A)') NL,' flow connection nodes from '//TRIM(FNAME(ITP))

 END SUBROUTINE MAIN1RDNETWORKTP

 !###======================================================================
 SUBROUTINE MAIN1RDNETWORKCP()
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=50) :: CPID
 INTEGER :: NC,I,J,IOS
 REAL(KIND=DP_KIND) :: D

 !## read number of curve-point from NETWORK.CP
 NC=0
 READ(IU(ICP),*)
 DO
  READ(IU(ICP),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  IF(INDEX(LINE,'BRCH').GT.0)THEN

   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) CPID

   READ(IU(ICP),'(A256)',IOSTAT=IOS) LINE 
   IF(IOS.NE.0)EXIT
   LINE=ADJUSTL(LINE)
   IF(LINE(1:4).EQ.'TBLE')THEN
    J=0
    DO
     READ(IU(ICP),'(A256)',IOSTAT=IOS) LINE
     IF(IOS.NE.0)EXIT
     IF(INDEX(LINE,'tble').GT.0)EXIT
     NC       =NC+1
     CP(NC)%ID=TRIM(CPID)
     READ(LINE,*) CP(NC)%SLEN,CP(NC)%ANGLE

     !## correct slen for end of branch
     IF(J.GT.0)THEN
      D          =CP(NC)%SLEN-CP(NC-1)%SLEN
      CP(NC)%SLEN=(D*2.0)+CP(NC-1)%SLEN
     ELSE
      CP(NC)%SLEN=CP(NC)%SLEN*2.0
     ENDIF

     J=J+1

    END DO

   ENDIF
  ENDIF
 ENDDO

 WRITE(IU(IOUT),'(8X,I10,A)') NC,' curve points from '//TRIM(FNAME(ICP))

 END SUBROUTINE MAIN1RDNETWORKCP

 !###======================================================================
 SUBROUTINE MAIN1RDNETWORKCR()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS,NC,I

 NC=0
 READ(IU(ICR),*)

 !## get profile location from NETWORK.CR
 DO
  READ(IU(ICR),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'CRSN')THEN

   NC=NC+1

   I=INDEX(LINE,' lc '); IF(I.LE.0)CALL MAINRDERROR('lc',LINE); READ(LINE(I+4:),*) CR(NC)%LC
   I=INDEX(LINE,' ci '); IF(I.LE.0)CALL MAINRDERROR('ci',LINE); READ(LINE(I+4:),*) CR(NC)%CI
   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) CR(NC)%ID

  ENDIF
 ENDDO

 WRITE(IU(IOUT),'(8X,I10,A)') NC,' cross-sections from '//TRIM(FNAME(ICR))

 END SUBROUTINE MAIN1RDNETWORKCR

 !###======================================================================
 SUBROUTINE MAIN1RDNETWORKGR()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS,NC,NQ,I,J,IB
 CHARACTER(LEN=50) :: CI,DUMMY1,DUMMY2,DUMMY3,DUMMY4

 !## read number of crosssection from NETWORK.GR
 NC=0
 NQ=0
 READ(IU(IGR),*)
 DO
  READ(IU(IGR),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'GRID')THEN

   I=INDEX(LINE,' ci ')
   IF(I.LE.0)CALL MAINRDERROR('ci',LINE)

   READ(LINE(I+4:),*) CI

   READ(IU(IGR),'(A256)',IOSTAT=IOS) LINE
   IF(IOS.NE.0)EXIT

   IF(INDEX(LINE,'GridPoint Table').GT.0)THEN

    READ(IU(IGR),*,IOSTAT=IOS) LINE
    LINE=ADJUSTL(LINE)
    IF(LINE(1:4).EQ.'TBLE')THEN
     J=0
     DO
      READ(IU(IGR),'(A256)',IOSTAT=IOS) LINE
      IF(IOS.NE.0)EXIT
      IF(INDEX(LINE,'tble').GT.0)EXIT

      J =J+1
      NC=NC+1

      GR(NC)%CI=TRIM(CI)
      READ(LINE,*) GR(NC)%LC,DUMMY1,DUMMY2,GR(NC)%ID,DUMMY3

      !## include qh-points
      IF(J.GT.1)THEN
       NQ       = NQ+1
       QH(NQ)%CI= TRIM(CI)
       QH(NQ)%LC= GR(NC-1)%LC+((GR(NC)%LC-GR(NC-1)%LC)/2.0D0)
       QH(NQ)%ID= DUMMY4
      ENDIF
      DUMMY4=DUMMY3

     END DO
    ENDIF

    DO IB=1,NBRCH
     !## adjust maximum length to actual length
     IF(TRIM(TP2(IB)%ID).EQ.TRIM(GR(NC)%CI))THEN
      GR(NC)%LC=TP2(IB)%AL
      QH(NQ)%LC=GR(NC-1)%LC+((GR(NC)%LC-GR(NC-1)%LC)/2.0d0)
     ENDIF
    ENDDO

   ENDIF
  ENDIF
 ENDDO

 NQHREL=NQ
 WRITE(IU(IOUT),'(8X,I10,A)') NC,' calculation points from '//TRIM(FNAME(IGR))
 WRITE(IU(IOUT),'(8X,I10,A)') NQ,' qh points from '//TRIM(FNAME(IGR))

 END SUBROUTINE MAIN1RDNETWORKGR

 !###======================================================================
 SUBROUTINE MAIN1RDNETWORKST()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS,NS,I

 !## read number of crosssection from NETWORK.CR
 NS=0
 READ(IU(IST),*)
 DO
  READ(IU(IST),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'STRU')THEN
   NS=NS+1

   I=INDEX(LINE,' lc '); IF(I.LE.0)CALL MAINRDERROR('lc',LINE); READ(LINE(I+4:),*) ST(NS)%LC
   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) ST(NS)%ID
   I=INDEX(LINE,' ci '); IF(I.LE.0)CALL MAINRDERROR('ci',LINE); READ(LINE(I+4:),*) ST(NS)%CI

  ENDIF
 ENDDO

 WRITE(IU(IOUT),'(8X,I10,A)') NS,' weirs points from '//TRIM(FNAME(IST))

 END SUBROUTINE MAIN1RDNETWORKST

 !###======================================================================
 SUBROUTINE MAIN1RDPROFILEDAT()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS,ND,I,J,K

 !## get profile id's from PROFILE.DAT
 ND=0
 DO
  READ(IU(IAT),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'CRSN')THEN

   I=INDEX(LINE,' id ')
   J=INDEX(LINE,' di ')
   K=INDEX(LINE,' rl ')

   IF(I.GT.0.AND.J.GT.0)THEN
    ND=ND+1
    !## get cross-section id
    READ(LINE(I+4:),*) PDAT(ND)%ID
    !## get cross-section definition-id
    READ(LINE(J+4:),*) PDAT(ND)%DI
    !## get cross-section definition-id
    IF(K.GT.0)READ(LINE(K+4:),*) PDAT(ND)%RL
   ENDIF

  ENDIF
 ENDDO

 WRITE(IU(IOUT),'(8X,I10,A)') ND,' cross-sections from '//TRIM(FNAME(IAT))

 END SUBROUTINE MAIN1RDPROFILEDAT

 !###======================================================================
 SUBROUTINE MAIN1RDPROFILEDEF()
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND) :: X,BL,BW,BS,OR,DOR,RD
 INTEGER :: IOS,ND,I,J,N
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:) :: XPROF,YPROF,XPROF_DUM,YPROF_DUM
 CHARACTER(LEN=256) :: LINE
 
 !## get profile definitions from PROFILES.DEF
 ND=0
 DO
  READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT; LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'CRDS')THEN

   ND=ND+1

   I=INDEX(LINE,' id '); IF(I.LE.0)CALL MAINRDERROR('id',LINE); READ(LINE(I+4:),*) PDEF(ND)%ID
   I=INDEX(LINE,' ty '); IF(I.LE.0)CALL MAINRDERROR('ty',LINE); READ(LINE(I+4:),*) PDEF(ND)%PTYPE

   SELECT CASE (PDEF(ND)%PTYPE)

    !## trapeziodal
    CASE (1)
     I=INDEX(LINE,' bl ')
     IF(I.LE.0)THEN
      BL=0.0D0
     ELSE
      READ(LINE(I+4:),*) BL
     ENDIF
     I=INDEX(LINE,' bw ')
     IF(I.LE.0)CALL MAINRDERROR('bw',LINE)
     READ(LINE(I+4:),*) BW
     I=INDEX(LINE,' bs ')
     IF(I.LE.0)CALL MAINRDERROR('bs',LINE)
     READ(LINE(I+4:),*) BS

     PDEF(ND)%NNXY    =4

     ALLOCATE(PDEF(ND)%XPROF(PDEF(ND)%NNXY),PDEF(ND)%YPROF(PDEF(ND)%NNXY))

     PDEF(ND)%XPROF(2)=-BW/2.0D0   
     PDEF(ND)%XPROF(1)= PDEF(ND)%XPROF(2)-5.0D0   
     PDEF(ND)%XPROF(3)= BW/2.0
     PDEF(ND)%XPROF(4)= PDEF(ND)%XPROF(3)+5.0D0   
     PDEF(ND)%YPROF(1)= BL+ABS(PDEF(ND)%XPROF(1))*BS
     PDEF(ND)%YPROF(2)= BL
     PDEF(ND)%YPROF(3)= BL
     PDEF(ND)%YPROF(4)= BL+PDEF(ND)%XPROF(4)*BS

    !## tabulated(0) --- heigth and total width for that heigth
    CASE (0)

     DO
      READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
      IF(IOS.NE.0)EXIT; LINE=ADJUSTL(LINE); IF(LINE(1:4).EQ.'TBLE')EXIT
     ENDDO
     IF(IOS.EQ.0)THEN

      PDEF(ND)%NNXY=0; ALLOCATE(XPROF(MAXPROF),YPROF(MAXPROF))

      DO
       READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
       IF(IOS.NE.0)EXIT; IF(INDEX(LINE,'tble').GT.0)EXIT
       PDEF(ND)%NNXY=PDEF(ND)%NNXY+1
       N=PDEF(ND)%NNXY
       IF(N.GT.SIZE(XPROF))THEN
        ALLOCATE(XPROF_DUM(N*2),YPROF_DUM(N*2))
        XPROF_DUM(1:N-1)=XPROF(1:N-1)
        YPROF_DUM(1:N-1)=YPROF(1:N-1)
        DEALLOCATE(XPROF,YPROF)
        XPROF=>XPROF_DUM; YPROF=>YPROF_DUM     
       ENDIF
       !## read heigth and flowing width
       READ(LINE,*) YPROF(PDEF(ND)%NNXY),X,XPROF(PDEF(ND)%NNXY)
      ENDDO

      N=PDEF(ND)%NNXY
      ALLOCATE(PDEF(ND)%XPROF(N*2),PDEF(ND)%YPROF(N*2))
      PDEF(ND)%XPROF(1:N)=XPROF(1:N)
      PDEF(ND)%YPROF(1:N)=YPROF(1:N)

      !## divide by 2
      PDEF(ND)%XPROF=PDEF(ND)%XPROF/2.0

      J=PDEF(ND)%NNXY
      DO I=1,PDEF(ND)%NNXY
       J                = J+1
       PDEF(ND)%XPROF(J)=-PDEF(ND)%XPROF(I)
       PDEF(ND)%YPROF(J)= PDEF(ND)%YPROF(I)
      END DO
      PDEF(ND)%NNXY=J
      CALL QKSORT(PDEF(ND)%NNXY,PDEF(ND)%XPROF,V2=PDEF(ND)%YPROF)
     
     ENDIF
     
    !## closed circle
    CASE (4)
    
     !## bedlevel
     I=INDEX(LINE,' bl ')
     IF(I.LE.0)THEN
      BL=0.0D0
     ELSE
      READ(LINE(I+4:),*) BL
     ENDIF

     I=INDEX(LINE,' rd ')
     IF(I.LE.0)CALL MAINRDERROR('rd',LINE)
     READ(LINE(I+4:),*) RD

     N=11; PDEF(ND)%NNXY=N

     ALLOCATE(PDEF(ND)%XPROF(PDEF(ND)%NNXY),PDEF(ND)%YPROF(PDEF(ND)%NNXY))
     OR =0.5D0*(2.0D0*PI)
     DOR=0.05D0*(2.0D0*PI)
     OR=OR+DOR
     DO I=1,N
      OR=OR-DOR
      PDEF(ND)%XPROF(I)=0.0D0+RD*COS(OR)
      PDEF(ND)%YPROF(I)= RD-RD*SIN(OR)
     ENDDO
     
    !## xz(10,11)
    CASE (10,11)

     I=INDEX(LINE,' ty ')
     !## skip tabulated storage widths
     IF(I.LE.0)THEN
      DO; READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
       IF(IOS.NE.0)EXIT; LINE=ADJUSTL(LINE); IF(LINE(1:4).EQ.'tble')EXIT
      ENDDO
     ENDIF

     DO
      READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
      IF(IOS.NE.0)EXIT; LINE=ADJUSTL(LINE); IF(LINE(1:4).EQ.'TBLE')EXIT
     ENDDO
     IF(IOS.EQ.0)THEN

      PDEF(ND)%NNXY=0; ALLOCATE(XPROF(MAXPROF),YPROF(MAXPROF))

      DO
       READ(IU(IEF),'(A256)',IOSTAT=IOS) LINE
       IF(IOS.NE.0)EXIT; IF(INDEX(LINE,'tble').GT.0)EXIT
       PDEF(ND)%NNXY=PDEF(ND)%NNXY+1
       N=PDEF(ND)%NNXY
       IF(N.GT.SIZE(XPROF))THEN
        ALLOCATE(XPROF_DUM(N*2),YPROF_DUM(N*2))
        XPROF_DUM(1:N-1)=XPROF(1:N-1)
        YPROF_DUM(1:N-1)=YPROF(1:N-1)
        DEALLOCATE(XPROF,YPROF)
        XPROF=>XPROF_DUM; YPROF=>YPROF_DUM     
       ENDIF
       !## reading horizontal and vertical distances
       READ(LINE,*) XPROF(PDEF(ND)%NNXY),YPROF(PDEF(ND)%NNXY)
      ENDDO

      N=PDEF(ND)%NNXY; ALLOCATE(PDEF(ND)%XPROF(N),PDEF(ND)%YPROF(N))
      PDEF(ND)%XPROF(1:N)=XPROF(1:N)
      PDEF(ND)%YPROF(1:N)=YPROF(1:N)

     ENDIF

    CASE DEFAULT
     LINE='Do not recognize profile type '//TRIM(ITOS(PDEF(ND)%PTYPE))
     WRITE(IU(IOUT),'(A)') TRIM(LINE)
     PDEF(ND)%NNXY=0

   END SELECT

  ENDIF
 ENDDO

 DEALLOCATE(XPROF,YPROF)

 WRITE(IU(IOUT),'(8X,I10,A)') ND,' cross-section definitions from '//TRIM(FNAME(IEF))
 
 END SUBROUTINE MAIN1RDPROFILEDEF

 !###======================================================================
 SUBROUTINE MAIN1RDFRICTIONDAT()
 !###======================================================================
 IMPLICIT NONE
 INTEGER :: IOS,NB,I,J

 !## get friction id's from FRICTION.DAT
 NB=0
 DO
  READ(IU(IFR),'(A256)',IOSTAT=IOS) LINE
  IF(IOS.NE.0)EXIT
  LINE=ADJUSTL(LINE)
  IF(LINE(1:4).EQ.'BDFR')THEN
   I=INDEX(LINE,' id ')
   J=INDEX(LINE,' ci ')

   IF(I.GT.0.AND.J.GT.0)THEN
    NB=NB+1
    !## get cross-section id
    READ(LINE(I+4:),*) BDFR(NB)%ID
    !## get cross-section definition-id
    READ(LINE(J+4:),*) BDFR(NB)%CI
   ENDIF

 !#nog iets met kd doen!
   BDFR(NB)%MRC=25.0

  ENDIF
 ENDDO

 !BDFR id '64' ci '64' mf 3 mt cp 0 25 0 mr cp 0 25 0 s1 6 s2 6     bdfr

 WRITE(IU(IOUT),'(8X,I10,A)') NB,' bedfriction definitions from '//TRIM(FNAME(IFR))

 END SUBROUTINE MAIN1RDFRICTIONDAT

 !###======================================================================
 SUBROUTINE MAINRDERROR(CERROR,LINE)
 !###======================================================================
 IMPLICIT NONE
 CHARACTER(LEN=*),INTENT(IN) :: CERROR,LINE

 WRITE(IU(IOUT),'(/1X,A/)') 'Cannot find keyword ['//TRIM(CERROR)//'] in line:'
 WRITE(IU(IOUT),'(1X,A)') TRIM(LINE)
 IF(IBATCH.EQ.0)THEN
  CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot find keyword ['//TRIM(CERROR)//'] in line:'//CHAR(13)// &
   TRIM(LINE),'Error')
 ELSE
  WRITE(*,'(/1X,A/)') 'Cannot find keyword ['//TRIM(CERROR)//'] in line:'
  WRITE(*,'(1X,A)') TRIM(LINE)
  STOP
 ENDIF
 
 END SUBROUTINE
 
END MODULE MOD_SOBEK