The TELEMAC-MASCARET system  trunk
lecbreach.f
Go to the documentation of this file.
1 ! ********************
2  SUBROUTINE lecbreach
3 ! ********************
4 !
5  &(ific)
6 !
7 !***********************************************************************
8 ! TELEMAC2D V7P2
9 !***********************************************************************
10 !
11 !brief READ THE BREACHES DATA FILE, ALLOCATE THE DEDICATED ARRAY
12 !+ AND IDENTIFY THE NODES
13 !
14 !
15 !history P. CHASSE (CETMEF) / C.COULET (ARTELIA)
16 !+ 03/08/2012
17 !+ V6P2
18 !+ Creation
19 !
20 !history Y.B. TADESSE (TUHH, INSTITUTE OF RIVER AND COASTAL ENGINEERING)
21 !+ 14/02/2014
22 !+ V6P3R2
23 !+ Addition of later breach growth option
24 !
25 !history J-M HERVOUET (EDF LAB, LNHE)
26 !+ 01/07/2016
27 !+ V7P2
28 !+ ITMP now an allocatable array, not automatic array ITMP(NPOIN).
29 !
30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31 !| IFIC |-->| LOGICAL UNIT OF BREACHES DATA FILE
32 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33 !
34  USE bief
36 !
38  IMPLICIT NONE
39 !
40 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
41 !
42  INTEGER , INTENT(IN) :: IFIC
43 !
44 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
45 !
46  INTEGER N, M, NUM, NBL, ISTAT
47  DOUBLE PRECISION X1, X2, Y1, Y2, DX, DY
48  DOUBLE PRECISION U1, U2, V1, V2, DELS, MIDDIS
49  DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: XL, YL, XP, YP
50  DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: DS
51 !
52  CHARACTER(LEN=6) :: NOM, NOMX, NOMY, NOMDS
53  CHARACTER(LEN=1), PARAMETER :: CHIFFRE(0:9) =
54  & (/'0','1','2','3','4','5','6','7','8','9'/)
55 !
56  INTEGER, ALLOCATABLE :: ITMP(:)
57  CHARACTER(LEN=1) :: COMMENT = "#"
58  INTEGER :: IERR
59 !
60 !-----------------------------------------------------------------------
61 !-----------------------------------------------------------------------
62 !
63  ALLOCATE(itmp(npoin))
64 !
65  CALL skip_comment_line(ific, comment, ierr)
66  IF(ierr.NE.0) GOTO 900
67  READ(ific,*,err=999) nbrech
68 !
69 ! ALLOCATION OF SPECIFIC ARRAYS
70 !
71  IF(nbrech.GT.0) THEN
72  CALL bief_allvec(2,optnbr,'OPTNBR',nbrech,1,0,mesh)
73  CALL bief_allvec(2,optero,'OPTERO',nbrech,1,0,mesh)
74  CALL bief_allvec(1,tdecbr,'TDECBR',nbrech,1,0,mesh)
75  CALL bief_allvec(1,durbr ,'DURBR ',nbrech,1,0,mesh)
76  CALL bief_allvec(1,zfinbr,'ZFINBR',nbrech,1,0,mesh)
77  CALL bief_allvec(1,zdecbr,'ZDECBR',nbrech,1,0,mesh)
78  CALL bief_allvec(1,zcrbr ,'ZCRBR ',nbrech,1,0,mesh)
79  CALL bief_allvec(1,polwdt,'POLWDT',nbrech,1,0,mesh)
80  CALL bief_allvec(2,numpsd,'NUMPSD',nbrech,1,0,mesh)
81  CALL bief_allvec(2,nbndbr,'NBNDBR',nbrech,1,0,mesh)
82  CALL bief_allvec(2,nponbr,'NPONBR',nbrech,1,0,mesh)
83  CALL bief_allvec(1,curbrw,'CURBRW',nbrech,1,0,mesh)
84  CALL bief_allvec(1,finbrw,'FINBRW',nbrech,1,0,mesh)
85  CALL bief_allvec(1,inibrw,'INIBRW',nbrech,1,0,mesh)
86  ELSE
87  CALL bief_allvec(2,optnbr,'OPTNBR',0,1,0,mesh)
88  CALL bief_allvec(2,optero,'OPTERO',0,1,0,mesh)
89  CALL bief_allvec(1,tdecbr,'TDECBR',0,1,0,mesh)
90  CALL bief_allvec(1,durbr ,'DURBR ',0,1,0,mesh)
91  CALL bief_allvec(1,zfinbr,'ZFINBR',0,1,0,mesh)
92  CALL bief_allvec(1,zdecbr,'ZDECBR',0,1,0,mesh)
93  CALL bief_allvec(1,zcrbr ,'ZCRBR ',0,1,0,mesh)
94  CALL bief_allvec(1,polwdt,'POLWDT',0,1,0,mesh)
95  CALL bief_allvec(2,numpsd,'NUMPSD',0,1,0,mesh)
96  CALL bief_allvec(2,nbndbr,'NBNDBR',0,1,0,mesh)
97  CALL bief_allvec(2,nponbr,'NPONBR',0,1,0,mesh)
98  CALL bief_allvec(1,curbrw,'CURBRW',0,1,0,mesh)
99  CALL bief_allvec(1,finbrw,'FINBRW',0,1,0,mesh)
100  CALL bief_allvec(1,inibrw,'INIBRW',0,1,0,mesh)
101  ENDIF
102  CALL allblo(indbr ,'INDBR ')
103  CALL allblo(dkaxcr,'DKAXCR')
104  CALL allblo(dkaycr,'DKAYCR')
105  CALL allblo(pondsb,'PONDSB')
106 !
107 !
108  DO n = 1, nbrech
109  CALL skip_comment_line(ific, comment, ierr)
110  IF(ierr.NE.0) GOTO 900
111  READ(ific,*,err=998) polwdt%R(n)
112  CALL skip_comment_line(ific, comment, ierr)
113  IF(ierr.NE.0) GOTO 900
114  READ(ific,*,err=997) optnbr%I(n)
115  CALL skip_comment_line(ific, comment, ierr)
116  IF(ierr.NE.0) GOTO 900
117  IF(optnbr%I(n).EQ.1) THEN
118  READ(ific,*,err=996) tdecbr%R(n)
119  CALL skip_comment_line(ific, comment, ierr)
120  IF(ierr.NE.0) GOTO 900
121  ELSE
122  tdecbr%R(n) = -9999.d0
123  ENDIF
124  READ(ific,*,err=995) durbr%R(n)
125  CALL skip_comment_line(ific, comment, ierr)
126  IF(ierr.NE.0) GOTO 900
127  READ(ific,*,err=810) optero%I(n)
128  CALL skip_comment_line(ific, comment, ierr)
129  IF(ierr.NE.0) GOTO 900
130  READ(ific,*,err=994) zfinbr%R(n)
131  CALL skip_comment_line(ific, comment, ierr)
132  IF(ierr.NE.0) GOTO 900
133  IF(optnbr%I(n).EQ.3) THEN
134  READ(ific,*,err=993) numpsd%I(n)
135  CALL skip_comment_line(ific, comment, ierr)
136  IF(ierr.NE.0) GOTO 900
137  IF(ncsize.GT.1) THEN
138  num = numpsd%I(n)
139  numpsd%I(n) = 0
140  DO m=1,mesh%NPOIN
141  IF(num.EQ.mesh%KNOLG%I(m)) THEN
142  numpsd%I(n) = m
143  ENDIF
144  ENDDO
145  ENDIF
146  ENDIF
147  IF(optnbr%I(n).NE.1) THEN
148  READ(ific,*,err=992) zdecbr%R(n)
149  CALL skip_comment_line(ific, comment, ierr)
150  IF(ierr.NE.0) GOTO 900
151  ENDIF
152 !
153  READ(ific,*,err=991) nbl
154  CALL skip_comment_line(ific, comment, ierr)
155  IF(ierr.NE.0) GOTO 900
156 !
157 ! ALLOCATION OF LOCAL VARIABLE TO READ BREACH DEFINITION
158  istat = 0
159  ALLOCATE(xl(nbl), stat=istat)
160  CALL check_allocate(istat, "XL")
161  ALLOCATE(yl(nbl), stat=istat)
162  CALL check_allocate(istat, "YL")
163  ALLOCATE(ds(nbl-1), stat=istat)
164  CALL check_allocate(istat, "DS")
165 !
166  DO m = 1, nbl
167  READ(ific,*,err=990) xl(m), yl(m)
168  ENDDO
169 ! SEARCH MESH POINTS INSIDE THE BREACH DOMAIN
170  istat = 0
171  ALLOCATE(xp(2*nbl), stat=istat)
172  CALL check_allocate(istat, "XP")
173  ALLOCATE(yp(2*nbl), stat=istat)
174  CALL check_allocate(istat, "YP")
175 !
176  x1 = xl(1)
177  y1 = yl(1)
178  x2 = xl(2)
179  y2 = yl(2)
180  dx = x2 - x1
181  dy = y2 - y1
182  dels=sqrt(dx*dx+dy*dy)
183  IF(dels.GT.0.d0) THEN
184  u1 = dx/dels
185  u2 = dy/dels
186  ELSE
187  WRITE(lu,*) 'PROBLEM IN DEFINITION OF BREACH :',n
188  CALL plante(1)
189  ENDIF
190  v1 = -u2
191  v2 = u1
192  xp(1) = x1 + v1*polwdt%R(n)/2.d0
193  yp(1) = y1 + v2*polwdt%R(n)/2.d0
194  xp(2*nbl) = x1 - v1*polwdt%R(n)/2.d0
195  yp(2*nbl) = y1 - v2*polwdt%R(n)/2.d0
196 !
197  DO m = 2,nbl
198  x2 = xl(m)
199  y2 = yl(m)
200  dx = x2 - x1
201  dy = y2 - y1
202  ds(m-1)=sqrt(dx*dx+dy*dy)
203  IF(ds(m-1).GT.0.d0) THEN
204  u1 = dx/ds(m-1)
205  u2 = dy/ds(m-1)
206  ELSE
207  WRITE(lu,*) 'PROBLEM IN DEFINITION OF BREACH :',n
208  CALL plante(1)
209  ENDIF
210  v1 = -u2
211  v2 = u1
212  xp(m) = x2 + v1*polwdt%R(n)/2.d0
213  yp(m) = y2 + v2*polwdt%R(n)/2.d0
214  xp(2*nbl-m+1) = x2 - v1*polwdt%R(n)/2.d0
215  yp(2*nbl-m+1) = y2 - v2*polwdt%R(n)/2.d0
216  x1=x2
217  y1=y2
218  ENDDO
219 !
220  nbndbr%I(n) = 0
221  DO m = 1, npoin
222  IF(inpoly(mesh%X%R(m), mesh%Y%R(m), xp, yp, 2*nbl)) THEN
223  nbndbr%I(n) = nbndbr%I(n)+1
224  itmp(nbndbr%I(n)) = m
225  ENDIF
226  ENDDO
227 !
228  IF(n.LE.indbr%MAXBLOCK) THEN
229  nom='NBR '
230  IF(n.LT.10) THEN
231  nom(4:4) = chiffre(n)
232  ELSEIF(n.LT.100) THEN
233  nom(4:4) = chiffre(n/10)
234  nom(5:5) = chiffre(n-10*(n/10))
235  ELSEIF(n.LT.1000) THEN
236  nom(4:4) = chiffre(n/100)
237  nom(5:5) = chiffre((n-100*(n/100))/10)
238  nom(6:6) = chiffre((n-100*(n/100))-10*((n-100*(n/100))/10))
239  ELSE
240  WRITE(lu,*) 'MORE THAN 999 BREACHS ASKED IN LECBREACH'
241  CALL plante(1)
242  stop
243  ENDIF
244  CALL first_all_biefobj(indbr%ADR(n)%P)
245  CALL bief_allvec(2,indbr%ADR(n)%P,nom,nbndbr%I(n),1,0,mesh)
246  indbr%ADR(n)%P%FATHER = indbr%NAME
247  ELSE
248  WRITE(lu,*) 'LECBREACH:'
249  WRITE(lu,*) 'MORE THAN ',indbr%MAXBLOCK,'(',n,')'
250  WRITE(lu,*) 'VECTORS TO BE ALLOCATED'
251  WRITE(lu,*) 'CHANGE MAXBLOCK IN ALLBLO.'
252  CALL plante(1)
253  stop
254  ENDIF
255 !
256  DO m=1, nbndbr%I(n)
257  indbr%ADR(n)%P%I(m) = itmp(m)
258  ENDDO
259 
260  IF (optero%I(n).EQ.2) THEN
261  IF(n.LE.dkaxcr%MAXBLOCK) THEN
262  nomx='XCB '
263  nomy='YCB '
264  nomds='PDS '
265  IF(n.LT.10) THEN
266  nomx(4:4) = chiffre(n)
267  nomy(4:4) = chiffre(n)
268  nomds(4:4) = chiffre(n)
269  ELSEIF(n.LT.100) THEN
270  nomx(4:4) = chiffre(n/10)
271  nomx(5:5) = chiffre(n-10*(n/10))
272  nomy(4:4) = chiffre(n/10)
273  nomy(5:5) = chiffre(n-10*(n/10))
274  nomds(4:4) = chiffre(n/10)
275  nomds(5:5) = chiffre(n-10*(n/10))
276  ELSEIF(n.LT.1000) THEN
277  nomx(4:4) = chiffre(n/100)
278  nomx(5:5) = chiffre((n-100*(n/100))/10)
279  nomx(6:6) =
280  & chiffre((n-100*(n/100))-10*((n-100*(n/100))/10))
281  nomy(4:4) = chiffre(n/100)
282  nomy(5:5) = chiffre((n-100*(n/100))/10)
283  nomy(6:6) =
284  & chiffre((n-100*(n/100))-10*((n-100*(n/100))/10))
285  nomds(4:4) = chiffre(n/100)
286  nomds(5:5) = chiffre((n-100*(n/100))/10)
287  nomds(6:6) =
288  & chiffre((n-100*(n/100))-10*((n-100*(n/100))/10))
289  ELSE
290  WRITE(lu,*) 'MORE THAN 999 BREACHS ASKED IN LECBREACH'
291  CALL plante(1)
292  stop
293  ENDIF
294  CALL first_all_biefobj(dkaxcr%ADR(n)%P)
295  CALL first_all_biefobj(dkaycr%ADR(n)%P)
296  CALL first_all_biefobj(pondsb%ADR(n)%P)
297  CALL bief_allvec(1,dkaxcr%ADR(n)%P,nomx,nbl,1,0,mesh)
298  CALL bief_allvec(1,dkaycr%ADR(n)%P,nomy,nbl,1,0,mesh)
299  CALL bief_allvec(1,pondsb%ADR(n)%P,nomds,nbl-1,1,0,mesh)
300  dkaxcr%ADR(n)%P%FATHER = dkaxcr%NAME
301  dkaycr%ADR(n)%P%FATHER = dkaycr%NAME
302  pondsb%ADR(n)%P%FATHER = pondsb%NAME
303  ELSE
304  WRITE(lu,*) 'LECBREACH:'
305  WRITE(lu,*) 'MORE THAN ',dkaxcr%MAXBLOCK,'(',n,')'
306  WRITE(lu,*) 'VECTORS TO BE ALLOCATED'
307  WRITE(lu,*) 'CHANGE MAXBLOCK IN ALLBLO.'
308  CALL plante(1)
309  stop
310  ENDIF
311  finbrw%R(n)=0.d0
312  DO m=1,nbl
313  dkaxcr%ADR(n)%P%R(m)=xl(m)
314  dkaycr%ADR(n)%P%R(m)=yl(m)
315  IF (m .GT.1) THEN
316  pondsb%ADR(n)%P%R(m-1)=ds(m-1)
317  finbrw%R(n)=finbrw%R(n)+ds(m-1)
318  ENDIF
319  ENDDO
320  nponbr%I(n)=nbl
321  middis=0.d0
322  curbrw%R(n)=0.d0
323  inibrw%R(n)=finbrw%R(n)/10.d0
324  findmidnode: DO m=2,nbl
325  middis=middis+ds(m-1)
326  IF (middis .GE. (finbrw%R(n)/2.d0)) THEN
327  IF (m .GT. 2) THEN
328  inibrw%R(n)=max(ds(m-1)+ds(m),inibrw%R(n))
329  curbrw%R(n)=inibrw%R(n)
330  ELSE
331  WRITE(lu,*) 'INCREASE NUMBER OF POINTS OF BREACH :',n
332  GO TO 2000
333  ENDIF
334  EXIT findmidnode
335  ENDIF
336  ENDDO findmidnode
337  END IF
338  IF ((optero%I(n).NE.1) .AND. (optero%I(n).NE.2)) THEN
339  GO TO 800
340  END IF
341  DEALLOCATE(xl)
342  DEALLOCATE(yl)
343  DEALLOCATE(xp)
344  DEALLOCATE(yp)
345  DEALLOCATE(ds)
346 !
347  ENDDO
348 !
349  indbr%N = nbrech
350  dkaxcr%N = nbrech
351  dkaycr%N = nbrech
352  pondsb%N = nbrech
353  GOTO 1000
354 !
355 !-----------------------------------------------------------------------
356 ! ERROR MESSAGES
357 !-----------------------------------------------------------------------
358 !
359 810 CONTINUE
360  WRITE(lu,*) 'LECBREACH : READ ERROR ON THE'
361  WRITE(lu,*) ' BREACHES DATA FILE'
362  WRITE(lu,*) ' FOR THE BREACH ',n
363  WRITE(lu,*) ' EROSION OPTION CANNOT BE READ'
364  GO TO 2000
365 !
366 800 CONTINUE
367  WRITE(lu,*) 'LECBREACH: UNAVAILABLE EROSION OPTION'
368  GO TO 2000
369 !
370 999 CONTINUE
371  WRITE(lu,*) 'LECBREACH: READ ERROR ON THE BREACHES DATA FILE'
372  WRITE(lu,*) ' AT LINE 2'
373  GO TO 2000
374 !
375 998 CONTINUE
376  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
377  WRITE(lu,*) ' BREACHES DATA FILE'
378  WRITE(lu,*) ' FOR THE BREACH ',n
379  WRITE(lu,*) ' BREACH POLYGON WIDTH CANNOT BE READ'
380  GO TO 2000
381 !
382 997 CONTINUE
383  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
384  WRITE(lu,*) ' BREACHES DATA FILE'
385  WRITE(lu,*) ' FOR THE BREACH ',n
386  WRITE(lu,*) ' OPTION CANNOT BE READ'
387  GO TO 2000
388 !
389 996 CONTINUE
390  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
391  WRITE(lu,*) ' BREACHES DATA FILE'
392  WRITE(lu,*) ' FOR THE BREACH ',n
393  WRITE(lu,*) ' THE STARTING TIME CANNOT BE READ'
394  GO TO 2000
395 !
396 995 CONTINUE
397  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
398  WRITE(lu,*) ' BREACHES DATA FILE'
399  WRITE(lu,*) ' FOR THE BREACH ',n
400  WRITE(lu,*) ' THE OPENNING DURATION CANNOT BE READ'
401  GO TO 2000
402 !
403 994 CONTINUE
404  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
405  WRITE(lu,*) ' BREACHES DATA FILE'
406  WRITE(lu,*) ' FOR THE BREACH ',n
407  WRITE(lu,*) ' THE FINAL LEVEL CANNOT BE READ'
408  GO TO 2000
409 !
410 993 CONTINUE
411  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
412  WRITE(lu,*) ' BREACHES DATA FILE'
413  WRITE(lu,*) ' FOR THE BREACH ',n
414  WRITE(lu,*) ' THE NUMBER OF TEST POINT CANNOT BE READ'
415  GO TO 2000
416 !
417 992 CONTINUE
418  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
419  WRITE(lu,*) ' BREACHES DATA FILE'
420  WRITE(lu,*) ' FOR THE BREACH ',n
421  WRITE(lu,*) ' THE STARTING LEVEL CANNOT BE READ'
422  GO TO 2000
423 !
424 991 CONTINUE
425  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
426  WRITE(lu,*) ' BREACHES DATA FILE'
427  WRITE(lu,*) ' FOR THE BREACH ',n
428  WRITE(lu,*) ' THE POINT NUMBER OF LINE CANNOT BE READ'
429  GO TO 2000
430 !
431 990 CONTINUE
432  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
433  WRITE(lu,*) ' BREACHES DATA FILE'
434  WRITE(lu,*) ' FOR THE BREACH ',n
435  WRITE(lu,*) ' THE COORDINATE OF POINT ',m
436  WRITE(lu,*) ' CANNOT BE READ'
437  GO TO 2000
438 !
439 900 CONTINUE
440  WRITE(lu,*) 'BRECHE : READ ERROR ON THE'
441  WRITE(lu,*) ' BREACHES DATA FILE'
442  WRITE(lu,*) ' UNEXPECTED END OF FILE'
443 !
444 2000 CONTINUE
445 !
446  CALL plante(1)
447  stop
448 !
449 1000 CONTINUE
450 !
451 !-----------------------------------------------------------------------
452 !
453  DEALLOCATE(itmp)
454 !
455 !-----------------------------------------------------------------------
456 !
457  RETURN
458  END
subroutine allblo(BLO, NOM)
Definition: allblo.f:7
type(bief_obj), target dkaxcr
subroutine bief_allvec(NAT, VEC, NOM, IELM, DIM2, STATUT, MESH)
Definition: bief_allvec.f:7
type(bief_obj), target zfinbr
type(bief_obj), target finbrw
type(bief_obj), target pondsb
logical function inpoly(X, Y, XSOM, YSOM, NSOM)
Definition: inpoly.f:7
type(bief_obj), target polwdt
type(bief_obj), target nponbr
type(bief_obj), target zdecbr
type(bief_obj), target durbr
type(bief_obj), target optnbr
type(bief_obj), target inibrw
type(bief_mesh), target mesh
type(bief_obj), target optero
subroutine lecbreach(IFIC)
Definition: lecbreach.f:7
type(bief_obj), target dkaycr
type(bief_obj), target curbrw
subroutine first_all_biefobj(OBJ)
type(bief_obj), target zcrbr
type(bief_obj), target nbndbr
type(bief_obj), target tdecbr
type(bief_obj), target indbr
type(bief_obj), target numpsd
Definition: bief.f:3