!! 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 . !! !! Contact: imod.support@deltares.nl !! Stichting Deltares !! P.O. Box 177 !! 2600 MH Delft, The Netherlands. !! MODULE MOD_ISG_GRID USE WINTERACTER USE RESOURCE USE MOD_PREF_PAR USE MOD_INTERSECT USE MOD_INTERSECT_PAR USE MOD_IDF_PAR USE MOD_IDF USE MODPLOT USE MOD_ISG_PAR USE MOD_PMANAGER_PAR USE MOD_UTL USE MOD_OSD USE MOD_IDFEDIT_TRACE USE MOD_ISG_TRAPEZIUM USE MOD_ISG_ADJ USE MOD_ISG_UTL USE MOD_POLYGON_UTL USE MOD_POLYGON_PAR USE MOD_IPF_PAR USE MOD_IPF USE MOD_PMANAGER_PAR USE MOD_QKSORT CONTAINS !##===================================================================== SUBROUTINE ISG_IMPORT(ISGFILE,EXPORTFNAME) !##===================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(INOUT) :: ISGFILE,EXPORTFNAME INTEGER :: I,IOS,IS,IDUMMY INTEGER,ALLOCATABLE,DIMENSION(:) :: IU CHARACTER(LEN=256) :: FNAME ALLOCATE(IU(SIZE(EXT))) DO I=1,SIZE(EXT) FNAME=EXPORTFNAME(:INDEX(EXPORTFNAME,'_',.TRUE.)-1)//'_'//TRIM(EXT(I))//'.TXT' IU(I)=UTL_GETUNIT(); CALL OSD_OPEN(IU(I),FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='READ,DENYWRITE',ACCESS='SEQUENTIAL') IF(IU(I).LE.0)RETURN ENDDO !## read ISG READ(IU(1),'(A256)') LINE READ(LINE,*,IOSTAT=IOS) NISG,ISFR IF(IOS.NE.0)THEN ISFR=0 IF(.NOT.ALLOCATED(ISDLABELS))THEN IF(ISGDOUBLE.EQ.4)THEN ALLOCATE(ISDLABELS(SIZE(TATTRIB1))); DO I=1,SIZE(ISDLABELS); ISDLABELS(I)=TATTRIB1(I); ENDDO ELSEIF(ISGDOUBLE.EQ.8)THEN ALLOCATE(ISDLABELS(SIZE(TATTRIB3))); DO I=1,SIZE(ISDLABELS); ISDLABELS(I)=TATTRIB3(I); ENDDO ENDIF ENDIF ELSE IF(.NOT.ALLOCATED(ISDLABELS))THEN IF(ISFR.EQ.0)THEN IF(ISGDOUBLE.EQ.4)THEN ALLOCATE(ISDLABELS(SIZE(TATTRIB1))); DO I=1,SIZE(ISDLABELS); ISDLABELS(I)=TATTRIB1(I); ENDDO ELSEIF(ISGDOUBLE.EQ.8)THEN ALLOCATE(ISDLABELS(SIZE(TATTRIB3))); DO I=1,SIZE(ISDLABELS); ISDLABELS(I)=TATTRIB3(I); ENDDO ENDIF ELSE ALLOCATE(ISDLABELS(SIZE(TATTRIB2))); DO I=1,SIZE(ISDLABELS); ISDLABELS(I)=TATTRIB2(I); ENDDO ENDIF ENDIF NV=UTL_COUNT_COLUMNS(LINE,','); NV=NV-2 !## overwrite default labels READ(LINE,*) NISG,ISFR,(ISDLABELS(I),I=1,NV) ENDIF ALLOCATE(ISG(NISG)); DIMISG=NISG DO IS=1,NISG READ(IU(1),*,IOSTAT=IOS) ISG(IS)%SNAME,ISG(IS)%ISEG,ISG(IS)%NSEG,ISG(IS)%ICLC,ISG(IS)%NCLC,ISG(IS)%ICRS, & ISG(IS)%NCRS, ISG(IS)%ISTW,ISG(IS)%NSTW,ISG(IS)%IQHR,ISG(IS)%NQHR IF(IOS.NE.0)THEN; WRITE(*,'(/A/)') 'Error reading '//ISGFILE(:INDEX(ISGFILE,'.',.TRUE.)-1)//'_'//TRIM(EXT(I))//'.TXT'; STOP; ENDIF ENDDO !## determine amount of memory consumed by *.isp DIMISP=SUM(ISG(1:NISG)%NSEG); ALLOCATE(ISP(DIMISP),STAT=IOS) !## determine amount of memory consumed by *.isd part I DIMISD=SUM(ISG(1:NISG)%NCLC); ALLOCATE(ISD(DIMISD),STAT=IOS) !## determine amount of memory consumed by *.isc part I DIMISC=SUM(ISG(1:NISG)%NCRS); ALLOCATE(ISC(DIMISC),STAT=IOS) !## determine amount of memory consumed by *.ist part I DIMIST=SUM(ISG(1:NISG)%NSTW); ALLOCATE(IST(DIMIST),STAT=IOS) !## determine amount of memory consumed by *.isq part I DIMISQ=SUM(ISG(1:NISG)%NQHR); ALLOCATE(ISQ(DIMISQ),STAT=IOS) READ(IU(2),*); DO NISP=1,SIZE(ISP) READ(IU(2),*) IDUMMY,ISP(NISP)%X,ISP(NISP)%Y ENDDO; CLOSE(IU(2)); !## read isd1/isd2 READ(IU(3),*); DIMDATISD=0 DO NISD=1,SIZE(ISD) READ(IU(3),*) IDUMMY,ISD(NISD)%N,ISD(NISD)%IREF,ISD(NISD)%DIST,ISD(NISD)%CNAME DIMDATISD=DIMDATISD+ISD(NISD)%N ENDDO CLOSE(IU(3)); ALLOCATE(DATISD(DIMDATISD)) READ(IU(4),*) IF(ISFR.EQ.0)THEN DO I=1,SIZE(DATISD) READ(IU(4),*) IDUMMY,DATISD(I)%IDATE,DATISD(I)%WLVL,DATISD(I)%BTML,DATISD(I)%RESIS,DATISD(I)%INFF ENDDO ELSE DO I=1,SIZE(DATISD) READ(IU(4),*) IDUMMY,DATISD(I)%IDATE,DATISD(I)%CTIME,DATISD(I)%WLVL, & DATISD(I)%BTML, DATISD(I)%WIDTH,DATISD(I)%THCK ,DATISD(I)%HCND,DATISD(I)%UPSG, DATISD(I)%DWNS, DATISD(I)%ICLC, & DATISD(I)%IPRI ,DATISD(I)%QFLW, DATISD(I)%QROF,DATISD(I)%PPTSW, DATISD(I)%ETSW ENDDO ENDIF CLOSE(IU(4)) !## read isc1/isc2 READ(IU(5),*); DIMDATISC=0 DO NISC=1,SIZE(ISC) READ(IU(5),*) IDUMMY,ISC(NISC)%N,ISC(NISC)%IREF,ISC(NISC)%DIST,ISC(NISC)%CNAME DIMDATISC=DIMDATISC+ISC(NISC)%N ENDDO CLOSE(IU(5)); ALLOCATE(DATISC(DIMDATISC)) READ(IU(6),*) DO I=1,SIZE(DATISC) READ(IU(6),*) IDUMMY,DATISC(I)%DISTANCE,DATISC(I)%BOTTOM,DATISC(I)%MRC,DATISC(I)%ZP ENDDO CLOSE(IU(6)) !## read ist1/ist2 READ(IU(7),*); DIMDATIST=0 DO NIST=1,SIZE(IST) READ(IU(7),*) IDUMMY,IST(NIST)%N,IST(NIST)%IREF,IST(NIST)%DIST,IST(NIST)%CNAME DIMDATIST=DIMDATIST+IST(NIST)%N ENDDO CLOSE(IU(7)); ALLOCATE(DATIST(DIMDATIST)) READ(IU(8),*) DO I=1,SIZE(DATIST) READ(IU(8),*) IDUMMY,DATIST(I)%IDATE,DATIST(I)%WLVL_UP,DATIST(I)%WLVL_DOWN ENDDO CLOSE(IU(8)) !## read isq1/isq2 READ(IU(9),*); DIMDATISQ=0 DO NISQ=1,SIZE(ISQ) READ(IU(9),*) IDUMMY,ISQ(NISQ)%N,ISQ(NISQ)%IREF,ISQ(NISQ)%DIST,ISQ(NISQ)%CNAME DIMDATISQ=DIMDATISQ+ISQ(NISQ)%N ENDDO CLOSE(IU(9)); ALLOCATE(DATISQ(DIMDATISQ)) READ(IU(10),*) DO I=1,SIZE(DATISQ) READ(IU(10),*) IDUMMY,DATISQ(I)%Q,DATISQ(I)%W,DATISQ(I)%D,DATISQ(I)%F ENDDO CLOSE(IU(10)) DO I=1,SIZE(IU); CLOSE(IU(I)); ENDDO CALL ISGSAVE(ISGFILE,2) END SUBROUTINE ISG_IMPORT !##===================================================================== SUBROUTINE ISG_EXPORT(ISGFILE,EXPORTFNAME,IEXPORT,IBATCH) !##===================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: ISGFILE CHARACTER(LEN=*),INTENT(IN) :: EXPORTFNAME INTEGER,INTENT(IN) :: IEXPORT,IBATCH INTEGER :: I,IU,NPNT,IPNT,NCLC,ICLC,J,ICRS,NCRS REAL(KIND=DP_KIND) :: DIST,TD CHARACTER(LEN=256) :: FNAME !## read entire ISG file IF(ISGREAD((/ISGFILE/),IBATCH))THEN; ENDIF IF(IEXPORT.EQ.0)THEN FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISG.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) CALL IOSCOPYFILE(ISGFILE,FNAME) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISP.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,2(1X,A15))') 'RECORD','X','Y' DO I=1,SIZE(ISP); WRITE(IU,'(I10,2(1X,F15.3))') I,ISP(I)%X,ISP(I)%Y; ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISD1.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,2(1X,A10),1X,A15,1X,A32)') 'RECORD','N','IREF','DIST','CNAME' DO I=1,SIZE(ISD); WRITE(IU,'(I10,2(1X,I10),1X,F15.3,1X,A)') I,ISD(I)%N,ISD(I)%IREF,ISD(I)%DIST,'"'//TRIM(ISD(I)%CNAME)//'"'; ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISD2.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') IF(ISFR.EQ.0)THEN WRITE(IU,'(A10,1X,A10,4(1X,A15))') 'RECORD','IDATE','WLVL','BTML','RESIS','INFF' DO I=1,SIZE(DATISD) WRITE(IU,'(I10,1X,I10,4(1X,F15.3))') I,DATISD(I)%IDATE,DATISD(I)%WLVL,DATISD(I)%BTML,DATISD(I)%RESIS,DATISD(I)%INFF ENDDO ELSEIF(ISFR.EQ.1)THEN WRITE(IU,'(A10,2(1X,A10),5(1X,A15),4(1X,A10),4(1X,A15))') 'RECORD','IDATE','CTIME','WLVL','BTML','WIDTH','THCK','HCND', & 'UPSG','DWNS','ICLC','IPRI','QFLW','QROF','PPTSW','ETSW' DO I=1,SIZE(DATISD) WRITE(IU,'(I10,1X,I10,1X,A10,5(1X,F15.3),4(1X,I10),4(1X,F15.3))') I,DATISD(I)%IDATE,DATISD(I)%CTIME,DATISD(I)%WLVL, & DATISD(I)%BTML, DATISD(I)%WIDTH,DATISD(I)%THCK ,DATISD(I)%HCND,DATISD(I)%UPSG, DATISD(I)%DWNS, DATISD(I)%ICLC, & DATISD(I)%IPRI ,DATISD(I)%QFLW, DATISD(I)%QROF,DATISD(I)%PPTSW, DATISD(I)%ETSW ENDDO ENDIF CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISC1.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,2(1X,A10),1X,A15,1X,A32)') 'RECORD','N','IREF','DIST','CNAME' DO I=1,SIZE(ISC); WRITE(IU,'(I10,2(1X,I10),1X,F15.3,1X,A)') I,ISC(I)%N,ISC(I)%IREF,ISC(I)%DIST,'"'//TRIM(ISC(I)%CNAME)//'"'; ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISC2.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,4(1X,A15))') 'RECORD','DISTANCE','BOTTOM','MRC','ZP' DO I=1,SIZE(DATISC) WRITE(IU,'(I10,4(1X,F15.3))') I,DATISC(I)%DISTANCE,DATISC(I)%BOTTOM,DATISC(I)%MRC,DATISC(I)%ZP ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_IST1.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,2(1X,A10),1X,A15,1X,A32)') 'RECORD','N','IREF','DIST','CNAME' DO I=1,SIZE(IST); WRITE(IU,'(I10,2(1X,I10),1X,F15.3,1X,A)') I,IST(I)%N,IST(I)%IREF,IST(I)%DIST,'"'//TRIM(IST(I)%CNAME)//'"'; ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_IST2.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,1X,A10,2(1X,A15))') 'RECORD','IDATE','WLVL_UP','WLVL_DOWN' DO I=1,SIZE(DATIST) WRITE(IU,'(I10,1X,I10,2(1X,F15.3))') I,DATIST(I)%IDATE,DATIST(I)%WLVL_UP,DATIST(I)%WLVL_DOWN ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISQ1.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,2(1X,A10),1X,A15,1X,A32)') 'RECORD','N','IREF','DIST','CNAME' DO I=1,SIZE(ISQ); WRITE(IU,'(I10,2(1X,I10),1X,F15.3,1X,A)') I,ISQ(I)%N,ISQ(I)%IREF,ISQ(I)%DIST,'"'//TRIM(ISQ(I)%CNAME)//'"'; ENDDO CLOSE(IU) FNAME=ISGFILE(INDEX(ISGFILE,'\',.TRUE.)+1:INDEX(ISGFILE,'.',.TRUE.)-1)//'_ISQ2.TXT' FNAME=TRIM(EXPORTFNAME)//'\'//TRIM(FNAME) IU=UTL_GETUNIT(); CALL OSD_OPEN(IU,FILE=FNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IU,'(A10,4(1X,A15))') 'RECORD','Q','W','D','F' DO I=1,SIZE(DATISQ) WRITE(IU,'(I10,4(1X,F15.3))') I,DATISQ(I)%Q,DATISQ(I)%W,DATISQ(I)%D,DATISQ(I)%F ENDDO CLOSE(IU) ELSE !## read associated textfile for selected location IU=UTL_GETUNIT() CALL OSD_OPEN(IU,FILE=EXPORTFNAME,STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') IF(IEXPORT.EQ.1)THEN WRITE(IU,'(A11,A51,4A11,A32)') 'ISEGMENT,','SEGNAME,','ICALC,','DISTANCE,','X,','Y,','CALC.NAME' ELSEIF(IEXPORT.EQ.2)THEN WRITE(IU,'(A11,A51,4A11,A32)') 'ISEGMENT,','SEGNAME,','ICROSS,','DISTANCE,','X,','Y,','CROSS.NAME' ENDIF !## proces each calculation point in the ISG DO ISELISG=1,NISG !## number of nodes on segment NPNT=ISG(ISELISG)%NSEG; IPNT=ISG(ISELISG)%ISEG IF(IEXPORT.EQ.1)THEN NCLC=ISG(ISELISG)%NCLC; ICLC=ISG(ISELISG)%ICLC DO J=1,NCLC DIST=ISD(ICLC+J-1)%DIST !## compute correct x/y coordinate of current computational node CALL ISGADJUSTCOMPUTEXY(IPNT,NPNT,DIST,TD) WRITE(IU,'(I10,A1,A50,A1,I10,3(A1,F10.2),A1,A32)') ISELISG,',',ADJUSTR(ISG(ISELISG)%SNAME),',',J,',', & DIST,',',ISGX,',',ISGY,',',ADJUSTR(ISD(ICLC+J-1)%CNAME) ENDDO ELSEIF(IEXPORT.EQ.2)THEN NCRS=ISG(ISELISG)%NCRS; ICRS=ISG(ISELISG)%ICRS DO J=1,NCRS DIST=ISC(ICRS+J-1)%DIST !## compute correct x/y coordinate of current computational node CALL ISGADJUSTCOMPUTEXY(IPNT,NPNT,DIST,TD) WRITE(IU,'(I10,A1,A50,A1,I10,3(A1,F10.2),A1,A32)') ISELISG,',',ADJUSTR(ISG(ISELISG)%SNAME),',',J,',', & DIST,',',ISGX,',',ISGY,',',ADJUSTR(ISC(ICRS+J-1)%CNAME) ENDDO ENDIF ENDDO ENDIF END SUBROUTINE ISG_EXPORT !##===================================================================== SUBROUTINE ISG_ADDSTAGES(ISGFILE,IPFFILE,IBATCH,STAGETYPE,ICLEAN) !##===================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: ISGFILE,IPFFILE INTEGER,INTENT(IN) :: IBATCH,STAGETYPE,ICLEAN INTEGER :: IPNT,NPNT,ICLC,NCLC,I,J,K,IPOS,ISEG,N,ID,NTXT,MTXT,IU,IREF INTEGER,ALLOCATABLE,DIMENSION(:) :: ISORT REAL(KIND=DP_KIND) :: MD,DIST,TD,DX,INFF,BTML,RESIS,STAGE CHARACTER(LEN=256) :: FNAME CHARACTER(LEN=14) :: CDATE INTEGER,ALLOCATABLE,DIMENSION(:) :: FIDATE,FITYPE REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: FSTAGE,FNODATA CHARACTER(LEN=52),ALLOCATABLE,DIMENSION(:) :: FLABEL !## allocate memory for ipf-plotting, they will be read in memory and drawn from that NIPF=1; CALL IPFALLOCATE() IPF(1)%FNAME=IPFFILE IPF(1)%XCOL=1; IPF(1)%YCOL=2; IPF(1)%ZCOL=2; IPF(1)%Z2COL=2; IPF(1)%QCOL=2; IPF(1)%ITYPE=0 IF(.NOT.IPFREAD2(1,1,1))RETURN !## read entire ISG file IF(ISGREAD((/ISGFILE/),IBATCH))THEN; ENDIF !## proces each calculation point in the ISG DO ISELISG=1,NISG !## number of nodes on segment NPNT=ISG(ISELISG)%NSEG; IPNT=ISG(ISELISG)%ISEG NCLC=ISG(ISELISG)%NCLC; ICLC=ISG(ISELISG)%ICLC DO J=1,NCLC DIST=ISD(ICLC+J-1)%DIST !## compute correct x/y coordinate of current computational node CALL ISGADJUSTCOMPUTEXY(IPNT,NPNT,DIST,TD) !## get closest point from ipf ID=0 DO I=1,IPF(1)%NROW DX=SQRT((ISGX-IPF(1)%XYZ(1,I))**2.0D0+(ISGY-IPF(1)%XYZ(2,I))**2.0D0) IF(I.EQ.1)THEN; MD=DX; ID=I; ENDIF IF(DX.LT.MD)THEN MD=DX; ID=I ENDIF ENDDO WRITE(*,'(A,F10.2)') 'Location found for calculation node '//TRIM(ISD(ICLC+J-1)%CNAME)//';distance ',MD FNAME=IPFFILE(:INDEX(IPFFILE,'\',.TRUE.)-1)//'\'//TRIM(IPF(1)%INFO(IPF(1)%ACOL,ID))//'.txt' WRITE(*,'(A)') ' - Reading associated file: '//TRIM(FNAME) !## read associated textfile for selected location IU=UTL_GETUNIT() CALL OSD_OPEN(IU,FILE=FNAME,STATUS='OLD',FORM='FORMATTED',ACTION='READ,DENYWRITE',ACCESS='SEQUENTIAL') IF(IU.LE.0)STOP READ(IU,*) NTXT; READ(IU,*) MTXT !## get location of calculation node IPOS=ISG(ISELISG)%ICLC-1+J !## store total content K=NTXT+ISD(IPOS)%N ALLOCATE(FITYPE(K),FIDATE(K),FSTAGE(K),FLABEL(MTXT),FNODATA(MTXT),ISORT(K)) DO I=1,MTXT; READ(IU,*) FLABEL(I),FNODATA(I); ENDDO K=0 DO I=1,NTXT READ(IU,*) CDATE,STAGE IF(STAGE.NE.FNODATA(2))THEN K=K+1 READ(CDATE,'(I8)') FIDATE(K) FSTAGE(K)=STAGE !## flag is 1, to set priority FITYPE(K)=1 ENDIF ENDDO NTXT=K CLOSE(IU) !## get available date TO start insertion IREF=ISD(IPOS)%IREF-1 DO I=1,ISD(IPOS)%N IREF=IREF+1 K=K+1 FIDATE(K)=DATISD(IREF)%IDATE FSTAGE(K)=DATISD(IREF)%WLVL FITYPE(K)=0 ENDDO !## store latest information BTML =DATISD(IREF)%BTML RESIS=DATISD(IREF)%RESIS INFF =DATISD(IREF)%INFF !## correct stage in case relative stages are entered IF(STAGETYPE.EQ.1)THEN; DO I=1,K; IF(FITYPE(I).EQ.0)FSTAGE(I)=FSTAGE(I)-BTML; ENDDO; ENDIF !## sort date - DO I=1,K; ISORT(I)=I; ENDDO CALL QKSORT_INT(K,FIDATE,V2=ISORT) !## get number of unique dates N=0 DO I=1,K !## check flag to remove all existing data as well IF(ICLEAN.EQ.1)THEN IF(FITYPE(ISORT(I)).EQ.0)CYCLE ENDIF !## for doubles only use itype.eq.1 IF(I.LT.K)THEN IF(FIDATE(I+1).EQ.FIDATE(I))THEN IF(FITYPE(ISORT(I)).EQ.0)CYCLE ENDIF ENDIF IF(I.GT.1)THEN IF(FIDATE(I-1).EQ.FIDATE(I))THEN IF(FITYPE(ISORT(I)).EQ.0)CYCLE ENDIF ENDIF N=N+1 ENDDO !## increase/decrease memory data calculation node - data will be replaced to the back CALL ISGMEMORYDATISD(N-ISD(IPOS)%N,IPOS,ISEG) !## rewrite data IREF=ISEG-1 DO I=1,K !## check flag to remove all existing data as well IF(ICLEAN.EQ.1)THEN IF(FITYPE(ISORT(I)).EQ.0)CYCLE ENDIF !## for doubles only use itype.eq.1 IF(I.LT.K)THEN IF(FIDATE(I+1).EQ.FIDATE(I))THEN IF(FITYPE(ISORT(I)).EQ.0)CYCLE ENDIF ENDIF IF(I.GT.1)THEN IF(FIDATE(I-1).EQ.FIDATE(I))THEN IF(FITYPE(ISORT(I)).EQ.0)CYCLE ENDIF ENDIF IREF=IREF+1 DATISD(IREF)%IDATE=FIDATE(I) IF(STAGETYPE.EQ.0)DATISD(IREF)%WLVL =FSTAGE(ISORT(I)) IF(STAGETYPE.EQ.1)DATISD(IREF)%WLVL =BTML+FSTAGE(ISORT(I)) DATISD(IREF)%BTML =BTML DATISD(IREF)%RESIS=RESIS DATISD(IREF)%INFF =INFF ENDDO DEALLOCATE(FLABEL,FITYPE,FNODATA,FSTAGE,FIDATE,ISORT) ENDDO ENDDO CALL IPFDEALLOCATE() END SUBROUTINE ISG_ADDSTAGES !##===================================================================== SUBROUTINE ISG_ADDCROSSSECTION(ISGFILE,FNAME,WIDTHFNAME,MAXDIST,CROSS_PNTR, & CROSS_BATH,CROSS_ZCHK,CROSS_CVAL,IBATCH,ICLEAN) !##===================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: ISGFILE,FNAME,WIDTHFNAME,CROSS_PNTR,CROSS_BATH,CROSS_ZCHK,CROSS_CVAL REAL(KIND=DP_KIND),INTENT(IN) :: MAXDIST INTEGER,INTENT(IN) :: IBATCH,ICLEAN INTEGER :: IU,IOS,N,I,J,IISG,IPOS,ISEG,NSEG,ICRS,NCRS,IPNT,NPNT,IROW,ICOL,IP,ZCHK,CVAL,NI,IN,ICLC,IREF,IT,NT INTEGER(KIND=1) :: CF CHARACTER(LEN=7500) :: STRING REAL(KIND=DP_KIND),DIMENSION(1000) :: X,Y REAL(KIND=DP_KIND) :: XCRD,YCRD,TDIST,DIST,W,TD,XC,YC,TH CHARACTER(LEN=MAXLENISG) :: LABEL TYPE(IDFOBJ) :: IDF TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: ICROSS REAL(KIND=DP_KIND),DIMENSION(:,:),ALLOCATABLE :: PISG !## read entire ISG file IF(ISGREAD((/ISGFILE/),IBATCH))THEN; ENDIF IF(CROSS_PNTR.EQ.'')THEN IF(ALLOCATED(PISG))DEALLOCATE(PISG); ALLOCATE(PISG(NISG,3)) IF(TRIM(WIDTHFNAME).NE.'')THEN IF(.NOT.IDFREAD(IDF,WIDTHFNAME,0))STOP 'Cannot read width fname' ENDIF !## remove all cross-sections on segment IF(ICLEAN.EQ.1)THEN WRITE(*,'(/A/)') 'Removing all cross-sections ...' DO ISELISG=1,NISG ICRS=ISG(ISELISG)%ICRS; NCRS=ISG(ISELISG)%NCRS DO I=1,NCRS; CALL ISGDELISC(ISELISG,ICRS); END DO ENDDO ISC%N=0; ISC%IREF=0; ISC%DIST=0.0D0; ISC%CNAME='' ENDIF IU=UTL_GETUNIT() CALL OSD_OPEN(IU,FILE=FNAME,STATUS='OLD',ACTION='READ') NT=0; DO READ(IU,*,IOSTAT=IOS) STRING; IF(IOS.NE.0)EXIT; NT=NT+1 ENDDO REWIND(IU) WRITE(*,'(/A,I10,A/)') 'Adding ',NT,' cross-sections from '//TRIM(FNAME) IT=0; DO IT=IT+1 READ(IU,'(A7500)',IOSTAT=IOS) STRING; IF(IOS.NE.0)EXIT; IF(LEN_TRIM(STRING).EQ.0)EXIT IF(STRING(LEN(STRING):LEN(STRING)).NE.'')THEN WRITE(*,'(A)') 'Reduce number of cross-section definitions to be less than '//TRIM(ITOS(LEN(STRING)))//' characters' ENDIF READ(STRING,*,IOSTAT=IOS) XCRD,YCRD,LABEL,N,(X(I),I=1,N),(Y(I),I=1,N) !## not enough points IF(N.LE.2)THEN CYCLE ENDIF !## check cross-section DO I=2,N IF(X(I).LT.X(I-1))THEN WRITE(*,'(/A)') 'Wrong cross-section read' WRITE(*,'(A/)') 'The x-values need to be increasing from left to right' WRITE(*,*) X(I),X(I-1) WRITE(*,*) TRIM(LABEL) STOP ENDIF ENDDO !## find distance of point on segment iisg that is closest CALL ISGSTUWEN_INTERSECT(MAXDIST,XCRD,YCRD,PISG,NI) !## put cross-section on all segments within distance DO IN=1,NI TDIST=PISG(IN,2) IISG =INT(PISG(IN,1)) !## remove cross-sections if exists on current segment IF(ICLEAN.EQ.2)THEN ICRS=ISG(IISG)%ICRS; NCRS=ISG(IISG)%NCRS DO I=1,NCRS IF(ABS(TDIST-ISC(ICRS)%DIST).LE.MAXDIST)THEN CALL ISGDELISC(IISG,ICRS) ELSE ICRS=ICRS+1 ENDIF ENDDO ENDIF ISELISG=IISG !## if sfr-isg file - change calc option to 2 IF(ISFR.EQ.1)THEN ICLC=ISG(IISG)%ICLC; IREF=ISD(ICLC)%IREF; DATISD(IREF)%ICLC=3 ENDIF !## increase memory location cross-section CALL ISGMEMORYISC(1,ISELISG,IPOS) ISC(IPOS)%N =0 ISC(IPOS)%IREF =NDISC+1 ISC(IPOS)%DIST =TDIST ISC(IPOS)%CNAME=TRIM(LABEL) !## increase memory data cross-section CALL ISGMEMORYDATISC(N,IPOS,ISEG) DO I=1,N DATISC(ISEG+I-1)%DISTANCE=X(I) DATISC(ISEG+I-1)%BOTTOM =Y(I) DATISC(ISEG+I-1)%MRC =0.03D0 ENDDO ENDDO WRITE(6,'(A,F10.2,A)') '+Progress ',REAL(IT*100)/REAL(NT),'% ' ENDDO CLOSE(IU) IF(TRIM(WIDTHFNAME).NE.'')THEN WRITE(*,'(/A/)') 'Adding default cross-sections for remaining segments' !## put cross-sections on segment that did not got a cross-section DO ISELISG=1,NISG ICRS=ISG(ISELISG)%ICRS; NCRS=ISG(ISELISG)%NCRS IF(NCRS.EQ.0)THEN !## get distance of segment, put cross-section in mid ISEG =ISG(ISELISG)%ISEG; NSEG =ISG(ISELISG)%NSEG TDIST=0.0D0 DO I=2,NSEG ISEG=ISEG+1; DIST=(ISP(ISEG)%X-ISP(ISEG-1)%X)**2.0D0+(ISP(ISEG)%Y-ISP(ISEG-1)%Y)**2.0D0 IF(DIST.GT.0.0D0)DIST=SQRT(DIST); TDIST=TDIST+DIST END DO TDIST=TDIST/2.0D0 !## get x,y coordinates for current cross-section location ISEG=ISG(ISELISG)%ISEG CALL ISGGETXY(ISP(ISEG:)%X,ISP(ISEG:)%Y,NSEG,TDIST,XCRD,YCRD) W=ABS(IDFGETXYVAL(IDF,XCRD,YCRD)) N=4 !## increase memory CALL ISGMEMORYISC(1,ISELISG,IPOS) ISC(IPOS)%N =0 ISC(IPOS)%IREF =NDISC+1 ISC(IPOS)%DIST =TDIST ISC(IPOS)%CNAME='Cross-Section '//TRIM(ITOS(ISELISG)) CALL ISGMEMORYDATISC(N,IPOS,ISEG) DATISC(ISEG)%DISTANCE =-W/2.0D0 DATISC(ISEG)%BOTTOM = 5.0D0 DATISC(ISEG)%MRC =0.03D0 DATISC(ISEG+1)%DISTANCE=-W/2.0D0 DATISC(ISEG+1)%BOTTOM = 0.0D0 DATISC(ISEG+1)%MRC =0.03D0 DATISC(ISEG+2)%DISTANCE= W/2.0D0 DATISC(ISEG+2)%BOTTOM = 0.0D0 DATISC(ISEG+2)%MRC =0.03D0 DATISC(ISEG+3)%DISTANCE= W/2.0D0 DATISC(ISEG+3)%BOTTOM = 5.0D0 DATISC(ISEG+3)%MRC =0.03D0 ENDIF ENDDO ENDIF IF(ALLOCATED(PISG))DEALLOCATE(PISG) ELSE ALLOCATE(ICROSS(4)); DO I=1,SIZE(ICROSS); CALL IDFNULLIFY(ICROSS(I)); ENDDO !## read additional reference heights ZCHK=0; IF(TRIM(CROSS_ZCHK).NE.'')ZCHK=1 !## read additional resistances CVAL=0; IF(TRIM(CROSS_CVAL).NE.'')CVAL=1 IF(.NOT.IDFREAD(ICROSS(1),CROSS_PNTR,1))RETURN !## read pointer !## read zval at pointer scale CALL IDFCOPY(ICROSS(1),ICROSS(2)) IF(.NOT.IDFREADSCALE(CROSS_BATH,ICROSS(2),2,1,0.0D0,0))RETURN !## average value IF(ZCHK.EQ.1)THEN !## read zchk values CALL IDFCOPY(ICROSS(1),ICROSS(3)) IF(.NOT.IDFREADSCALE(CROSS_ZCHK,ICROSS(3),2,1,0.0D0,0))RETURN !## average value ENDIF IF(CVAL.EQ.1)THEN !## read cval values CALL IDFCOPY(ICROSS(1),ICROSS(4)) IF(.NOT.IDFREADSCALE(CROSS_CVAL,ICROSS(4),2,1,0.0D0,0))RETURN !## average value ENDIF DO ISELISG=1,NISG !## number of nodes on segment NPNT=ISG(ISELISG)%NSEG; IPNT=ISG(ISELISG)%ISEG ICRS=ISG(ISELISG)%ICRS; NCRS=ISG(ISELISG)%NCRS WRITE(6,'(2(A,I10),A)') '+Busy with segment ',ISELISG,' adding ',NCRS,' cross-sections' DO J=1,NCRS IF(J.EQ.NCRS)THEN WRITE(*,*) 'DSDS' ENDIF DIST=ISC(ICRS+J-1)%DIST !## compute correct x/y coordinate of current cross-section CALL ISGADJUSTCOMPUTEXY(IPNT,NPNT,DIST,TD) !## get irow/icol for current location of cross-section CALL IDFIROWICOL(ICROSS(1),IROW,ICOL,ISGX,ISGY) IF(IROW.NE.0.AND.ICOL.NE.0)THEN IP=ICROSS(1)%X(ICOL,IROW) IF(IP.NE.ICROSS(1)%NODATA.AND.ABS(IP).GT.0)THEN TH=ICROSS(3)%X(ICOL,IROW) !## get number of cross-section locations N=0; DO IROW=1,ICROSS(1)%NROW; DO ICOL=1,ICROSS(1)%NCOL !## location of gridcell equal to pointer value at location of cross-section IF(ABS(ICROSS(1)%X(ICOL,IROW)).EQ.ABS(IP))THEN IF(ICROSS(2)%X(ICOL,IROW).NE.ICROSS(2)%NODATA)N=N+1 ENDIF ENDDO; ENDDO !## add extra record to store dx,dy N=N+1 !## get location of cross-sections IPOS=ISG(ISELISG)%ICRS-1+J !## increase/decrease memory data cross-section N=N-ABS(ISC(IPOS)%N) ! if(n.ne.0)then ! write(*,*) j,n ! endif CALL ISGMEMORYDATISC(N,IPOS,ISEG) ISC(IPOS)%N=-1.0D0*ABS(ISC(IPOS)%N) ! !## check references ! do ii=1,ncrs ! iii=isg(iselisg)%icrs-1+ii ! iiii=isc(iii)%iref ! if(datisc(iiii)%distance.gt.0.0D0)then ! write(*,*) 'wrong' ! endif ! enddo IF(ZCHK.EQ.0)THEN DATISC(ISEG)%DISTANCE= ICROSS(1)%DX DATISC(ISEG)%BOTTOM = ICROSS(1)%DY DATISC(ISEG)%MRC = 0.0D0 !## empty, not to be used (yet) ELSE DATISC(ISEG)%DISTANCE=-ICROSS(1)%DX DATISC(ISEG)%BOTTOM =-ICROSS(1)%DY DATISC(ISEG)%MRC = TH !## threshold ENDIF N=0; DO IROW=1,ICROSS(1)%NROW; DO ICOL=1,ICROSS(1)%NCOL !## location of gridcell equal to pointer value at location of cross-section IF(ABS(ICROSS(1)%X(ICOL,IROW)).EQ.ABS(IP))THEN IF(ICROSS(2)%X(ICOL,IROW).NE.ICROSS(2)%NODATA)THEN N=N+1 CALL IDFGETLOC(ICROSS(1),IROW,ICOL,XC,YC) DATISC(ISEG+N)%DISTANCE=XC DATISC(ISEG+N)%BOTTOM =YC DATISC(ISEG+N)%MRC =ICROSS(2)%X(ICOL,IROW) IF(ZCHK.EQ.1)THEN CF=INT(1,1) IF(CVAL.EQ.1)THEN IF(INT(ICROSS(4)%X(ICOL,IROW)).LE.HUGE(CF))CF=INT(ICROSS(4)%X(ICOL,IROW)) ENDIF IF(ICROSS(1)%X(ICOL,IROW).GT.0)DATISC(ISEG+N)%ZP=INT( CF,1) IF(ICROSS(1)%X(ICOL,IROW).LT.0)DATISC(ISEG+N)%ZP=INT(-CF,1) ENDIF ENDIF ENDIF ENDDO; ENDDO ELSE WRITE(*,'(/A,I10)') 'Pointer value is equal to nodata value/or le 0, IP=',IP WRITE(*,'(2(A,I10)/)') 'For cross-section ',J,' on segment ',ISELISG ENDIF ELSE WRITE(*,'(/2(A,I10),A/)') 'Position of cross-section ',J,' on segment ',ISELISG,' is outside Pointer IDF' ENDIF ENDDO ENDDO CALL IDFDEALLOCATE(ICROSS,SIZE(ICROSS)); DEALLOCATE(ICROSS) ENDIF END SUBROUTINE ISG_ADDCROSSSECTION !##===================================================================== SUBROUTINE ISG_SIMPLIFYMAIN(ISGFILE,ZTOLERANCE,NODATA,IBATCH) !##===================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: ISGFILE REAL(KIND=DP_KIND),INTENT(IN) :: ZTOLERANCE,NODATA INTEGER,INTENT(IN) :: IBATCH INTEGER :: I,J,K,ITYPE,ISS,IREF,NCLC,IAVERAGE REAL(KIND=DP_KIND),DIMENSION(4) :: RVAL,XNR REAL(KIND=DP_KIND),DIMENSION(1,4) :: QSORT,NDATA REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: ZDIST,XDIST,GCODE INTEGER,ALLOCATABLE,DIMENSION(:) :: ILIST !## read entire ISG file IF(ISGREAD((/ISGFILE/),IBATCH))THEN; ENDIF ITYPE=1 !## itype=1: isd ISS=1 !## mean value NCLC=MAXVAL(ISG%NCLC); ALLOCATE(ZDIST(NCLC),XDIST(NCLC),GCODE(NCLC),ILIST(NCLC)) NDATA=NODATA !## arithmetic mean IAVERAGE=1 !## process each segment in ISG file to simplify calculation points DO I=1,NISG !## get mean waterlevels in segment IREF=ISG(I)%ICLC-1; NCLC=ISG(I)%NCLC DO J=1,NCLC IREF=IREF+1 CALL ISG2GRIDGETDATA(0,0,1,QSORT,XNR,4,RVAL,ISD(IREF)%N,ISD(IREF)%IREF,ISS,ITYPE,NDATA,IAVERAGE) XDIST(J)=ISD(IREF)%DIST; ZDIST(J)=RVAL(1) IF(SUM(XNR)/4.NE.XNR(1))ZDIST(J)=NODATA ENDDO !## remove nodata --- interpolate ILIST=0; K=0; DO J=1,NCLC IF(ZDIST(J).NE.NODATA)THEN K=K+1; ILIST(J)=K; ZDIST(K)=ZDIST(J); XDIST(K)=XDIST(J) ENDIF ENDDO !## process line CALL PEUCKER_SIMPLIFYLINE(XDIST,ZDIST,GCODE,K) !## reset GCODE DO J=NCLC,1,-1 IF(ILIST(J).NE.0)THEN GCODE(J)=ABS(GCODE(ILIST(J))) ELSE GCODE(J)=0.0D0 ENDIF ENDDO IREF=ISG(I)%ICLC-1 !## never remove the first or last GCODE(1) =ZTOLERANCE+1.0D0 GCODE(NCLC)=ZTOLERANCE+1.0D0 K=1 DO J=2,NCLC K=K+1 !## remove point from Urs-Douglas-Peucker algorithm (less then given tolerance) IF(GCODE(J).LT.ZTOLERANCE)THEN CALL ISGDELISD(I,IREF+K) !## reset pointer one backwards K=K-1 ENDIF ENDDO WRITE(*,'(A,F10.2,A)') 'Progress ',REAL(I*100)/REAL(NISG),'%' ENDDO DEALLOCATE(ZDIST,XDIST,GCODE,ILIST) END SUBROUTINE ISG_SIMPLIFYMAIN !###==================================================================== LOGICAL FUNCTION ISG2GRIDMAIN(ISGFILE,IBATCH,NLAY,MLAY,TOP,BOT,KHV,BND,GRIDISG) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IBATCH,NLAY,MLAY TYPE(GRIDISGOBJ),INTENT(INOUT) :: GRIDISG CHARACTER(LEN=*),INTENT(IN) :: ISGFILE TYPE(IDFOBJ),DIMENSION(NLAY),INTENT(INOUT) :: TOP,BOT,KHV,BND TYPE(IDFOBJ),DIMENSION(1) :: SFT,VCW INTEGER :: NROW,NCOL,SIMDATE1,SIMDATE2,SIMDATE3 INTEGER,DIMENSION(1) :: MP CHARACTER(LEN=52) :: PPOSTFIX ISG2GRIDMAIN=.FALSE. IF(IBATCH.EQ.1)THEN !## deallocate memory IF(ISGREAD((/ISGFILE/),IBATCH))THEN; ENDIF ENDIF CALL ISG2GRIDGETDIMENSION(GRIDISG%IDIM,GRIDISG%XMIN,GRIDISG%YMIN,GRIDISG%XMAX,GRIDISG%YMAX,GRIDISG%CS) !,NROW,NCOL GRIDISG%IDIM=ABS(GRIDISG%IDIM) CALL UTL_IDFSNAPTOGRID_LLC(GRIDISG%XMIN,GRIDISG%XMAX,GRIDISG%YMIN,GRIDISG%YMAX,GRIDISG%CS,GRIDISG%CS,NCOL,NROW,LLC=.TRUE.) !CALL UTL_IDFSNAPTOGRID(GRIDISG%XMIN,GRIDISG%XMAX,GRIDISG%YMIN,GRIDISG%YMAX,GRIDISG%CS,NCOL,NROW) !## fill bnd(1) BND(1)%NCOL=NCOL BND(1)%NROW=NROW BND(1)%XMIN=GRIDISG%XMIN BND(1)%XMAX=GRIDISG%XMAX BND(1)%YMIN=GRIDISG%YMIN BND(1)%YMAX=GRIDISG%YMAX BND(1)%DX=GRIDISG%CS BND(1)%DY=GRIDISG%CS !## single precision BND(1)%ITYPE=4 SIMDATE1=UTL_IDATETOJDATE(GRIDISG%SDATE) !## begin juliandate SIMDATE3=UTL_IDATETOJDATE(GRIDISG%EDATE) !## eind juliandate DO IF(GRIDISG%DDATE.EQ.0)THEN SIMDATE2=SIMDATE3 ELSE SIMDATE2=SIMDATE1+MAX(0,GRIDISG%DDATE-1) ENDIF IF(SIMDATE2.GT.SIMDATE3.AND.GRIDISG%ISTEADY.EQ.2)EXIT GRIDISG%SDATE=UTL_JDATETOIDATE(SIMDATE1) GRIDISG%EDATE=UTL_JDATETOIDATE(SIMDATE2) IF(GRIDISG%DDATE.EQ.0)PPOSTFIX=GRIDISG%POSTFIX IF(GRIDISG%DDATE.NE.0)PPOSTFIX=TRIM(GRIDISG%POSTFIX)//'_'//TRIM(ITOS(GRIDISG%SDATE)) IF(IBATCH.EQ.1.AND.GRIDISG%DDATE.NE.0)WRITE(*,'(/A,I8,A,I8,A/)') 'Gridding between ',GRIDISG%SDATE,' and ',GRIDISG%EDATE,' using ppostfix:'//TRIM(PPOSTFIX) IF(GRIDISG%ISTEADY.EQ.2.AND.GRIDISG%EDATE.LT.GRIDISG%SDATE)THEN IF(IBATCH.EQ.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Enddate '//TRIM(ITOS(GRIDISG%EDATE))//' is less than startdate '// & TRIM(ITOS(GRIDISG%SDATE)),'Error') RETURN ELSE STOP 'Enddate if less than start date' ENDIF ENDIF ! CALCULATING IDF ISG2GRIDMAIN=ISG2GRID(PPOSTFIX,NROW,NCOL,NLAY,MLAY,TOP,BOT,KHV,BND,VCW,IBATCH,MP,0,GRIDISG,SFT,.FALSE.,1,1.0D0,0.0D0) IF(.NOT.ISG2GRIDMAIN)RETURN IF(GRIDISG%ISIMGRO.EQ.1)CALL ISG2GRIDMAIN_SVAT(GRIDISG) !## no multiply griddings IF(GRIDISG%DDATE.EQ.0)EXIT SIMDATE1=SIMDATE2+1 END DO END FUNCTION ISG2GRIDMAIN !###==================================================================== SUBROUTINE ISG2GRIDMAIN_SVAT(GRIDISG) !###==================================================================== IMPLICIT NONE TYPE(GRIDISGOBJ),INTENT(INOUT) :: GRIDISG TYPE(IDFOBJ) :: AHN,THIESSEN TYPE TRAPTYPE INTEGER :: IROW,ICOL REAL(KIND=DP_KIND) :: X,Y REAL(KIND=DP_KIND) :: LN,BH,BW,CT,COND_IN,COND_OUT CHARACTER(LEN=30) :: CSEGMENT END TYPE TRAPTYPE TYPE(TRAPTYPE),DIMENSION(2,20) :: TRAP INTEGER,DIMENSION(2) :: NETTRAP INTEGER :: IU,JU,ISWNR,I,J,IOS,ITRAP,ISVAT LOGICAL :: LEX REAL(KIND=DP_KIND) :: MV !## reading ahn WRITE(*,'(A)') 'Opening '//TRIM(GRIDISG%AHNFNAME)//' ...' CALL IDFNULLIFY(AHN); IF(.NOT.IDFREAD(AHN,GRIDISG%AHNFNAME,0))THEN; ENDIF WRITE(*,'(A)') 'Opening '//TRIM(GRIDISG%THIESSENFNAME)//' ...' CALL IDFNULLIFY(THIESSEN); IF(.NOT.IDFREAD(THIESSEN,GRIDISG%THIESSENFNAME,0))THEN; ENDIF CALL ISG2GRID_SEGREAD(GRIDISG) IU=UTL_GETUNIT() CALL OSD_OPEN(IU,FILE=TRIM(GRIDISG%ROOT)//'\'//TRIM(GRIDISG%SVAT2SWNR_DRNG),STATUS='OLD',FORM='FORMATTED', & ACTION='READ,DENYWRITE',ACCESS='SEQUENTIAL') JU=UTL_GETUNIT() CALL OSD_OPEN(JU,FILE=TRIM(GRIDISG%ROOT)//'\'//GRIDISG%SVAT2SWNR_DRNG(:INDEX(GRIDISG%SVAT2SWNR_DRNG,'.',.TRUE.))//'_tmp', & STATUS='UNKNOWN',FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') ISWNR=0 DO J=1,GRIDISG%NSVATS READ(IU,'(I10)',IOSTAT=IOS) NETTRAP(2) IF(IOS.NE.0)THEN; WRITE(*,*) TRIM(GRIDISG%SVAT2SWNR_DRNG)//' file does contain probably variable nettrap'; STOP; ENDIF DO ITRAP=1,NETTRAP(2) READ(IU,'(2F10.0,2I10,6F10.0,A30)',IOSTAT=IOS) TRAP(2,ITRAP)%X,TRAP(2,ITRAP)%Y,TRAP(2,ITRAP)%IROW,TRAP(2,ITRAP)%ICOL, & TRAP(2,ITRAP)%LN,TRAP(2,ITRAP)%BH,TRAP(2,ITRAP)%BW,TRAP(2,ITRAP)%CT, & TRAP(2,ITRAP)%COND_IN,TRAP(2,ITRAP)%COND_OUT,TRAP(2,ITRAP)%CSEGMENT IF(IOS.NE.0)THEN; WRITE(*,*) TRIM(GRIDISG%SVAT2SWNR_DRNG)//' file does contain probably ***"s'; STOP; ENDIF TRAP(2,ITRAP)%CSEGMENT=UTL_CAP(TRAP(2,ITRAP)%CSEGMENT,'U') TRAP(2,ITRAP)%CSEGMENT=ADJUSTL(TRAP(2,ITRAP)%CSEGMENT) END DO IF(J.GT.1)THEN LEX=.FALSE. IF(NETTRAP(2).EQ.NETTRAP(1))THEN LEX=.TRUE. DO ITRAP=1,NETTRAP(1) IF(TRAP(1,ITRAP)%IROW.NE.TRAP(2,ITRAP)%IROW) LEX=.FALSE. IF(TRAP(1,ITRAP)%ICOL.NE.TRAP(2,ITRAP)%ICOL) LEX=.FALSE. IF(TRAP(1,ITRAP)%BH.NE.TRAP(2,ITRAP)%BH) LEX=.FALSE. IF(TRAP(1,ITRAP)%BW.NE.TRAP(2,ITRAP)%BW) LEX=.FALSE. IF(TRAP(1,ITRAP)%CSEGMENT.NE.TRAP(2,ITRAP)%CSEGMENT)LEX=.FALSE. END DO IF(LEX)THEN TRAP(1,1:NETTRAP(1))%LN =TRAP(1,1:NETTRAP(1))%LN +TRAP(2,1:NETTRAP(1))%LN TRAP(1,1:NETTRAP(1))%COND_IN =TRAP(1,1:NETTRAP(1))%COND_IN +TRAP(2,1:NETTRAP(1))%COND_IN TRAP(1,1:NETTRAP(1))%COND_OUT=TRAP(1,1:NETTRAP(1))%COND_OUT+TRAP(2,1:NETTRAP(1))%COND_OUT ENDIF ENDIF IF(.NOT.LEX)THEN I=0 DO ISWNR=ISWNR+1 IF(ISWNR.GT.GRIDISG%NSWNR)THEN ISWNR=1 I =I+1 ENDIF IF(TRIM(GRIDISG%CSOBEK(ISWNR)).EQ.TRIM(TRAP(1,1)%CSEGMENT))EXIT IF(I.GT.2)EXIT END DO IF(I.GT.2)THEN WRITE(*,'(A)') 'Cannot find label '//TRIM(TRAP(1,1)%CSEGMENT)//' in '//TRIM(GRIDISG%SEGMENTCSVFNAME) ELSE DO ITRAP=1,NETTRAP(1) !## determine svat unit ISVAT=INT(IDFGETXYVAL(THIESSEN,TRAP(1,ITRAP)%X,TRAP(1,ITRAP)%Y)) !## determine surface level MV=IDFGETXYVAL(AHN,TRAP(1,ITRAP)%X,TRAP(1,ITRAP)%Y) !## translate m2/day -> day TRAP(1,ITRAP)%COND_IN =1.0D0/(TRAP(1,ITRAP)%COND_IN/(GRIDISG%CS**2.0D0)) TRAP(1,ITRAP)%COND_OUT =1.0D0/(TRAP(1,ITRAP)%COND_OUT/(GRIDISG%CS**2.0D0)) TRAP(1,ITRAP)%COND_IN =MIN(TRAP(1,ITRAP)%COND_IN,99999.99) TRAP(1,ITRAP)%COND_OUT =MIN(TRAP(1,ITRAP)%COND_OUT,99999.99) WRITE(JU,'(I10,I6,3F8.2,8X,2F8.2,8X,F8.2,8X,I10,8X,A30)') ISVAT,GRIDISG%SYSID,MV-TRAP(1,ITRAP)%BH,TRAP(1,ITRAP)%BW,TRAP(1,ITRAP)%CT, & TRAP(1,ITRAP)%LN,TRAP(1,ITRAP)%COND_IN,TRAP(1,ITRAP)%COND_OUT,ISWNR,TRAP(1,ITRAP)%CSEGMENT ENDDO ENDIF !## big change that next will be neighbourhood ISWNR=ISWNR-1 TRAP(1,:) =TRAP(2,:) NETTRAP(1)=NETTRAP(2) ENDIF ELSE NETTRAP(1)=NETTRAP(2) TRAP(1,:) =TRAP(2,:) ENDIF ENDDO CLOSE(IU,STATUS='DELETE'); CLOSE(JU) CALL IOSRENAMEFILE(TRIM(GRIDISG%ROOT)//'\'//GRIDISG%SVAT2SWNR_DRNG(:INDEX(GRIDISG%SVAT2SWNR_DRNG,'.',.TRUE.))//'_tmp', & TRIM(GRIDISG%ROOT)//'\'//TRIM(GRIDISG%SVAT2SWNR_DRNG)) END SUBROUTINE ISG2GRIDMAIN_SVAT !###====================================================================== SUBROUTINE ISG2GRID_SEGREAD(GRIDISG) !###====================================================================== IMPLICIT NONE TYPE(GRIDISGOBJ),INTENT(INOUT) :: GRIDISG INTEGER :: IU,I,IOS IU=UTL_GETUNIT(); OPEN(IU,FILE=GRIDISG%SEGMENTCSVFNAME,STATUS='OLD',ACTION='READ') READ(IU,*); GRIDISG%NSWNR=0; DO; GRIDISG%NSWNR=GRIDISG%NSWNR+1; READ(IU,*,IOSTAT=IOS); IF(IOS.NE.0)EXIT; ENDDO ALLOCATE(GRIDISG%CSOBEK(GRIDISG%NSWNR),GRIDISG%SWNR(GRIDISG%NSWNR)) REWIND(IU); READ(IU,*) DO I=1,GRIDISG%NSWNR READ(IU,*,IOSTAT=IOS) GRIDISG%SWNR(I),GRIDISG%CSOBEK(I) GRIDISG%CSOBEK(I)=UTL_CAP(GRIDISG%CSOBEK(I),'U'); GRIDISG%CSOBEK(I)=ADJUSTL(GRIDISG%CSOBEK(I)) ENDDO CLOSE(IU) END SUBROUTINE ISG2GRID_SEGREAD !###==================================================================== LOGICAL FUNCTION ISG2GRID(PPOSTFIX,NROW,NCOL,NLAY,MLAY,TOP,BOT,KHV,BND,VCW,IBATCH, & MP,JU,GRIDISG,SFT,LSFT,ISYS,FCT,IMP) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: NROW,NCOL,NLAY,MLAY,IBATCH,JU,ISYS INTEGER,INTENT(INOUT),DIMENSION(:) :: MP REAL(KIND=DP_KIND),INTENT(IN) :: FCT,IMP TYPE(IDFOBJ),DIMENSION(NLAY),INTENT(INOUT) :: TOP,BOT,KHV,BND TYPE(IDFOBJ),DIMENSION(NLAY-1),INTENT(INOUT) :: VCW TYPE(IDFOBJ),DIMENSION(:),INTENT(INOUT) :: SFT TYPE(GRIDISGOBJ),INTENT(INOUT) :: GRIDISG CHARACTER(LEN=*),INTENT(IN) :: PPOSTFIX LOGICAL,INTENT(IN) :: LSFT TYPE(WIN_MESSAGE) :: MESSAGE INTEGER :: I,J,K,II,JJ,TTIME,IROW,ICOL,NETTRAP,ITYPE,N,ISTW,IR,IC,MDIM,I1,I2,IZONE INTEGER :: JCRS,MAXNSEG,IRAT,IRAT1 INTEGER,DIMENSION(4) :: IUDMM REAL(KIND=DP_KIND) :: C,INFF,DXY,RWIDTH,WETPER,ISGLEN,AORG,ATRAP,XSTW,YSTW,GSTW,ZCHK,D,W,BEDT,BEDK,DD,CDXY REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: QSORT,RVAL REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: DIST,XNR,NDATA REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: X,Y INTEGER,ALLOCATABLE,DIMENSION(:) :: IPOS INTEGER,ALLOCATABLE,DIMENSION(:,:) :: IPSPOS REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: XIN,YIN,XSYM,YSYM,CUMDIST REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: XTRAP,YTRAP INTEGER :: ISEG,JSEG,NSEG,IREF,MAXDIM,NDIM,NSYM,NTRAP,ITRAP,JQHR,NITEMS,LU,IOS,ICRS,NCRS,IUSIMGRO,NDMM,INODE REAL(KIND=DP_KIND) :: X1,X2,Y1,Y2,CT,BH,BW,NETWD,COND,H1,H2,Z,WL,TD,TC,XC,YC LOGICAL :: LEX,LNODAT CHARACTER(LEN=256) :: LINE,TMPFNAME TYPE(IDFOBJ),ALLOCATABLE,DIMENSION(:) :: IDF TYPE(IDFOBJ) :: ICROSS,PCROSS REAL(KIND=DP_KIND),PARAMETER :: NODATAIDF=0.0D0 !## do not change !!! ISG2GRID=.FALSE. GRIDISG%NSVATS=0 IF(IBATCH.EQ.0)THEN CALL WINDOWSELECT(0); CALL WINDOWOUTSTATUSBAR(2,''); CALL WINDOWOUTSTATUSBAR(4,'Initialisation 0%') ENDIF IF(PBMAN%DMMFILE.EQ.1)THEN IUDMM(1)=UTL_GETUNIT() CALL OSD_OPEN(IUDMM(1),FILE=TRIM(GRIDISG%ROOT)//'\MFRIVTODFM1D_Q.DMM',STATUS='UNKNOWN', & FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IUDMM(1),'(3A15)') 'FM-X','FM-Y','RIV' IUDMM(2)=UTL_GETUNIT(); NDMM=0 CALL OSD_OPEN(IUDMM(2),FILE=TRIM(GRIDISG%ROOT)//'\DFM1DWATLEVTOMFRIV_H.DMM',STATUS='UNKNOWN', & FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IUDMM(2),'(4A15)') 'RIV','FM-X','FM-Y','WEIGHT' IUDMM(3)=UTL_GETUNIT(); NDMM=0 CALL OSD_OPEN(IUDMM(3),FILE=TRIM(GRIDISG%ROOT)//'\DFLOWFM_POINTS.DAT',STATUS='UNKNOWN', & FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') WRITE(IUDMM(3),'(5A15)') 'ISEG','INODE','IZONE','FM-X','FM-Y' !## read afwatidf IF(PBMAN%AFWATIDF%FNAME.NE.'')THEN CALL IDFCOPY(BND(1),PBMAN%AFWATIDF) IF(.NOT.IDFREADSCALE(PBMAN%AFWATIDF%FNAME,PBMAN%AFWATIDF,7,0,0.0D0,0))THEN WRITE(*,'(/1X,A/)') 'CANNOT READ '//TRIM(PBMAN%AFWATIDF%FNAME); STOP ENDIF ENDIF IF(JU.EQ.0)THEN IUDMM(4)=UTL_GETUNIT() IF(GRIDISG%DDATE.EQ.0)THEN CALL OSD_OPEN(IUDMM(4),FILE=TRIM(GRIDISG%ROOT)//'\modflow.riv_',STATUS='UNKNOWN',ACTION='WRITE') ELSE CALL OSD_OPEN(IUDMM(4),FILE=TRIM(GRIDISG%ROOT)//'\modflow_'//TRIM(JDATETOGDATE(GRIDISG%SDATE,2))//'.riv_', & STATUS='UNKNOWN',ACTION='WRITE') ENDIF WRITE(IUDMM(4),'(A)') 'NaN1#' ELSE IUDMM(4)=JU ENDIF ENDIF IF(GRIDISG%ISIMGRO.EQ.1)THEN CALL OSD_OPEN(IUSIMGRO,FILE=TRIM(GRIDISG%ROOT)//'\'//TRIM(GRIDISG%SVAT2SWNR_DRNG),STATUS='UNKNOWN', & FORM='FORMATTED',ACTION='WRITE,DENYREAD',ACCESS='SEQUENTIAL') ENDIF !## compute structure influences between sdate and edate IF(GRIDISG%ICDIST.EQ.1)DATISD%WL_STW=DATISD%WLVL NITEMS=MAXITEMS; IF(GRIDISG%ICDIST.EQ.0)NITEMS=9 !## open idf-filename ALLOCATE(IDF(NITEMS)) DO I=1,NITEMS CALL IDFNULLIFY(IDF(I)) CALL IDFCOPY(BND(1),IDF(I)) IDF(I)%NODATA=NODATAIDF IF(.NOT.IDFALLOCATEX(IDF(I)))THEN IF(IBATCH.EQ.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Cannot allocate memory for IDF '//TRIM(FNAME(I))//TRIM(PPOSTFIX)//'.IDF'//CHAR(13)// & 'Row/Column='//TRIM(ITOS(IDF(I)%NROW))//'-'//TRIM(ITOS(IDF(I)%NCOL)),'Error') ELSEIF(IBATCH.EQ.1)THEN WRITE(*,*) 'Cannot allocate memory for IDF '//TRIM(FNAME(I))//TRIM(PPOSTFIX)//'.IDF' WRITE(*,*) 'NROW/NCOLUMN=',IDF(I)%NCOL,IDF(I)%NROW ENDIF CALL IDFDEALLOCATE(IDF,SIZE(IDF)); DEALLOCATE(IDF) RETURN ENDIF IDF(I)%X=0.0D0 END DO !## translate cdate in to julian date - for transient simulations only! IF(GRIDISG%ISTEADY.EQ.2)THEN GRIDISG%SDATE=UTL_IDATETOJDATE(GRIDISG%SDATE) GRIDISG%EDATE=UTL_IDATETOJDATE(GRIDISG%EDATE)+1 !## take average value IF(GRIDISG%IAVERAGE.EQ.1)THEN TTIME=1 !## take median value ELSE TTIME=GRIDISG%EDATE-GRIDISG%SDATE ENDIF ELSEIF(GRIDISG%ISTEADY.EQ.1)THEN TTIME=1 ENDIF IF(ALLOCATED(QSORT))DEALLOCATE(QSORT); IF(ALLOCATED(XNR))DEALLOCATE(XNR); IF(ALLOCATED(NDATA))DEALLOCATE(NDATA) MDIM=4; IF(ISFR.EQ.1)MDIM=13; ALLOCATE(QSORT(TTIME,MDIM),XNR(MDIM),NDATA(MDIM)) NDATA=GRIDISG%NODATA 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); IF(ALLOCATED(CUMDIST))DEALLOCATE(CUMDIST) IF(ALLOCATED(IPSPOS))DEALLOCATE(IPSPOS) !## max. numbers of coordinates AND number of calculation points AND number of structures MAXNSEG=MAXVAL(ISG(1:NISG)%NSEG)+MAXVAL(ISG(1:NISG)%NCLC)+2*MAXVAL(ISG(1:NISG)%NSTW) ALLOCATE(DIST(MAXNSEG),IPOS(MAXNSEG),IPSPOS(2,MAXNSEG),RVAL(MDIM,0:MAXNSEG),X(MAXNSEG),Y(MAXNSEG),CUMDIST(MAXNSEG)) MAXDIM=0 !## compute structures IF(GRIDISG%ICDIST.EQ.1)THEN LU=UTL_GETUNIT() I =0 DO I =I+1 TMPFNAME=TRIM(PREFVAL(1))//'\tmp\tmp_stuwen_'//TRIM(ITOS(I))//'.ipf' INQUIRE(FILE=TMPFNAME,EXIST=LEX) IF(.NOT.LEX)THEN OPEN(LU,FILE=TMPFNAME,STATUS='UNKNOWN',FORM='FORMATTED',IOSTAT=IOS,ACTION='WRITE') IF(IOS.EQ.0)EXIT ENDIF ENDDO WRITE(LU,*) SUM(ISG(1:NISG)%NSTW) WRITE(LU,'(A)') '7' WRITE(LU,'(A)') 'X-COORD.' WRITE(LU,'(A)') 'Y-COORD.' WRITE(LU,'(A)') 'Z-COORD.' WRITE(LU,'(A)') 'ORIENTATION_STUWING' WRITE(LU,'(A)') 'H1-H2' WRITE(LU,'(A)') 'NUMBER_IDENT' WRITE(LU,'(A)') 'STUWNAME' WRITE(LU,'(A)') '0,TXT' ENDIF IF(IBATCH.EQ.0)THEN CALL WINDOWSELECT(0) CALL WINDOWOUTSTATUSBAR(2,'Press Escape to stop!') CALL WINDOWOUTSTATUSBAR(4,'Progress Gridding 0%') ENDIF IRAT1=0 ISTW =0 ISGLOOP: DO I=1,NISG IF(IBATCH.EQ.0)THEN CALL WMESSAGEPEEK(ITYPE,MESSAGE) IF(ITYPE.EQ.KEYDOWN.AND.MESSAGE%VALUE1.EQ.KEYESCAPE)EXIT ENDIF LEX=.TRUE. IF(GRIDISG%IDIM.EQ.3)THEN LEX=.FALSE. IF(ISG(I)%ILIST.EQ.1)LEX=.TRUE. ENDIF !## could be existing, an element with 1 coordinate (point) IF(ISG(I)%NSEG.LE.1)LEX=.FALSE. !## write nodes for coupling with dflowfmm IF(PBMAN%DMMFILE.EQ.1)THEN IZONE=0; DO INODE=1,ISG(I)%NCLC ISEG=ISG(I)%ISEG NSEG=ISG(I)%NSEG D=ISD(ISG(I)%ICLC+INODE-1)%DIST CALL ISGGETXY(ISP(ISEG:)%X,ISP(ISEG:)%Y,NSEG,D,XCRD,YCRD) !## get zone number for dfflow-fm node IF(PBMAN%AFWATIDF%FNAME.NE.'')THEN CALL IDFIROWICOL(PBMAN%AFWATIDF,IROW,ICOL,XCRD,YCRD) IZONE=PBMAN%AFWATIDF%X(ICOL,IROW) ENDIF WRITE(IUDMM(3),'(3I15,2F15.3)') I,INODE,IZONE,XCRD,YCRD ENDDO CALL IDFDEALLOCATEX(PBMAN%AFWATIDF) ENDIF IF(LEX)THEN NSEG=ISG(I)%NSEG ISEG=ISG(I)%ISEG JSEG=ISEG+NSEG-1 !## copy coordinates K=0; DO J=ISEG,JSEG; K=K+1; X(K)=ISP(J)%X; Y(K)=ISP(J)%Y; END DO CALL ISGGRIDGETSTREAMDATA(X,Y,DIST,IPOS,RVAL,ISG(I)%ICLC,ISG(I)%NCLC,ISG(I)%ISTW,ISG(I)%NSTW,NSEG,MAXNSEG, & QSORT,XNR,NDATA,TTIME,MDIM,GRIDISG%ISTEADY,GRIDISG%SDATE,GRIDISG%EDATE,GRIDISG%IAVERAGE) !## determine flow-direction LNODAT=.TRUE. DO JJ=1,4; IF(RVAL(JJ,1) .EQ.GRIDISG%NODATA)THEN; LNODAT=.FALSE. ; EXIT; ENDIF; ENDDO DO JJ=1,4; IF(RVAL(JJ,NSEG).EQ.GRIDISG%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 !## write ipf for structures IF(GRIDISG%ICDIST.EQ.1)THEN IREF=ISG(I)%ISTW+ABS(IPOS(J))-1 CALL ISG2GRIDCOMPUTE_GETXYORIENT(I,IREF,XSTW,YSTW,GSTW) IF(H2.GT.H1)GSTW=GSTW+180.0D0 GSTW=ISG2GRIDSTUWEN_GETORIENT(GSTW) ISTW=ISTW+1 LINE=TRIM(RTOS(XSTW,'G',7))//','//TRIM(RTOS(YSTW,'G',7)) //','//TRIM(RTOS(Z,'G',7))//','// & TRIM(RTOS(GSTW,'G',7))//','//TRIM(RTOS(H1-H2,'G',7))//','//TRIM(ITOS(ISTW))//',"'//TRIM(IST(IREF)%CNAME)//'"' WRITE(LU,'(A)') TRIM(LINE) ENDIF J=J+1 ENDIF ENDDO !## interpolate waterlevel,waterbottom,inf.factor,c-value - do not interupt it by structures (ipos().gt.0)! CALL ISGGRIDINTSTREAMDATA(DIST,IPOS,RVAL,NSEG,MAXNSEG,MDIM,GRIDISG%NODATA,ISG(I)%NSTW,GRIDISG%ICDIST) !## add cummulative distances CUMDIST=0.0D0 DO ISEG=2,NSEG CUMDIST(ISEG)=CUMDIST(ISEG-1)+UTL_DIST(X(ISEG-1),Y(ISEG-1),X(ISEG),Y(ISEG)) ENDDO !## get next calculation point DO ISEG=1,NSEG IF(IPOS(ISEG).GT.0)J=ISEG IPSPOS(1,ISEG)=J ENDDO !## get previous calculation point J=NSEG; IPSPOS(2,NSEG)=J DO ISEG=NSEG-1,1,-1 IF(IPSPOS(1,ISEG).NE.IPSPOS(1,ISEG+1))J=IPSPOS(1,ISEG+1) IPSPOS(2,ISEG)=J ENDDO !## start to intersect all segment/segmentpoints to the model-grid ISGLEN=0.0D0; CDXY=0.0D0 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.0D0)+((Y(ISEG)-Y(ISEG-1))**2.0D0) IF(DXY.GT.0.0D0)THEN DXY=SQRT(DXY) DO J=1,4; RVAL(J,0)=(RVAL(J,ISEG)-RVAL(J,ISEG-1))/DXY; ENDDO !## intersect line with rectangular-regular-equidistantial-grid N=0 IF(GRIDISG%NCOL.EQ.0)THEN CALL INTERSECT_EQUI(GRIDISG%XMIN,GRIDISG%XMAX,GRIDISG%YMIN,GRIDISG%YMAX, & GRIDISG%CS,GRIDISG%CS,X1,X2,Y1,Y2,N,.FALSE.) ELSE CALL INTERSECT_NONEQUI(GRIDISG%DELR,GRIDISG%DELC,GRIDISG%NROW,GRIDISG%NCOL,X1,X2,Y1,Y2,N,.FALSE.) ENDIF !## fill result array DXY=0.0D0 !## previous calculation point IF(ISEG.GT.2)THEN IF(IPSPOS(1,ISEG-1).NE.IPSPOS(1,ISEG-2))CDXY=0.0D0 ENDIF DO J=1,N ISGLEN=ISGLEN+LN(J) !## which cross-section is active within current segment CALL ISG2GRIDGETCROSS(JCRS,ISG(I)%ICRS,ISG(I)%NCRS,ISGLEN) IF(JCRS.GT.0)THEN !## start of cross-section JSEG=ISG(I)%ICRS+JCRS-1 NDIM=ABS(ISC(JSEG)%N) ELSE !## use default cross-section for current segment since no 1d cross-section or other has been found NDIM=4 ENDIF !## there will be no cross-section IF(NDIM.LE.0)CYCLE IF(NDIM*2.GT.MAXDIM)THEN MAXDIM=NDIM*2 IF(ALLOCATED(XIN))DEALLOCATE(XIN,YIN,XSYM,YSYM,XTRAP,YTRAP) ALLOCATE(XIN(MAXDIM),YIN(MAXDIM),XSYM(MAXDIM),YSYM(MAXDIM)) ALLOCATE(XTRAP(4,MAXDIM),YTRAP(4,MAXDIM)) ENDIF IF(JCRS.GT.0)THEN XIN(1:NDIM)=DATISC(ISC(JSEG)%IREF:ISC(JSEG)%IREF+NDIM-1)%DISTANCE YIN(1:NDIM)=DATISC(ISC(JSEG)%IREF:ISC(JSEG)%IREF+NDIM-1)%BOTTOM NSYM =NDIM ELSE XIN(1)=-5.0D0; XIN(2)=-5.0D0; XIN(3)=5.0D0; XIN(4)=5.0D0 YIN(1)= 5.0D0; YIN(2)= 0.0D0; YIN(3)=0.0D0; YIN(4)=0.0D0 NSYM = 4 ENDIF !## make sure cross-section its minimal z-coordinate is zero! Z=MINVAL(YIN(1:NSYM)); YIN(1:NSYM)=YIN(1:NSYM)-Z !## compute trapezia IF(GRIDISG%ISIMGRO.EQ.1)THEN CALL ISGCOMPUTETRAPEZIUM(XIN,YIN,XSYM,YSYM,XTRAP,YTRAP,NTRAP,MAXDIM,NSYM,AORG,ATRAP) IF(NSYM.LE.0)THEN IF(IBATCH.EQ.0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Crosssection: ['//TRIM(ISC(JSEG)%CNAME)//'] within: ['// & TRIM(ISG(I)%SNAME)//'] incorrect','Error') ELSE WRITE(*,*) 'Crosssection: ['//TRIM(ISC(JSEG)%CNAME)//'] within: ['// & TRIM(ISG(I)%SNAME)//'] incorrect' ENDIF EXIT ISGLOOP ENDIF ENDIF IF(LN(J).GT.0.0D0)THEN ICOL=CA(J); IROW=RA(J) !## interpolate to centre of line ... DXY = DXY+(LN(J)/2.0D0) CDXY=CDXY+(LN(J)/2.0D0) !## within model-domain IF(ICOL.GE.1 .AND.IROW.GE.1.AND. & ICOL.LE.NCOL.AND.IROW.LE.NROW)THEN ISGVALUE(1,2)=RVAL(1,ISEG-1)+DXY*RVAL(1,0) !## waterlevel ISGVALUE(1,3)=RVAL(2,ISEG-1)+DXY*RVAL(2,0) !## waterbottom IF(ISFR.EQ.0)THEN C =RVAL(3,ISEG-1)+DXY*RVAL(3,0) !## c-value ISGVALUE(1,4)=RVAL(4,ISEG-1)+DXY*RVAL(4,0) !## inf.factors ELSEIF(ISFR.EQ.1)THEN BEDT =RVAL(4,ISEG-1)+DXY*RVAL(4,0) !## river bed thickness BEDK =RVAL(5,ISEG-1)+DXY*RVAL(5,0) !## river bed permeability C =BEDT/BEDK !## c-value ISGVALUE(1,4)=1.0D0 !## inf.factors ENDIF !## factor effects conductance, therefore apply it here as reciproke IF(FCT.EQ.0.0D0)THEN; C=IMP; ELSE; C=C/FCT+IMP; ENDIF IF(C.GT.0.0D0)THEN !## translate from local to global coordinates and get proper wetted perimeter and width of channel! CALL ISG2GRIDGETPARAM(XIN,YIN,NDIM,ISGVALUE(1,3),ISGVALUE(1,2),RWIDTH,WETPER,GRIDISG%MINDEPTH) ISGVALUE(1,1)=(LN(J)*WETPER)/C !conductances ISGVALUE(1,5)= LN(J) ISGVALUE(1,6)= WETPER ISGVALUE(1,7)= RWIDTH ISGVALUE(1,8)= C IF(GRIDISG%ICDIST.EQ.1)THEN ISGVALUE(1,10)=RVAL(1,ISEG-1)+DXY*RVAL(1,0) !waterlevels ENDIF !## export in case dmmfiles are concerned IF(PBMAN%DMMFILE.EQ.1)THEN MP(1)=MP(1)+1 I1=IPSPOS(1,ISEG-1) !## next calculation point I2=IPSPOS(2,ISEG-1) !## get most close location IF(ABS(CUMDIST(I1)-CUMDIST(ISEG)).LE.ABS(CUMDIST(I2)-CUMDIST(ISEG)))THEN WRITE(IUDMM(1),'(2F15.3,I15)') X(I1),Y(I1),MP(1) ELSE WRITE(IUDMM(1),'(2F15.3,I15)') X(I2),Y(I2),MP(1) ENDIF !## fractie naar boven- en benedenstrooms knooppunt D=CUMDIST(I2)-CUMDIST(I1) DD=CDXY W=DD/D !## ignore small parts IF(W.GT.0.99D0)THEN WRITE(IUDMM(2),'(I15,3F15.3)') MP(1),X(I2),Y(I2),1.0D0 ; NDMM=NDMM+1 ELSEIF(W.LT.0.01D0)THEN WRITE(IUDMM(2),'(I15,3F15.3)') MP(1),X(I1),Y(I1),1.0D0 ; NDMM=NDMM+1 ELSE WRITE(IUDMM(2),'(I15,3F15.3)') MP(1),X(I1),Y(I1),1.0D0-W; NDMM=NDMM+1 WRITE(IUDMM(2),'(I15,3F15.3)') MP(1),X(I2),Y(I2),W ; NDMM=NDMM+1 ENDIF !## write conventional river package WRITE(IUDMM(4),'(3I10,4F10.3,I10)') 1,IROW,ICOL,ISGVALUE(1,2),ISGVALUE(1,1),ISGVALUE(1,3),ISGVALUE(1,4),ISYS ENDIF !## export data for simgro IF(GRIDISG%ISIMGRO.EQ.1)THEN !## which qhpoint is active within current segment CALL ISG2GRIDGETQHR(JQHR,ISG(I)%IQHR,ISG(I)%NQHR,ISGLEN) !## start of qhr JSEG =ISG(I)%IQHR+JQHR-1 !## get number of trapezia TC=0.0D0 DO ITRAP=1,NTRAP !## compute nett waterdepth for current part of trapezium IF(ITRAP.LT.NTRAP)THEN; NETWD=YTRAP(3,ITRAP+1)-YTRAP(3,ITRAP) ELSE; NETWD=GRIDISG%WDEPTH-YTRAP(3,ITRAP); ENDIF IF(NETWD.LE.0.0D0)EXIT !## compute wetted perimeter,bottomwidth,cotanges CALL ISG2GRIDPERIMETERTRAPEZIUM(XTRAP(:,ITRAP),YTRAP(:,ITRAP),WETPER,CT,BW,NETWD) !## conductance(m2/dag) TC=TC+(LN(J)*WETPER)/C ENDDO NETTRAP=ITRAP-1; WRITE(IUSIMGRO,'(I10)') NETTRAP !## correction factor TC=ISGVALUE(1,1)/TC DO ITRAP=1,NTRAP !## compute nett waterdepth for current part of trapezium IF(ITRAP.LT.NTRAP)THEN NETWD=YTRAP(3,ITRAP+1)-YTRAP(3,ITRAP) ELSE NETWD=GRIDISG%WDEPTH-YTRAP(3,ITRAP) ENDIF !## still water left ... IF(NETWD.GT.0.0D0)THEN !## compute wetted perimeter,bottomwidth,cotanges CALL ISG2GRIDPERIMETERTRAPEZIUM(XTRAP(:,ITRAP),YTRAP(:,ITRAP),WETPER,CT,BW,NETWD) !## get bottom height (bh) BH = ISGVALUE(1,3)+YTRAP(3,ITRAP) COND=(LN(J)*WETPER)/C !conductance(m2/dag) COND= TC*COND !## minimal value cond=0.01D0 CALL IDFGETLOC(IDF(1),IROW,ICOL,XC,YC) WRITE(IUSIMGRO,'(2F10.2,2I10,6F10.2,A30,3F10.2)') XC,YC,IROW,ICOL,LN(J),BH,BW,CT,MAX(0.01D0,COND),MAX(0.01D0,COND*ISGVALUE(1,4)), & ISQ(JSEG)%CNAME,AORG-ATRAP,AORG,ATRAP ENDIF ENDDO GRIDISG%NSVATS=GRIDISG%NSVATS+1 ENDIF !## level=conductance weighed mean IF(IDF(1)%X(ICOL,IROW).EQ.IDF(1)%NODATA)IDF(1)%X(ICOL,IROW)=0.0D0 DO II=2,NITEMS IF(IDF(II)%X(ICOL,IROW).EQ.IDF(II)%NODATA)IDF(II)%X(ICOL,IROW)=0.0D0 !## read previous level/bottom/inf.factor ISGVALUE(2,II)=IDF(II)%X(ICOL,IROW) !## multiply with conductance (5-total length) IF(II.NE.5)ISGVALUE(1,II)=ISGVALUE(1,II)*ISGVALUE(1,1) !## conductance weighed IF(ISGVALUE(2,II).NE.IDF(II)%NODATA)ISGVALUE(1,II)=ISGVALUE(2,II)+ISGVALUE(1,II) IDF(II)%X(ICOL,IROW)=ISGVALUE(1,II) END DO !## sum conductance ISGVALUE(2,1)=IDF(1)%X(ICOL,IROW) IF(ISGVALUE(2,1).NE.IDF(1)%NODATA)ISGVALUE(1,1)=ISGVALUE(2,1)+ISGVALUE(1,1) IDF(1)%X(ICOL,IROW)=ISGVALUE(1,1) ENDIF ENDIF DXY =DXY +(LN(J)/2.0D0) CDXY=CDXY+(LN(J)/2.0D0) ENDIF ENDDO ENDIF ENDDO ENDIF ENDIF IF(IBATCH.EQ.0)CALL UTL_WAITMESSAGE(IRAT,IRAT1,I,NISG,'Progress gridding ') ENDDO ISGLOOP IF(IBATCH.EQ.1)WRITE(*,'(A)') 'Finished gridding' !## divide by conductance again DO IROW=1,NROW DO ICOL=1,NCOL IF(IDF(1)%X(ICOL,IROW).NE.IDF(1)%NODATA)THEN ISGVALUE(1,1)=IDF(1)%X(ICOL,IROW) DO II=2,NITEMS !## skip length IF(II.EQ.5)CYCLE ISGVALUE(2,II)=IDF(II)%X(ICOL,IROW) ISGVALUE(1,II)=ISGVALUE(2,II)/ISGVALUE(1,1) IDF(II)%X(ICOL,IROW)=ISGVALUE(1,II) ENDDO ENDIF ENDDO ENDDO IF(PBMAN%DMMFILE.EQ.1)THEN CALL IDFWRITEFREE_HEADER(IUDMM(4),BND(1)) DO I=1,SIZE(IUDMM); IF(IUDMM(I).GT.0)CLOSE(IUDMM(I)); ENDDO ! CALL UTL_MF2005_MAXNO(TRIM(GRIDISG%ROOT)//'\MFRIVTODFM1D_Q.DMM_',(/MP(1)/)) ! CALL UTL_MF2005_MAXNO(TRIM(GRIDISG%ROOT)//'\DFM1DTOMFRIV_WL.DMM_',(/NDMM/)) ELSE !## add 2d cross-section IRAT=0; IRAT1=0 DO ISELISG=1,NISG LEX=.TRUE.; IF(GRIDISG%IDIM.EQ.3)THEN; LEX=.FALSE.; IF(ISG(ISELISG)%ILIST.EQ.1)LEX=.TRUE.; ENDIF IF(.NOT.LEX)CYCLE IF(IBATCH.EQ.0)THEN CALL WMESSAGEPEEK(ITYPE,MESSAGE) IF(ITYPE.EQ.KEYDOWN.AND.MESSAGE%VALUE1.EQ.KEYESCAPE)EXIT ENDIF NSEG=ISG(ISELISG)%NSEG; ISEG=ISG(ISELISG)%ISEG; ICRS=ISG(ISELISG)%ICRS-1; NCRS=ISG(ISELISG)%NCRS !## cross-sections IF(.NOT.ISGATTRIBUTESREADISCVALUE(0))EXIT DO I=1,NCRS ICRS=ICRS+1 !## found 2d cross-section IF(ISC(ICRS)%N.GE.0)CYCLE IF(ISC(ICRS)%N.GT.0)THEN IF(ISC(ICRS)%N.LE.2)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'iMODFLOW cannot apply a 1D cross-section with less or equal 2 points'//CHAR(13)// & 'Cross-section name '//TRIM(ISC(ICRS)%CNAME),'Error') ENDIF ELSE IF(ABS(ISC(ICRS)%N).LE.1)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'iMODFLOW cannot apply a 2D cross-section with no locations'//CHAR(13)// & 'Cross-section name '//TRIM(ISC(ICRS)%CNAME),'Error') ENDIF ENDIF CALL ISGADJUSTCOMPUTEXY(ISEG,NSEG,ISC(ICRS)%DIST,TD) !## compute correct x/y coordinate of current cross-section CALL IDFIROWICOL(IDF(2),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(I,ICROSS,PCROSS,ZCHK))THEN !## read bathymetry current cross-section WL=IDF(2)%X(ICOL,IROW) !## waterlevel at cross-section !## infiltration factor at location of cross-section INFF=IDF(4)%X(ICOL,IROW) C =IDF(8)%X(ICOL,IROW) !## resistance at location of cross-section !## intersection might miss the cell IF(C.LE.0.0D0)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(8)%X(IC,IR) IF(C.NE.0.0D0)EXIT IRLOOP ENDDO ENDDO IRLOOP ENDIF CALL ISG2GRID_BATHEMETRY(IDF,SIZE(IDF),ICROSS,PCROSS,ZCHK,WL,C,INFF) !## adjust stage grid for bathymetry CALL IDFDEALLOCATEX(ICROSS); CALL IDFDEALLOCATEX(PCROSS) ENDIF IF(IBATCH.EQ.0)CALL UTL_WAITMESSAGE(IRAT,IRAT1,ISELISG,NISG,'Progress gridding 2d cross-sections') ENDDO ENDDO IF(IBATCH.EQ.1)WRITE(*,'(A)') 'Finished gridding 2d cross-sections' ENDIF IF(ALLOCATED(DATISC2))DEALLOCATE(DATISC2); IF(ALLOCATED(TISC))DEALLOCATE(TISC); IF(ALLOCATED(ISCN))DEALLOCATE(ISCN) IF(GRIDISG%ICDIST.EQ.1)CLOSE(LU) DO IROW=1,NROW; DO ICOL=1,NCOL !## turn into nodata for conductances equal to zero IF(IDF(1)%X(ICOL,IROW).EQ.IDF(1)%NODATA)THEN IDF(2)%X(ICOL,IROW)=IDF(2)%NODATA IDF(3)%X(ICOL,IROW)=IDF(3)%NODATA IF(GRIDISG%ICDIST.EQ.1)IDF(10)%X(ICOL,IROW)=IDF(10)%NODATA ENDIF ENDDO; ENDDO IF(ALLOCATED(RVAL))DEALLOCATE(RVAL) IF(ALLOCATED(DIST))DEALLOCATE(DIST) IF(ALLOCATED(IPOS))DEALLOCATE(IPOS) IF(ALLOCATED(IPSPOS))DEALLOCATE(IPSPOS) IF(ALLOCATED(CUMDIST))DEALLOCATE(CUMDIST) IF(ALLOCATED(X))DEALLOCATE(X); IF(ALLOCATED(Y))DEALLOCATE(Y) IF(ALLOCATED(QSORT))DEALLOCATE(QSORT) IF(ALLOCATED(XNR))DEALLOCATE(XNR) IF(ALLOCATED(NDATA))DEALLOCATE(NDATA) CALL INTERSECT_DEALLOCATE() IF(GRIDISG%ISIMGRO.EQ.1)THEN IF(ALLOCATED(XIN))DEALLOCATE(XIN) IF(ALLOCATED(YIN))DEALLOCATE(YIN) IF(ALLOCATED(XSYM))DEALLOCATE(XSYM) IF(ALLOCATED(YSYM))DEALLOCATE(YSYM) IF(ALLOCATED(XTRAP))DEALLOCATE(XTRAP) IF(ALLOCATED(YTRAP))DEALLOCATE(YTRAP) GRIDISG%EDATE=GRIDISG%NSVATS ENDIF IF(PBMAN%DMMFILE.EQ.0)THEN IF(GRIDISG%ICDIST.EQ.1)THEN IF(SUM(ISG%NSTW).GT.0)THEN IF(.NOT.IDFWRITE(IDF(10),TRIM(GRIDISG%ROOT)//'\'//TRIM(FNAME(10))//TRIM(PPOSTFIX)//'.IDF',1))THEN; RETURN; ENDIF IF(.NOT.IDFWRITE(IDF(11),TRIM(GRIDISG%ROOT)//'\'//TRIM(FNAME(11))//TRIM(PPOSTFIX)//'.IDF',1))THEN; RETURN; ENDIF CALL ISG2GRIDCOMPUTESTUWEN(TRIM(TMPFNAME),TRIM(GRIDISG%ROOT)//'\'//TRIM(FNAME(10))//TRIM(PPOSTFIX)//'.IDF', & !## effect TRIM(GRIDISG%ROOT)//'\'//TRIM(FNAME(11))//TRIM(PPOSTFIX)//'.IDF') !## current_id ENDIF ENDIF !## clean data if one of the entries is nodata DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL DO I=1,4; IF(IDF(I)%X(ICOL,IROW).EQ.IDF(I)%NODATA)IDF(I)%X(ICOL,IROW)=0.0D0; ENDDO ENDDO; ENDDO !## extent grids based upon their width CALL ISG2GRID_EXTENT_WITH_WIDTH(SIZE(IDF),IDF,IBATCH,GRIDISG%MAXWIDTH) !## clean all DO IROW=1,NROW; DO ICOL=1,NCOL !## turn into nodata for conductances equal to zero IF(IDF(1)%X(ICOL,IROW).EQ.IDF(1)%NODATA)THEN DO I=1,SIZE(IDF); IDF(I)%X(ICOL,IROW)=-9999.00D0; ENDDO ENDIF ENDDO; ENDDO IDF%NODATA=-9999.0D0 IF(GRIDISG%IEXPORT.EQ.0)THEN DO I=1,9 IF(GRIDISG%ISAVE(I).EQ.0)CYCLE IF(IBATCH.EQ.1)WRITE(*,*) 'Saving '//TRIM(GRIDISG%ROOT)//'\'//TRIM(FNAME(I))//TRIM(PPOSTFIX)//'.IDF ...' IF(.NOT.IDFWRITE(IDF(I),TRIM(GRIDISG%ROOT)//'\'//TRIM(FNAME(I))//TRIM(PPOSTFIX)//'.IDF',1))THEN !##error IF(IBATCH.EQ.1)WRITE(*,*) '---- ERROR saving file ----' ENDIF ENDDO ELSEIF(GRIDISG%IEXPORT.EQ.1)THEN CALL ISG2GRID_EXPORTRIVER(JU,IDF,NLAY,MLAY,TOP,BOT,KHV,BND,VCW,MP,GRIDISG,SFT,LSFT,ISYS) ENDIF ELSE IF(JU.EQ.0)THEN IF(GRIDISG%DDATE.EQ.0)THEN CALL UTL_MF2005_MAXNO(TRIM(GRIDISG%ROOT)//'\modflow.riv_',(/MP(1)/)) ELSE CALL UTL_MF2005_MAXNO(TRIM(GRIDISG%ROOT)//'\modflow_'//TRIM(JDATETOGDATE(GRIDISG%SDATE,2))//'.riv_',(/MP(1)/)) ENDIF ENDIF ENDIF CALL IDFDEALLOCATE(IDF,SIZE(IDF)); DEALLOCATE(IDF) IF(GRIDISG%ISIMGRO.EQ.1)CLOSE(IUSIMGRO) ISG2GRID=.TRUE. END FUNCTION ISG2GRID !###==================================================================== LOGICAL FUNCTION ISG2SFR(NROW,NCOL,NLAY,ILAY,IPER,NPER,MP,JU,GRIDISG,EXFNAME,TOP,BOT,FCT,IMP) !###==================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),PARAMETER :: MAXQFLOW=1.0D5 !## error perhaps in entering discharge REAL(KIND=DP_KIND),PARAMETER :: FAREA=0.50D0 !## error in area of cross-section REAL(KIND=DP_KIND),INTENT(IN) :: FCT,IMP CHARACTER(LEN=*),INTENT(IN) :: EXFNAME INTEGER,INTENT(IN) :: NROW,NCOL,NLAY,ILAY,JU,IPER,NPER TYPE(GRIDISGOBJ),INTENT(INOUT) :: GRIDISG TYPE(IDFOBJ),DIMENSION(NLAY),INTENT(IN) :: TOP,BOT INTEGER,INTENT(INOUT),DIMENSION(:) :: MP INTEGER :: I,J,K,II,JJ,KK,KKK,TTIME,IROW,ICOL,N,ISEG,JSEG,NSEG,IREF,NDIM,KSEG,MSEG,IDATE,IISG, & ICALC,OUTSEG,IUPSEG,IPRIOR,NSTRPTS,ICRS,ICLC,IQHR,NREACH,NSTREAM,KCRS,CRSREF,KCLC,CLCREF,& NSIM,JLAY REAL(KIND=DP_KIND) :: DXY,X1,X2,Y1,Y2,QFLOW,QROFF,EVT,PREC,ROUGHCH,ROUGHBK,CDPTH,FDPTH,AWDTH,BWDTH,DIST, & HC1FCT,THICKM1,ELEVUP,WIDTH1,DEPTH1,HC2FCT,THICKM2,ELEVDN,WIDTH2,DEPTH2,WLVLUP,WLVLDN,Z, & AORG,ASIMPLE,GRAD_WLVL,GRAD_ELEV,TDIST,T,B,RB,GRAD_THCK,TH REAL(KIND=DP_KIND),PARAMETER :: MOFFSET=0.001D0 !#3 minimal offset for gradients less or equal zero REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:,:) :: QSORT,RVAL REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: XNR,NDATA REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: XCRS,ZCRS,MCRS,QCRS,WCRS,DCRS REAL(KIND=DP_KIND),DIMENSION(8) :: EXCRS,EZCRS LOGICAL :: LWARNING1,LWARNING2,LWARNING3 CHARACTER(LEN=512) :: LINE CHARACTER(LEN=MAXLENISG) :: CNAME ISG2SFR =.FALSE. IF(IPER.EQ.1)THEN LWARNING1=.TRUE.; LWARNING2=.TRUE.; LWARNING3=.TRUE. ELSE LWARNING1=.FALSE.; LWARNING2=.FALSE.; LWARNING3=.FALSE. ENDIF !## translate cdate in to julian date - for transient simulations only! IF(GRIDISG%ISTEADY.EQ.2)THEN GRIDISG%SDATE=UTL_IDATETOJDATE(GRIDISG%SDATE) GRIDISG%EDATE=UTL_IDATETOJDATE(GRIDISG%EDATE)+1 !## take average value IF(GRIDISG%IAVERAGE.EQ.1)THEN TTIME=1 !## take median value ELSE TTIME=GRIDISG%EDATE-GRIDISG%SDATE ENDIF ELSEIF(GRIDISG%ISTEADY.EQ.1)THEN TTIME=1 ENDIF IF(ALLOCATED(QSORT))DEALLOCATE(QSORT); IF(ALLOCATED(XNR))DEALLOCATE(XNR); IF(ALLOCATED(NDATA))DEALLOCATE(NDATA) NDIM=13; ALLOCATE(QSORT(TTIME,NDIM),XNR(NDIM),NDATA(NDIM)) !## allocate storage for cross-sectional data ALLOCATE(XCRS(NDIM),ZCRS(NDIM),MCRS(NDIM),QCRS(NDIM),WCRS(NDIM),DCRS(NDIM)) !## variable used to store segment information IF(ALLOCATED(RVAL))DEALLOCATE(RVAL); ALLOCATE(RVAL(NDIM,2)); RVAL=0.0 !## only specify for first stress-period - write output to regular ISG as well IF(IPER.EQ.1)THEN DO KK=1,2 IF(KK.EQ.2)THEN !## single precision (for now) ISGDOUBLE=4 !## create copy to store results ISFR=0; IF(.NOT.ISGOPENFILES(EXFNAME,'REPLACE'))THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'iMOD cannot (re)write ISG file: '//CHAR(13)//TRIM(EXFNAME),'Error') RETURN ENDIF !## write header *.ISG file in result-ISG CALL ISGSAVEHEADERS() LINE=TRIM(ITOS(IISG))//','//TRIM(ITOS(ISFR))//',Date,"StreamLevel (m+MSL)","StreamDepth (m)","StreamWidth (m)","StreamDischarge (m3/s)"' WRITE(ISGIU(1,1),'(A)') TRIM(LINE) ISFR=1 ENDIF IISG=0 KSEG=1; KCRS=1; CRSREF=1; KCLC=1; CLCREF=1 !## variable used to stored number of reaches per segment IF(ALLOCATED(ISTR))DEALLOCATE(ISTR); ALLOCATE(ISTR(NISG)); ISTR=0 !## count number of streams IF(ALLOCATED(IACTSTREAM))DEALLOCATE(IACTSTREAM) ALLOCATE(IACTSTREAM(NISG)); IACTSTREAM=0; NSTREAM=0 DO I=1,NISG NSEG=ISG(I)%NSEG; ISEG=ISG(I)%ISEG; JSEG=ISEG+NSEG-1 !## start to intersect all segment/segmentpoints to the model-grid ISTR(I)=0 N=0 !## intersect complete line DO J=ISEG+1,JSEG !## get coordinates of current reach in segment X1 =ISP(J-1)%X; Y1=ISP(J-1)%Y; X2=ISP(J)%X; Y2=ISP(J)%Y !## distance between two points with information DXY=UTL_DIST(X1,Y1,X2,Y2) !## intersect line with rectangular-regular-equidistantial-grid CALL INTERSECT_EQUI(GRIDISG%XMIN,GRIDISG%XMAX,GRIDISG%YMIN,GRIDISG%YMAX,GRIDISG%CS, & GRIDISG%CS,X1,X2,Y1,Y2,N,.FALSE.) ENDDO IF(ILAY.EQ.0)THEN !## read data for start-point IREF=ISG(I)%ICLC CALL ISG2GRIDGETDATA(GRIDISG%SDATE,GRIDISG%EDATE,TTIME,QSORT,XNR,SIZE(XNR),RVAL(1,1),ISD(IREF)%N, & ISD(IREF)%IREF,GRIDISG%ISTEADY,1,NDATA,GRIDISG%IAVERAGE) !## read data for end-point IREF=ISG(I)%ICLC+1 CALL ISG2GRIDGETDATA(GRIDISG%SDATE,GRIDISG%EDATE,TTIME,QSORT,XNR,SIZE(XNR),RVAL(1,2),ISD(IREF)%N, & ISD(IREF)%IREF,GRIDISG%ISTEADY,1,NDATA,GRIDISG%IAVERAGE) WLVLUP =RVAL(1,1); WLVLDN =RVAL(1,2) ELEVUP =RVAL(2,1); ELEVDN =RVAL(2,2) THICKM1=RVAL(4,1); THICKM2=RVAL(4,2) !## get total distance DIST=0.0D0; DO K=1,N; DIST=DIST+LN(K); ENDDO GRAD_WLVL=(WLVLDN-WLVLUP)/DIST GRAD_ELEV=(ELEVDN-ELEVUP)/DIST GRAD_THCK=(THICKM2-THICKM1)/DIST ELSE JLAY=ILAY ENDIF !## fill result array MSEG=0; TDIST=0.0D0 K=1; DO !## skip outside model domain IF(CA(K).LT.1.OR.CA(K).GT.NCOL.OR.RA(K).LT.1.OR.RA(K).GT.NROW)THEN K=K+1; TDIST=TDIST+LN(K); IF(K.GT.N)EXIT; CYCLE ENDIF !## found segment inside model IF(LN(K).GT.0.0D0)THEN !## increase number of stream reaches ISTR(I)=ISTR(I)+1 !## increase number of streams IF(ISTR(I).EQ.1)NSTREAM=NSTREAM+1 !## total number of reaches per segment IACTSTREAM(I)=IACTSTREAM(I)+1 !## total segments IISG=IISG+1 DIST=0.0D0; KKK=K; DO IF(CA(KKK).NE.CA(K).OR.RA(KKK).NE.RA(K))EXIT DIST=DIST+LN(KKK); KKK=KKK+1 ENDDO IF(KK.EQ.2)THEN ICOL=CA(K); IROW=RA(K) !## what modellayer need the SFR to be in JLAY=ILAY IF(ILAY.EQ.0)THEN RB=ELEVUP +GRAD_ELEV*(TDIST+0.5D0*DIST) TH=THICKM1+GRAD_THCK*(TDIST+0.5D0*DIST) RB=RB-TH JJ=0 DO JLAY=1,NLAY T=TOP(JLAY)%X(ICOL,IROW); B=BOT(JLAY)%X(ICOL,IROW) IF(T.NE.TOP(JLAY)%NODATA.AND.B.NE.BOT(JLAY)%NODATA)THEN JJ=JJ+1 IF(RB.LE.T.AND.RB.GE.B)EXIT ENDIF ENDDO !## put in modelayer 1 for nodes above top elevation IF(JLAY.GT.NLAY)JLAY=1 IF(JLAY.GT.NLAY)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'iMOD cannot assign element to correct modellayer'//CHAR(13)// & 'Segment '//TRIM(ITOS(I))//'; bottomlevel estimated on '//TRIM(RTOS(RB,'F',2)),'Error') RETURN ENDIF ENDIF !## write into sfr package LINE=TRIM(ITOS(JLAY)) //','//TRIM(ITOS(IROW)) //','//TRIM(ITOS(ICOL))//','// & TRIM(ITOS(NSTREAM))//','//TRIM(ITOS(ISTR(I)))//','//TRIM(RTOS(DIST,'G',7)) WRITE(JU,'(A)') TRIM(LINE) !## write coordinate-couple in isp (use all coordinates) DO JJ=K,KKK-1 KSEG=KSEG+1; MSEG=MSEG+1 IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(2,1),REC=KSEG) REAL(XA(JJ),4) ,REAL(YA(JJ) ,4) IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(2,1),REC=KSEG) XA(JJ) , YA(JJ) KSEG=KSEG+1; MSEG=MSEG+1 IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(2,1),REC=KSEG) REAL(XA(JJ+1),4),REAL(YA(JJ+1),4) IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(2,1),REC=KSEG) XA(JJ+1) , YA(JJ+1) ENDDO !## add cross-section ICRS=ISG(I)%ICRS !## position in isc that starts cross-section NDIM=ABS(ISC(ICRS)%N) !## number of cross-sectional data-points KCRS=KCRS+1 IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(5,1),REC=KCRS) ISC(ICRS)%N,CRSREF,REAL(0.5*DIST,4),ISC(ICRS)%CNAME IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(5,1),REC=KCRS) ISC(ICRS)%N,CRSREF, 0.5D0*DIST ,ISC(ICRS)%CNAME !## write cross-section II=ISC(ICRS)%IREF-1 ; DO JJ=1,NDIM II=II+1; CRSREF=CRSREF+1 IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(6,1),REC=CRSREF) REAL(DATISC(II)%DISTANCE,4),REAL(DATISC(II)%BOTTOM,4),REAL(DATISC(II)%MRC,4) IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(6,1),REC=CRSREF) DATISC(II)%DISTANCE , DATISC(II)%BOTTOM , DATISC(II)%MRC ENDDO ICLC=ISG(I)%ICLC !## position in isd that starts calculation node !## add space for time-variant data - put data into it DO II=1,2 KCLC=KCLC+1 IF(II.EQ.1)THEN CNAME='From_Node' IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(3,1),REC=KCLC) NPER,CLCREF,0.0 ,CNAME IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(3,1),REC=KCLC) NPER,CLCREF,0.0D0,CNAME ELSEIF(II.EQ.2)THEN CNAME='To_Node' IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(3,1),REC=KCLC) NPER,CLCREF,REAL(DIST,4),CNAME IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(3,1),REC=KCLC) NPER,CLCREF,DIST ,CNAME ENDIF DO JJ=1,NPER CLCREF=CLCREF+1 IF(SIM(JJ)%DELT.GT.0.0D0)THEN IDATE=SIM(JJ)%IYR*10000+SIM(JJ)%IMH*100+SIM(JJ)%IDY ELSE IF(NPER.EQ.1)THEN IDATE=UTL_GETCURRENTDATE() ELSE IDATE=JD(SIM(2)%IYR,SIM(2)%IMH,SIM(2)%IDY) IDATE=UTL_JDATETOIDATE(IDATE-1) ENDIF ENDIF IF(ISGDOUBLE.EQ.4)WRITE(ISGIU(4,1),REC=CLCREF) IDATE,10.0 ,10.0 ,10.0 IF(ISGDOUBLE.EQ.8)WRITE(ISGIU(4,1),REC=CLCREF) IDATE,10.0D0,10.0D0,10.0D0 ENDDO ENDDO !## each reach is a segment in the isg-file LINE='"'//TRIM(ISG(I)%SNAME)//'_'//TRIM(ITOS(NSTREAM))//'-'//TRIM(ITOS(ISTR(I)))//'",'//TRIM(ITOS(KSEG-MSEG))//','// & TRIM(ITOS(MSEG)) //','//TRIM(ITOS(KCLC-2))//','//TRIM(ITOS(2))//','// & TRIM(ITOS(KCRS-1)) //','//TRIM(ITOS(1))//',0,0,0,0' WRITE(ISGIU(1,1),'(A)') TRIM(LINE) MSEG=0 ENDIF TDIST=TDIST+DIST K=KKK ELSE K=K+1 ENDIF IF(K.GT.N)EXIT ENDDO ENDDO IF(KK.EQ.1)DEALLOCATE(ISTR,IACTSTREAM) ENDDO ENDIF !## total number of reaches - determines for stress-period 1, stays the same NREACH=SUM(ISTR) !## to be filled in later - number of streams LINE='NaN1#,'//TRIM(ITOS(IRDFLG))//','//TRIM(ITOS(IPTFLG))//',0' WRITE(JU,'(A)') TRIM(LINE) !## write cross-sectional data NSTREAM=0; DO I=1,NISG !## skip this if not in current model domain IF(IACTSTREAM(I).LE.0)CYCLE !## subsequent numbering of segments NSTREAM=NSTREAM+1 !## get this information from the cross-sectional data (only one cross-section available) ICRS=ISG(I)%ICRS !## position in isc that starts cross-section NDIM=ABS(ISC(ICRS)%N) !## number of cross-sectional data-points IF(NDIM.GT.SIZE(XCRS))THEN DEALLOCATE(XCRS,ZCRS,MCRS) ALLOCATE(XCRS(NDIM),ZCRS(NDIM),MCRS(NDIM)) ENDIF J=ISC(ICRS)%IREF-1 ; DO K=1,NDIM J=J+1; XCRS(K)=DATISC(J)%DISTANCE; ZCRS(K)=DATISC(J)%BOTTOM; MCRS(K)=DATISC(J)%MRC ENDDO RVAL=0.0D0 !## read data for start-point IREF=ISG(I)%ICLC CALL ISG2GRIDGETDATA(GRIDISG%SDATE,GRIDISG%EDATE,TTIME,QSORT,XNR,SIZE(XNR),RVAL(1,1),ISD(IREF)%N, & ISD(IREF)%IREF,GRIDISG%ISTEADY,1,NDATA,GRIDISG%IAVERAGE) !## read data for end-point IREF=ISG(I)%ICLC+1 CALL ISG2GRIDGETDATA(GRIDISG%SDATE,GRIDISG%EDATE,TTIME,QSORT,XNR,SIZE(XNR),RVAL(1,2),ISD(IREF)%N, & ISD(IREF)%IREF,GRIDISG%ISTEADY,1,NDATA,GRIDISG%IAVERAGE) WLVLUP =RVAL(1,1); WLVLDN =RVAL(1,2) ELEVUP =RVAL(2,1); ELEVDN =RVAL(2,2) WIDTH1 =RVAL(3,1); WIDTH2 =RVAL(3,2) THICKM1=RVAL(4,1); THICKM2=RVAL(4,2) HC1FCT =RVAL(5,1); HC2FCT =RVAL(5,2) !## factor effects conductance, therefore apply it here as reciproke HC1FCT=HC1FCT*FCT+IMP HC2FCT=HC2FCT*FCT+IMP ICALC =INT(RVAL(8,1)) !## calculation option streamdepth IUPSEG=INT(RVAL(6,1)) !## upstream segment OUTSEG=INT(RVAL(7,2)) !## downstream segment IPRIOR=INT(RVAL(9,1)) !## dividing option QFLOW =RVAL(10,1) !## inflow QROFF =RVAL(11,1) !## runoff flow !## bottom level up always larger than bottom level down IF(ELEVUP.LE.ELEVDN)THEN IF(LWARNING3)THEN CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'You entered an upstream bottom elevation ('//TRIM(RTOS(ELEVUP,'G',7))//') for segment '//trim(itos(I))//CHAR(13)// & 'which is lower/equal to the bottom elevation ('//TRIM(RTOS(ELEVDN,'G',7))//') downstream.'//CHAR(13)//'The SFR package does not allow this at any time.'//CHAR(13)// & 'Do you want to add an offset of '//TRIM(RTOS(MOFFSET,'F',4))//' for all that are equal and continue?','Question') IF(WINFODIALOG(4).NE.1)RETURN LWARNING3=.FALSE. ENDIF ELEVDN=ELEVUP-MOFFSET ENDIF !## check qflow bigger than maxqflow; m3/s might be entered as m3/d IF(QFLOW.GT.MAXQFLOW.AND.LWARNING1)THEN CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'You might entered the external discharge in m3/d instead of m3/s.'//CHAR(13)// & 'iMOD found a value of '//TRIM(RTOS(QFLOW,'*',5))//CHAR(13)// & 'That might be wrong'//CHAR(13)//'Do you want to continue ?','Question') IF(WINFODIALOG(4).NE.0)RETURN LWARNING1=.FALSE. ENDIF !## translate inflow from m3/s to m3/d - only if not a diversion IF(IUPSEG.EQ.0)QFLOW=QFLOW*86400.0D0 !## translate surface runoff from m3/s to m3/d QROFF=QROFF*86400.0D0 !## check if lake, whether the lake number is still correct CALL ISG2SFR_GETLAKENUMBER(IUPSEG,I) CALL ISG2SFR_GETLAKENUMBER(OUTSEG,I) !## corrections for reading out of a menu ICALC=ICALC-1; IPRIOR=-1*(IPRIOR-1) EVT =RVAL(12,1)/1000.0D0 PREC=RVAL(13,1)/1000.0D0 LINE=TRIM(ITOS(NSTREAM))//','//TRIM(ITOS(ICALC))//','//TRIM(ITOS(OUTSEG))//','//TRIM(ITOS(IUPSEG)) IF(IUPSEG.GT.0)LINE=TRIM(LINE)//','//TRIM(ITOS(IPRIOR)) !## IF(ICALC.EQ.4)THEN !## get this information from the q-width/depth relationships data (only one q-width/depth relationships available) IQHR=ISG(I)%IQHR !## position in isq that starts q-width/depth relationships NSTRPTS=ABS(ISQ(IQHR)%N) !## number of q-width/depth relationships LINE=TRIM(LINE)//','//TRIM(ITOS(NSTRPTS)) ENDIF LINE=TRIM(LINE)//','//TRIM(RTOS(QFLOW,'G',7))//','//TRIM(RTOS(QROFF,'G',7))//','// & TRIM(RTOS(EVT,'G',7))//','//TRIM(RTOS(PREC,'G',7)) !## riverbed mannings coefficient IF(ICALC.EQ.1.OR.ICALC.EQ.2)THEN ROUGHCH=MCRS(4) LINE=TRIM(LINE)//','//TRIM(RTOS(ROUGHCH,'G',7)) ENDIF !## riverbank mannings coefficient IF(ICALC.EQ.2)THEN ROUGHBK=MCRS(1) LINE=TRIM(LINE)//','//TRIM(RTOS(ROUGHBK,'G',7)) ENDIF !## function IF(ICALC.EQ.3)THEN !## not supported CDPTH=0.3D0 FDPTH=0.35D0 AWDTH=3.8D0 BWDTH=0.6D0 LINE=TRIM(LINE)//','//TRIM(RTOS(CDPTH,'G',7))//','//TRIM(RTOS(FDPTH,'G',7))//','// & TRIM(RTOS(AWDTH,'G',7))//','//TRIM(RTOS(BWDTH,'G',7)) ENDIF WRITE(JU,'(A)') TRIM(LINE) !## waterdepth up- and downstream DEPTH1 =WLVLUP-ELEVUP DEPTH2 =WLVLDN-ELEVDN LINE=TRIM(RTOS(HC1FCT,'G',7))//','//TRIM(RTOS(THICKM1,'G',7))//','//TRIM(RTOS(ELEVUP,'G',7)) IF(ICALC.LE.1)THEN IF(WIDTH1.LE.0.0D0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You cannot asign a width of '//TRIM(RTOS(WIDTH1,'F',2))//CHAR(13)// & ' for segment '//TRIM(ITOS(I))//' called'//CHAR(13)//TRIM(ISG(I)%SNAME),'Error') RETURN ENDIF LINE=TRIM(LINE)//','//TRIM(RTOS(WIDTH1,'G',7)) ENDIF IF(ICALC.EQ.0)LINE=TRIM(LINE)//','//TRIM(RTOS(DEPTH1,'G',7)) WRITE(JU,'(A)') TRIM(LINE) LINE=TRIM(RTOS(HC2FCT,'G',7))//','//TRIM(RTOS(THICKM2,'G',7))//','//TRIM(RTOS(ELEVDN,'G',7)) IF(ICALC.LE.1)THEN LINE=TRIM(LINE)//','//TRIM(RTOS(WIDTH2,'G',7)) IF(WIDTH2.LE.0.0D0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'You cannot asign a width of '//TRIM(RTOS(WIDTH2,'F',2))//CHAR(13)// & ' for segment '//TRIM(ITOS(I))//' called'//CHAR(13)//TRIM(ISG(I)%SNAME),'Error') RETURN ENDIF ENDIF IF(ICALC.EQ.0)LINE=TRIM(LINE)//','//TRIM(RTOS(DEPTH2,'G',7)) WRITE(JU,'(A)') TRIM(LINE) !## eight points cross-section IF(ICALC.EQ.2)THEN !## try to get eight-point cross-section NSIM=8; CALL ISGCOMPUTEEIGHTPOINTS(XCRS,ZCRS,NDIM,EXCRS,EZCRS,NSIM,AORG,ASIMPLE) IF(LWARNING2.AND.(ASIMPLE/AORG.LT.FAREA.OR.ASIMPLE/AORG.GT.1.0D0/FAREA))THEN CALL WMESSAGEBOX(YESNO,QUESTIONICON,COMMONNO,'For cross-section number '//TRIM(ITOS(I))//CHAR(13)// & 'the area of the estimated cross-section with eight points ('// & TRIM(RTOS(ASIMPLE,'*',5))//')'//CHAR(13)//'differs significantly with the area of the original'//CHAR(13)// & 'cross-section ('//TRIM(RTOS(AORG,'*',5))//').'//CHAR(13)//'Do you want to continue ?','Question') IF(WINFODIALOG(4).NE.1)RETURN LWARNING2=.FALSE. ENDIF !## make excrs relative Z=EXCRS(1); DO J=1,NSIM; EXCRS(J)=EXCRS(J)-Z; ENDDO DO J=1,NSIM IF(EXCRS(J).LT.0.0D0)THEN CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'Check your cross-section for segment '//TRIM(ITOS(I))//CHAR(13)// & 'iMOD founds that the construction of this cross-section is invalid','Error') RETURN ENDIF ENDDO !## make ezcrs relative Z=MINVAL(EZCRS); DO J=1,NSIM; EZCRS(J)=EZCRS(J)-Z; ENDDO LINE=TRIM(RTOS(EXCRS(1),'G',7)); DO J=2,NSIM; LINE=TRIM(LINE)//','//TRIM(RTOS(EXCRS(J),'G',7)); ENDDO WRITE(JU,'(A)') TRIM(LINE) LINE=TRIM(RTOS(EZCRS(1),'G',7)); DO J=2,NSIM; LINE=TRIM(LINE)//','//TRIM(RTOS(EZCRS(J),'G',7)); ENDDO WRITE(JU,'(A)') TRIM(LINE) !## q-width/depth relationships ELSEIF(ICALC.EQ.4)THEN !## get this information from the q-width/depth relationships data IQHR=ISG(I)%IQHR IF(NSTRPTS.GT.SIZE(QCRS))THEN DEALLOCATE(QCRS,DCRS,WCRS) ALLOCATE(QCRS(NSTRPTS),DCRS(NSTRPTS),WCRS(NSTRPTS)) ENDIF J=ISQ(IQHR)%IREF-1 ; DO K=1,NSTRPTS J=J+1; QCRS(K)=DATISQ(J)%Q*86400.0D0; WCRS(K)=DATISQ(J)%W; DCRS(K)=DATISQ(J)%D ENDDO LINE=TRIM(RTOS(QCRS(1),'G',7)); DO J=2,NSTRPTS; LINE=TRIM(LINE)//','//TRIM(RTOS(QCRS(J),'G',7)); ENDDO WRITE(JU,'(A)') TRIM(LINE) LINE=TRIM(RTOS(DCRS(1),'G',7)); DO J=2,NSTRPTS; LINE=TRIM(LINE)//','//TRIM(RTOS(DCRS(J),'G',7)); ENDDO WRITE(JU,'(A)') TRIM(LINE) LINE=TRIM(RTOS(WCRS(1),'G',7)); DO J=2,NSTRPTS; LINE=TRIM(LINE)//','//TRIM(RTOS(WCRS(J),'G',7)); ENDDO WRITE(JU,'(A)') TRIM(LINE) ENDIF ENDDO MP(1)=NSTREAM MP(2)=NREACH IF(ALLOCATED(XCRS))DEALLOCATE(XCRS,ZCRS,MCRS) IF(ALLOCATED(QCRS))DEALLOCATE(QCRS,WCRS,DCRS) IF(ALLOCATED(RVAL))DEALLOCATE(RVAL) IF(ALLOCATED(QSORT))DEALLOCATE(QSORT); IF(ALLOCATED(XNR))DEALLOCATE(XNR); IF(ALLOCATED(NDATA))DEALLOCATE(NDATA) !## clean up IF(IPER.EQ.NPER)THEN DEALLOCATE(ISTR,IACTSTREAM) ENDIF CALL INTERSECT_DEALLOCATE() ISG2SFR=.TRUE. END FUNCTION ISG2SFR !###==================================================================== SUBROUTINE ISG2SFR_GETLAKENUMBER(ISEG,ISTR) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(INOUT) :: ISEG INTEGER,INTENT(IN) :: ISTR INTEGER :: I !## nothing to do - no lake IF(ISEG.GE.0)RETURN !## remove connection if no lakes available IF(.NOT.ALLOCATED(ULAKES))THEN !## turn off connection ISEG=0; WRITE(*,'(/A)') 'Error cannot connect Stream '//TRIM(ITOS(ISTR))//' as Lake package is inactive' WRITE(*,'(A/)') 'iMOD removed the connection to continue' RETURN ENDIF !## find appropriate lake DO I=1,SIZE(ULAKES) IF(ULAKES(I).EQ.ABS(ISEG))THEN !## get appropriate lake number - negative as it is a lake ISEG=-I RETURN ENDIF ENDDO WRITE(*,'(/A)') 'Error cannot connect Stream '//TRIM(ITOS(ISTR))//' to lake '//TRIM(ITOS(ISEG)) WRITE(*,'(A)') 'Perhaps you used an updated LAKE-ID file and did not update connection ID in SFR file?' WRITE(*,'(A/)') 'iMOD removed the connection to continue' ISEG=0; RETURN END SUBROUTINE ISG2SFR_GETLAKENUMBER !###==================================================================== SUBROUTINE ISG2GRID_EXPORTRIVER(JU,IDF,NLAY,MLAY,TOP,BOT,KHV,BND,VCW,MP,GRIDISG,SFT,LSFT,ISYS) !###==================================================================== IMPLICIT NONE TYPE(GRIDISGOBJ),INTENT(INOUT) :: GRIDISG INTEGER,INTENT(IN) :: NLAY,MLAY,JU,ISYS INTEGER,INTENT(INOUT),DIMENSION(:) :: MP TYPE(IDFOBJ),DIMENSION(:),INTENT(INOUT) :: IDF,SFT,VCW TYPE(IDFOBJ),DIMENSION(NLAY),INTENT(INOUT) :: TOP,BOT,KHV,BND LOGICAL,INTENT(IN) :: LSFT INTEGER :: IROW,ICOL,N,IU,I,J REAL(KIND=DP_KIND) :: T,WL,BL,CD,F,TL,WP,LC,C1,C0,DXY,D,K CHARACTER(LEN=25) :: FRM REAL(KIND=DP_KIND),ALLOCATABLE,DIMENSION(:) :: TLP,KH,TP,BT REAL(KIND=DP_KIND),PARAMETER :: MINKHT=0.0D0 ! INTEGER,PARAMETER :: ICLAY=1 !## shift to nearest aquifer IF(MLAY.EQ.0)THEN !## read in all top/bottom layers DO I=1,NLAY IF(.NOT.ASSOCIATED(TOP(I)%X))THEN WRITE(*,'(A)') 'Reading '//TRIM(TOP(I)%FNAME)//' ...' CALL IDFCOPY(IDF(1),TOP(I)) IF(.NOT.IDFREADSCALE(TOP(I)%FNAME,TOP(I),2,1,0.0D0,0))RETURN ENDIF IF(.NOT.ASSOCIATED(BOT(I)%X))THEN WRITE(*,'(A)') 'Reading '//TRIM(BOT(I)%FNAME)//' ...' CALL IDFCOPY(IDF(1),BOT(I)) IF(.NOT.IDFREADSCALE(BOT(I)%FNAME,BOT(I),2,1,0.0D0,0))RETURN ENDIF ENDDO ENDIF WRITE(FRM,'(A9,I2.2,A14)') '(3(I5,1X),',4,'(G15.7,1X),I5)' !## idf(1)=cond !## idf(2)=stage !## idf(3)=bottom !## idf(4)=inffct ALLOCATE(TLP(NLAY),KH(NLAY),TP(NLAY),BT(NLAY)) DO J=1,2 N=0 DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL !## no element on this location IF(IDF(1)%X(ICOL,IROW).LE.0.0D0)CYCLE !## assign layers in between waterlevel and bottom elevation IF(PBMAN%IDEFLAYER.EQ.0)THEN WL=IDF(2)%X(ICOL,IROW) ELSE !## assign layers in between top modellayer 1 and bottom elevation WL=TOP(1)%X(ICOL,IROW) ENDIF BL=IDF(3)%X(ICOL,IROW); CD=IDF(1)%X(ICOL,IROW); F =IDF(4)%X(ICOL,IROW) TL=IDF(5)%X(ICOL,IROW); WP=IDF(6)%X(ICOL,IROW); C0=IDF(8)%X(ICOL,IROW) TLP=0.0D0 IF(MLAY.EQ.0)THEN DO I=1,NLAY TP(I)=TOP(I)%X(ICOL,IROW) BT(I)=BOT(I)%X(ICOL,IROW) KH(I)=KHV(I)%X(ICOL,IROW) T =TP(I)-BT(I) !## apply minimal transmissivity IF(BND(I)%X(ICOL,IROW).GT.0.AND.T.GT.0.0D0)THEN IF(KH(I)*T.LT.PBMAN%MINKD)KH(I)=PBMAN%MINKD/T ELSE KH(I)=0.0D0 ENDIF ENDDO CALL UTL_PCK_GETTLP(NLAY,TLP,KH,TP,BT,WL,BL,MINKHT) ELSE TLP(MLAY)=1.0D0; DO I=1,NLAY; KH(I)=KHV(I)%X(ICOL,IROW) ; ENDDO ENDIF !## assign to modellayer 1 by default, when no layers are read in DO I=1,NLAY IF(TLP(I).EQ.0.0D0)CYCLE !## skip inactive cells IF(BND(I)%X(ICOL,IROW).LE.0.0D0)CYCLE N=N+1 IF(J.EQ.2)THEN LC=1.0D0 IF(PBMAN%IFVDL.EQ.1)THEN DXY=IDFGETAREA(IDF(1),ICOL,IROW) C1 =VCW(MLAY)%X(ICOL,IROW) IF(C1.GT.0.0D0)THEN !## if sft package is active use that parameter IF(LSFT)THEN D=SFT(1)%X(ICOL,IROW) K=SFT(2)%X(ICOL,IROW) ELSE D=TOP(MLAY)%X(ICOL,IROW)-BOT(MLAY)%X(ICOL,IROW) K=KH(I) ENDIF LC=ISG_UTL_LEAKAGE(CD,DXY,F,TL,WP,C0,C1,D,K) ELSE LC=0.0D0 ENDIF ENDIF WRITE(IU,FRM) I,IROW,ICOL,WL,LC*CD*TLP(I),BL,F,ISYS ENDIF ENDDO ENDDO; ENDDO IF(J.EQ.1)THEN MP(1)=MP(1)+N IF(JU.EQ.0)THEN IU=UTL_GETUNIT() IF(GRIDISG%DDATE.EQ.0)THEN CALL OSD_OPEN(IU,FILE=TRIM(GRIDISG%ROOT)//'\modflow.riv',STATUS='UNKNOWN',ACTION='WRITE') WRITE(IU,'(A)') 'SCD 1 ' ELSE CALL OSD_OPEN(IU,FILE=TRIM(GRIDISG%ROOT)//'\modflow_'//TRIM(JDATETOGDATE(GRIDISG%SDATE,2))//'.riv', & STATUS='UNKNOWN',ACTION='WRITE') WRITE(IU,'(A)') 'SCD 1 '//TRIM(JDATETOGDATE(GRIDISG%SDATE,2))//'-'//TRIM(JDATETOGDATE(GRIDISG%EDATE,2)) ENDIF WRITE(IU,'(I10)') N ELSE IU=JU ENDIF ENDIF ENDDO !## not for sfr IF(PBMAN%IFORMAT.EQ.2)CALL IDFWRITEFREE_HEADER(IU,BND(1)) CLOSE(IU) DEALLOCATE(TLP,KH,TP,BT) END SUBROUTINE ISG2GRID_EXPORTRIVER !###==================================================================== SUBROUTINE ISG2GRID_BATHEMETRY(IDF,NIDF,ICROSS,PCROSS,ZCHK,WL,C,INFF) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: NIDF TYPE(IDFOBJ),DIMENSION(NIDF),INTENT(INOUT) :: IDF TYPE(IDFOBJ),INTENT(IN) :: ICROSS,PCROSS REAL(KIND=DP_KIND),INTENT(IN) :: WL,C,INFF,ZCHK INTEGER :: IROW,ICOL,JROW,JCOL REAL(KIND=DP_KIND) :: 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.0D0.AND.WL.LE.ZCHK)CYCLE !## all okay - continue CALL IDFGETLOC(ICROSS,IROW,ICOL,XC,YC) CALL IDFIROWICOL(IDF(1),JROW,JCOL,XC,YC) !## outside current model IF(JROW.EQ.0.OR.JCOL.EQ.0)CYCLE !## manipulate resistance for ... CR=C !## skip location with c=0.0 IF(CR.EQ.0.0D0)CYCLE IF(PCROSS%X(ICOL,IROW).NE.0.0D0)CR=CR*ABS(PCROSS%X(ICOL,IROW)) IF(IDF(9)%X(JCOL,JROW).EQ.0.0D0)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.0D0 !## 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,IBATCH,MAXWIDTH) !###==================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),PARAMETER :: FDX=1.5 !## apply extent whenever width is 1.5*dx cell INTEGER,INTENT(IN) :: NIDF,IBATCH REAL(KIND=DP_KIND),INTENT(IN) :: MAXWIDTH TYPE(IDFOBJ),INTENT(INOUT),DIMENSION(NIDF) :: IDF INTEGER :: IROW,ICOL,I,IRAT,IRAT1 REAL(KIND=DP_KIND) :: W,L,A DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL !## already visited by bathymetry routine - so skip it IF(IDF(9)%X(ICOL,IROW).NE.0.0D0)THEN IDF(9)%X(ICOL,IROW)=-1.0D0*IDF(9)%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(9)%X(ICOL,IROW)=1.0D0 !## 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(5)%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 IRAT=0; IRAT1=0 DO IROW=1,IDF(1)%NROW; DO ICOL=1,IDF(1)%NCOL !## already processed by bathymetry IF(IDF(9)%X(ICOL,IROW).LT.0.0D0)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; IF(IBATCH.EQ.0)CALL UTL_WAITMESSAGE(IRAT,IRAT1,IROW,IDF(1)%NROW,'Progress computing erosion matrix'); ENDDO IF(IBATCH.EQ.1)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.0D0)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) !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: ICOL,IROW,NIDF TYPE(IDFOBJ),INTENT(INOUT),DIMENSION(NIDF) :: IDF REAL(KIND=DP_KIND),INTENT(IN) :: W,A INTEGER :: I,K,IC,IR,IC1,IC2,IR1,IR2 REAL(KIND=DP_KIND) :: X,Y,X1,X2,Y1,Y2,TRAD,C,F,DX,DY REAL(KIND=DP_KIND),DIMENSION(9) :: XP,YP !## half of the river width TRAD=W/2.0 !## get current location CALL IDFGETLOC(IDF(1),IROW,ICOL,X,Y) !## get window for brush-function CALL IDFIROWICOL(IDF(1),IR2,IC1,X-TRAD,Y-TRAD) CALL IDFIROWICOL(IDF(1),IR1,IC2,X+TRAD,Y+TRAD) !## store conductance,stage,bottom,entry-restistance C=IDF(1)%X(ICOL,IROW)/A DO IR=MAX(1,IR1),MIN(IR2,IDF(1)%NROW) DO IC=MAX(1,IC1),MIN(IC2,IDF(1)%NCOL) !## skip already filled in bathymetry/or previous wide profile IF(IDF(9)%X(IC,IR).LT.0.0D0.OR.IDF(9)%X(IC,IR).GE.1.0D0)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.0D0 DY=(YP(I)-Y)**2.0D0 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.0D0.OR.F.GT.IDF(9)%X(IC,IR))THEN IDF(9)%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 !###==================================================================== SUBROUTINE ISG2GRIDGETDIMENSION(IDIM,XMIN,YMIN,XMAX,YMAX,CS) !,NROW,NCOL !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: IDIM REAL(KIND=DP_KIND),INTENT(IN) :: CS REAL(KIND=DP_KIND),INTENT(OUT) :: XMIN,YMIN,XMAX,YMAX ! INTEGER,INTENT(OUT) :: NROW,NCOL INTEGER :: I,J,NSEG,ISEG,JSEG LOGICAL :: LEX SELECT CASE (IDIM) CASE (-1,-2,-3,1) !## current domain XMIN=MPW%XMIN YMIN=MPW%YMIN XMAX=MPW%XMAX YMAX=MPW%YMAX CASE (2,3) !## entire isg, current selection within isg J=0 DO I=1,NISG LEX=.TRUE. IF(IDIM.EQ.3)THEN LEX=.FALSE. IF(ISG(I)%ILIST.EQ.1)LEX=.TRUE. ENDIF IF(LEX)THEN NSEG=ISG(I)%NSEG ISEG=ISG(I)%ISEG JSEG=ISEG+NSEG-1 J =J+1 IF(J.EQ.1)THEN XMIN=MINVAL(ISP(ISEG:JSEG)%X) XMAX=MAXVAL(ISP(ISEG:JSEG)%X) YMIN=MINVAL(ISP(ISEG:JSEG)%Y) YMAX=MAXVAL(ISP(ISEG:JSEG)%Y) ELSE XMIN=MIN(XMIN,MINVAL(ISP(ISEG:JSEG)%X)) XMAX=MAX(XMAX,MAXVAL(ISP(ISEG:JSEG)%X)) YMIN=MIN(YMIN,MINVAL(ISP(ISEG:JSEG)%Y)) YMAX=MAX(YMAX,MAXVAL(ISP(ISEG:JSEG)%Y)) ENDIF ENDIF ENDDO !## increase it a little bit to make sure everything is captured XMIN=XMIN-CS XMAX=XMAX+CS YMIN=YMIN-CS YMAX=YMAX+CS END SELECT END SUBROUTINE ISG2GRIDGETDIMENSION !###====================================================================== SUBROUTINE ISG2GRIDCOMPUTESTUWEN(IPFFNAME,IDFFNAME1,IDFFNAME2) !###====================================================================== IMPLICIT NONE CHARACTER(LEN=*),INTENT(IN) :: IPFFNAME,IDFFNAME1,IDFFNAME2 REAL(KIND=DP_KIND) :: XC,YC,Z,ORIENT,DH INTEGER :: I,NSTUW,NC,IROW,ICOL,NTHREAD,DTERM,MAXTHREAD,J,IMENU,IRAT,IRAT1,IU,IOS,MAXN INTEGER(KIND=1),POINTER,DIMENSION(:) :: ISPEC INTEGER(KIND=2),POINTER,DIMENSION(:,:) :: THREAD,YSEL TYPE(IDFOBJ),DIMENSION(:),ALLOCATABLE :: IDF IF(ALLOCATED(IDF))DEALLOCATE(IDF) ALLOCATE(IDF(2)); DO I=1,SIZE(IDF); CALL IDFNULLIFY(IDF(I)); ENDDO IF(.NOT.IDFREAD(IDF(1),IDFFNAME1,1))RETURN IF(.NOT.IDFREAD(IDF(2),IDFFNAME2,1))RETURN IU=UTL_GETUNIT() CALL OSD_OPEN(IU,FILE=IPFFNAME,STATUS='OLD',ACTION='READ') READ(IU,*) NSTUW READ(IU,*) NC DO I=1,NC; READ(IU,*); END DO READ(IU,*) CALL POLYGON1DEALLOCATE_SELIDF(); ALLOCATE(SELIDF(1)); CALL IDFNULLIFY(SELIDF(1)) CALL IDFCOPY(IDF(1),SELIDF(1)) IF(.NOT.ASSOCIATED(SELIDF(1)%X))THEN ALLOCATE(SELIDF(1)%X(SELIDF(1)%NCOL,SELIDF(1)%NROW)) ENDIF SELIDF(1)%X=SELIDF(1)%NODATA MAXTHREAD=1000; MAXN=MAXTHREAD; ALLOCATE(ISPEC(MAXTHREAD),THREAD(3,MAXTHREAD),YSEL(2,MAXTHREAD)) IMENU=3 !le ! IMENU=2 !lt DTERM=0 IRAT =0 DO I=1,NSTUW READ(IU,*,IOSTAT=IOS) XC,YC,Z,ORIENT,DH !## can be less, whenever current window/segment is gridded IF(IOS.NE.0)EXIT IF(DH.NE.0.0D0)THEN CALL IDFIROWICOL(IDF(1),IROW,ICOL,XC,YC) IF(IROW.GT.0.AND.ICOL.GT.0)THEN NTHREAD=0; CALL ISG2GRIDCOMPUTESTUWENANGLE(ORIENT,NTHREAD,XC,YC,YSEL,MAXTHREAD,IDF(1)) !## only whenever current waterlevel is less that structure level IF(IDF(1)%X(ICOL,IROW).LT.Z)THEN CALL IDFEDITTRACE(IDF(1),SELIDF(1),THREAD,YSEL,ISPEC,DTERM,IMENU,MAXTHREAD,MAXN,Z,NTHREAD,1, & X2CRIT=IDF(1)%X(ICOL,IROW)) !-0.1) ENDIF DO J=1,NTHREAD ICOL=INT(YSEL(1,J)); IROW=INT(YSEL(2,J)) IF(ICOL.GT.0.AND.ICOL.LE.IDF(1)%NCOL.AND. & IROW.GT.0.AND.IROW.LE.IDF(1)%NROW)THEN IF(SELIDF(1)%X(ICOL,IROW).GT.0)THEN IDF(1)%X(ICOL,IROW)=Z IDF(2)%X(ICOL,IROW)=REAL(I) ENDIF SELIDF(1)%X(ICOL,IROW)=SELIDF(1)%NODATA ENDIF END DO CALL UTL_WAITMESSAGE(IRAT,IRAT1,I,NSTUW,'Progress searching structures ') ENDIF ENDIF ENDDO IF(.NOT.IDFWRITE(IDF(1),IDFFNAME1,1))RETURN IF(.NOT.IDFWRITE(IDF(2),IDFFNAME2,1))RETURN CALL IDFDEALLOCATE(IDF,SIZE(IDF)) CALL POLYGON1DEALLOCATE_SELIDF() CLOSE(IU,STATUS='DELETE') DEALLOCATE(ISPEC,THREAD,YSEL,IDF) END SUBROUTINE ISG2GRIDCOMPUTESTUWEN !###=============================================================================== SUBROUTINE ISG2GRIDCOMPUTE_GETXYORIENT(IISG,IIST,X,Y,ANGLE) !###=============================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),PARAMETER :: R2G=360.0D0/(2.0*3.1415) !1 rad = 57.15 degrees REAL(KIND=DP_KIND),INTENT(OUT) :: X,Y,ANGLE INTEGER,INTENT(IN) :: IISG,IIST INTEGER :: I,J REAL(KIND=DP_KIND) :: DXY,F,TD TD=0.0D0 J =ISG(IISG)%ISEG DO I=2,ISG(IISG)%NSEG J=J+1 DXY=((ISP(J)%X-ISP(J-1)%X)**2.0D0)+((ISP(J)%Y-ISP(J-1)%Y)**2.0D0) IF(DXY.GT.0.0D0)DXY=SQRT(DXY) TD=TD+DXY IF(TD.GE.IST(IIST)%DIST)EXIT END DO !## distance current segment F = (IST(IIST)%DIST-(TD-DXY))/DXY X = ISP(J-1)%X+(ISP(J)%X-ISP(J-1)%X)*F Y = ISP(J-1)%Y+(ISP(J)%Y-ISP(J-1)%Y)*F ANGLE= ATAN2(ISP(J)%Y-ISP(J-1)%Y,ISP(J)%X-ISP(J-1)%X) ANGLE= R2G*ANGLE END SUBROUTINE ISG2GRIDCOMPUTE_GETXYORIENT !###=============================================================================== REAL(KIND=DP_KIND) FUNCTION ISG2GRIDSTUWEN_GETORIENT(OR) !###=============================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),INTENT(IN) :: OR ISG2GRIDSTUWEN_GETORIENT=OR DO IF(ISG2GRIDSTUWEN_GETORIENT.LT.0.0D0) ISG2GRIDSTUWEN_GETORIENT=OR+360.0D0 IF(ISG2GRIDSTUWEN_GETORIENT.GT.360.0D0)ISG2GRIDSTUWEN_GETORIENT=OR-360.0D0 IF(ISG2GRIDSTUWEN_GETORIENT.GE.0.0D0.AND.ISG2GRIDSTUWEN_GETORIENT.LE.360.0D0)EXIT ENDDO END FUNCTION ISG2GRIDSTUWEN_GETORIENT !###====================================================================== SUBROUTINE ISG2GRIDCOMPUTESTUWENANGLE(ORIENT,NTHREAD,XC,YC,YSEL,MAXTHREAD,IDF) !###====================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),PARAMETER :: GRAD =360.0D0/(2.0*3.1415) REAL(KIND=DP_KIND),INTENT(IN) :: ORIENT,XC,YC TYPE(IDFOBJ),INTENT(IN) :: IDF INTEGER,INTENT(INOUT) :: NTHREAD INTEGER,INTENT(IN) :: MAXTHREAD INTEGER(KIND=2),DIMENSION(2,MAXTHREAD) :: YSEL REAL(KIND=DP_KIND) :: OR,X,Y,DXY INTEGER :: ICOL,IROW DXY=0.5*(IDF%DX+IDF%DY) DXY=1.0D0*(IDF%DX+IDF%DY) !## block correct side of structure OR=90.0D0 DO OR=OR+5.0; IF(OR.GT.270.0D0)EXIT X=XC-COS((ORIENT-OR)/GRAD)*DXY Y=YC-SIN((ORIENT-OR)/GRAD)*DXY CALL IDFIROWICOL(IDF,IROW,ICOL,X,Y) IF(IROW.GT.0.AND.ICOL.GT.0)THEN !## been there ... artificially IF(IDF%X(ICOL,IROW).NE.IDF%NODATA)THEN SELIDF(1)%X(ICOL,IROW)=-1.0D0; NTHREAD=NTHREAD+1; YSEL(1,NTHREAD)=INT(ICOL,2); YSEL(2,NTHREAD)=INT(IROW,2) ENDIF ENDIF END DO !## begin position CALL IDFIROWICOL(IDF,IROW,ICOL,XC,YC) NTHREAD=NTHREAD+1; YSEL(1,NTHREAD)=INT(ICOL,2); YSEL(2,NTHREAD)=INT(IROW,2); SELIDF(1)%X(ICOL,IROW)=1.0D0 END SUBROUTINE ISG2GRIDCOMPUTESTUWENANGLE !###====================================================================== SUBROUTINE ISG2GRIDPERIMETERTRAPEZIUM(XTRAP,YTRAP,WETPER,CT,BW,WDEPTH) !###====================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),DIMENSION(4),INTENT(IN) :: XTRAP,YTRAP REAL(KIND=DP_KIND),INTENT(IN) :: WDEPTH REAL(KIND=DP_KIND),INTENT(OUT) :: WETPER,BW,CT REAL(KIND=DP_KIND) :: DX,DY BW=XTRAP(3)-XTRAP(4) DX=XTRAP(2)-XTRAP(3) DY=YTRAP(2)-YTRAP(3) IF(DY.LE.0.0D0)CALL WMESSAGEBOX(OKONLY,EXCLAMATIONICON,COMMONOK,'ERROR DY.LE.0.0D0 in ISG2GRIDPERIMETERTRAPEZIUM() !!!'//CHAR(13)// & 'iMOD will be terminated, call 030-2564766','Error') !## cotanges CT=DX/DY DX=0.0D0 IF(CT.NE.0.0D0)DX=WDEPTH*CT !/(1.0D0/CT) WETPER=SQRT((DX**2.0D0)+(WDEPTH**2.0D0)) WETPER=WETPER*2.0+BW END SUBROUTINE ISG2GRIDPERIMETERTRAPEZIUM !###==================================================================== SUBROUTINE ISG2GRIDGETPARAM(DISTANCE,BOTTOM,N,BH,WP,RWIDTH,WETPER,MINDEPTH) !###==================================================================== !RE.: In case waterlevel exceeds cross-section, rwidth and wetper will ! yield max. values based upon cross-section information only! !###==================================================================== IMPLICIT NONE INTEGER,INTENT(IN) :: N REAL(KIND=DP_KIND),INTENT(IN),DIMENSION(N) :: DISTANCE,BOTTOM REAL(KIND=DP_KIND),INTENT(IN) :: BH,WP,MINDEPTH REAL(KIND=DP_KIND),INTENT(OUT) :: RWIDTH,WETPER INTEGER :: I REAL(KIND=DP_KIND) :: DISTB,DISTW,FACTOR,X1,X2,WPCOR,X WPCOR =WP-BH !## waterdepth interpolated WETPER=0.0D0 RWIDTH=0.0D0 !## riverbed is dry - make it minimal equal to mindepth WPCOR=MAX(MINDEPTH,WPCOR) !## find left-x-coordinate - downgoing direction X1=DISTANCE(N) DO I=1,N-1 IF(BOTTOM(I).GE.WPCOR.AND.BOTTOM(I+1).LE.WPCOR)THEN DISTB = BOTTOM(I)-BOTTOM(I+1) DISTW = BOTTOM(I)-WPCOR IF(DISTB.NE.0.0D0)THEN FACTOR=DISTW/DISTB X =FACTOR*(DISTANCE(I+1)-DISTANCE(I)) X =DISTANCE(I)+X X1=MIN(X1,X) !## wetted perimeter FACTOR =(X-DISTANCE(I+1))**2.0D0+(WPCOR-BOTTOM(I+1))**2.0D0 IF(FACTOR.NE.0.0D0)FACTOR=SQRT(FACTOR) WETPER=WETPER+FACTOR ENDIF ENDIF END DO !## cross-section not high enough left IF(X1.EQ.DISTANCE(N))X1=DISTANCE(1) !## find right-x-coordinate - upwards direction X2=DISTANCE(1) DO I=N-1,1,-1 IF(BOTTOM(I+1).GE.WPCOR.AND.BOTTOM(I).LE.WPCOR)THEN DISTB = BOTTOM(I+1)-BOTTOM(I) DISTW = BOTTOM(I+1)-WPCOR IF(DISTB.NE.0.0D0)THEN FACTOR=DISTW/DISTB X =FACTOR*(DISTANCE(I+1)-DISTANCE(I)) X =DISTANCE(I+1)-X X2=MAX(X2,X) !## wetted perimeter FACTOR =(X-DISTANCE(I))**2.0D0+(WPCOR-BOTTOM(I))**2.0D0 IF(FACTOR.NE.0.0D0)FACTOR=SQRT(FACTOR) WETPER=WETPER+FACTOR ENDIF ENDIF END DO !## cross-section not high enough right IF(X2.EQ.DISTANCE(1))X2=DISTANCE(N) !## proces wetted perimeter for 'in-between' sections DO I=1,N-1 IF(BOTTOM(I).LT.WPCOR.AND.BOTTOM(I+1).LT.WPCOR)THEN !## wetted perimeter FACTOR=(DISTANCE(I)-DISTANCE(I+1))**2.0D0+(BOTTOM(I)-BOTTOM(I+1))**2.0D0 IF(FACTOR.NE.0.0D0)FACTOR=SQRT(FACTOR) WETPER=WETPER+FACTOR ENDIF END DO RWIDTH=X2-X1 END SUBROUTINE ISG2GRIDGETPARAM !###==================================================================== SUBROUTINE ISG2GRIDGETCROSS(JCRS,ICRS,NCRS,ISGLEN) !###==================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),INTENT(IN) :: ISGLEN INTEGER,INTENT(IN) :: NCRS,ICRS INTEGER,INTENT(OUT) :: JCRS INTEGER :: IC1,IC2,ISEG,I REAL(KIND=DP_KIND) :: DCL,DCR,D JCRS=0 IF(NCRS.EQ.1)THEN; JCRS=1; RETURN; ENDIF IC1 =0 !cross-section left IC2 =0 !cross-section right DCL =10.0D10 !left DCR =10.0D10 !right ISEG=ICRS-1 DO I=1,NCRS ISEG=ISEG+1 !## take 1d cross-sections only IF(ISC(ISEG)%N.GT.0)THEN D =ISC(ISEG)%DIST-ISGLEN IF(D.GT.0)THEN !## lies after !## take whenever closer than next cross-section in the front IF(D.LE.DCL)THEN DCL=D IC2=I ENDIF ELSE !## lies behind !## take whenever closer than previous cross-section in the back IF(ABS(D).LE.DCR)THEN DCR=ABS(D) IC1=I ENDIF ENDIF ENDIF ENDDO !## only use cross-section in front if nothing found in the back IF(IC1.EQ.0)IC1=IC2 JCRS=IC1 END SUBROUTINE ISG2GRIDGETCROSS !###==================================================================== SUBROUTINE ISG2GRIDGETQHR(JQHR,IQHR,NQHR,ISGLEN) !###==================================================================== IMPLICIT NONE REAL(KIND=DP_KIND),INTENT(IN) :: ISGLEN INTEGER,INTENT(IN) :: NQHR,IQHR INTEGER,INTENT(OUT) :: JQHR INTEGER :: ISEG,I REAL(KIND=DP_KIND) :: D,MAXD JQHR=1 IF(NQHR.EQ.1)RETURN !#take qh-point nearest by MAXD=10.0D10 ISEG=IQHR-1 DO I=1,NQHR ISEG=ISEG+1 D =ABS(ISQ(ISEG)%DIST-ISGLEN) IF(D.LT.MAXD)THEN MAXD=D JQHR=I ENDIF END DO END SUBROUTINE ISG2GRIDGETQHR END MODULE MOD_ISG_GRID