The TELEMAC-MASCARET system  trunk
breach.f
Go to the documentation of this file.
1 ! ******************
2  SUBROUTINE breach
3 ! ******************
4 !
5 !
6 !***********************************************************************
7 ! TELEMAC2D V6P2 03/08/2012
8 !***********************************************************************
9 !
10 !brief MODIFICATION OF THE BOTTOM TOPOGRAPHY FOR BREACHES
11 !
12 !
13 !history P. CHASSE (CETMEF) / C. COULET (ARTELIA)
14 !+ 03/08/2012
15 !+ V6P2
16 !+ Creation
17 !
18 !history Y.B. TADESSE (TUHH, INSTITUTE OF RIVER AND COASTAL ENGINEERING)
19 !+ 14/02/2014
20 !+ V6P3R2
21 !+ Addition of later breach growth option
22 !
23 !history J,RIEHME (ADJOINTWARE)
24 !+ November 2016
25 !+ V7P2
26 !+ Replaced EXTERNAL statements to parallel functions / subroutines
27 !+ by the INTERFACE_PARALLEL
28 !
29 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31 !
32  USE bief
34 !
36  USE interface_parallel, ONLY : p_sum,p_max,p_min,
37  & p_sum
38  IMPLICIT NONE
39 !
40 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
41 !
42 !
43  INTEGER I, J, K, N, M,END1, END2,CURNBR
44  INTEGER ISTAT,VECZ
45  INTEGER TEMPND(npoin)
46  DOUBLE PRECISION ZC, ZW, ZB
47  DOUBLE PRECISION AT1, AT2, AT3
48 !
49  DOUBLE PRECISION X1, X2, Y1, Y2, DX, DY,DS1,DS2,END1X,END1Y
50  DOUBLE PRECISION U1, U2, V1, V2, DELS,CURDIS,DIS,END2X,END2Y
51  DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: XL, YL, XP, YP
52 !
53 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
54 !
55  IF (.NOT.dejalu_breach) THEN
56  IF(debug.GT.0) WRITE(lu,*) 'CALLING LECBREACH'
57  CALL lecbreach(t2d_files(t2dbrc)%LU)
58  IF(debug.GT.0) WRITE(lu,*) 'BACK FROM LECBREACH'
59  dejalu_breach=.true.
60  WRITE (lu,*) 'READING BREACH DATA = OK'
61 !
62  DO i = 1, nbrech
63  zcrbr%R(i) = -huge(100.d0)
64  DO j = 1, nbndbr%I(i)
65  k = indbr%ADR(i)%P%I(j)
66  IF(zf%R(k).GT.zcrbr%R(i)) THEN
67  zcrbr%R(i) = zf%R(k)
68  ENDIF
69  ENDDO
70  IF(ncsize.GT.1) THEN
71  zcrbr%R(i) = p_max(zcrbr%R(i))
72  ENDIF
73  ENDDO
74  ENDIF
75 !
76  DO i = 1, nbrech
77  IF (optero%I(i).EQ.1) THEN
78  IF((optnbr%I(i).EQ.2).AND.(tdecbr%R(i).LT.0.d0)) THEN
79  zc = 0.d0
80  n = 0
81  DO j = 1, nbndbr%I(i)
82  k = indbr%ADR(i)%P%I(j)
83  IF(h%R(k).GT.0.d0) THEN
84  zw = zf%R(k)+h%R(k)
85  zc = zc + zw
86  n = n + 1
87  ENDIF
88  ENDDO
89  IF(ncsize.GT.1) THEN
90  n = p_sum(n)
91  zc = p_sum(zc)
92  ENDIF
93  IF(n.GT.1) zc = zc/n
94  IF(zc.GT.zdecbr%R(i)) THEN
95  WRITE(lu,20) i, at
96  tdecbr%R(i) = at
97  ENDIF
98  ENDIF
99  IF((optnbr%I(i).EQ.3).AND.(tdecbr%R(i).LT.0.d0)) THEN
100  IF(numpsd%I(i).GT.0) THEN
101  zw = zf%R(numpsd%I(i)) + h%R(numpsd%I(i))
102  ELSE
103  zw = 0.d0
104  ENDIF
105 ! CASE WHERE ONE OF THE ENDS IS NOT IN THE SUB-DOMAIN
106  IF(ncsize.GT.1) THEN
107  zw = p_max(zw)+p_min(zw)
108  ENDIF
109  IF(zw.GT.zdecbr%R(i)) THEN
110  WRITE(lu,20) i, at
111  tdecbr%R(i) = at
112  ENDIF
113  ENDIF
114  at1 = tdecbr%R(i)
115  at2 = at1 + durbr%R(i)
116  IF(at1.GT.0.d0) THEN
117  IF(at.GT.at1) THEN
118  IF(at.GT.at2) THEN
119  zb = zfinbr%R(i)
120  ELSE
121  zb = zcrbr%R(i)+(zfinbr%R(i)-zcrbr%R(i))/(at2-at1)
122  & *(at-at1)
123  ENDIF
124  DO j = 1, nbndbr%I(i)
125  k = indbr%ADR(i)%P%I(j)
126  zf%R(k)=min(zf%R(k), zb)
127  ENDDO
128  ENDIF
129  ENDIF
130  ELSEIF (optero%I(i).EQ.2) THEN
131 ! FIND NUMBER POINTS IN THE CURRENT BREACH REGION
132  curdis=0.d0
133  END1X=0.d0
134  END1Y=0.d0
135  END2X=0.d0
136  END2Y=0.d0
137  END1=0
138  END2=0
139  firstend: DO m=1,nponbr%I(i)-1
140  curdis=curdis+pondsb%ADR(i)%P%R(m)
141  IF (curdis .GT. ((finbrw%R(i)-curbrw%R(i))/2.d0)) THEN
142  END1=M+1
143  ds1= curdis - (finbrw%R(i)-curbrw%R(i))/2.d0
144  dx = dkaxcr%ADR(i)%P%R(m+1) - dkaxcr%ADR(i)%P%R(m)
145  dy = dkaycr%ADR(i)%P%R(m+1) - dkaycr%ADR(i)%P%R(m)
146  u1 = dx/pondsb%ADR(i)%P%R(m)
147  u2 = dy/pondsb%ADR(i)%P%R(m)
148  IF(curbrw%R(i).LT.finbrw%R(i)) THEN
149  END1X = DKAXCR%ADR(I)%P%R(M+1) - U1*ds1
150  END1Y = DKAYCR%ADR(I)%P%R(M+1) - U2*ds1
151  ELSEIF(curbrw%R(i).EQ.finbrw%R(i)) THEN
152  END1X = DKAXCR%ADR(I)%P%R(1)- U1*
153  & (pondsb%ADR(i)%P%R(1)/10.d0)
154  END1Y = DKAYCR%ADR(I)%P%R(1)- U2*
155  & (pondsb%ADR(i)%P%R(1)/10.d0)
156  ELSE
157  WRITE(lu,400)
158  stop
159  ENDIF
160  dis=ds1
161  secondend: DO j=m+1,nponbr%I(i)-1
162  dis=dis+pondsb%ADR(i)%P%R(j)
163  IF (dis .GE. curbrw%R(i)) THEN
164  END2=j
165  ds2=dis - curbrw%R(i)
166  dx = dkaxcr%ADR(i)%P%R(j+1) - dkaxcr%ADR(i)%P%R(j)
167  dy = dkaycr%ADR(i)%P%R(j+1) - dkaycr%ADR(i)%P%R(j)
168  u1 = dx/pondsb%ADR(i)%P%R(j)
169  u2 = dy/pondsb%ADR(i)%P%R(j)
170  IF(curbrw%R(i).LT.finbrw%R(i)) THEN
171  END2X = DKAXCR%ADR(I)%P%R(J+1) - U1*ds2
172  END2Y = DKAYCR%ADR(I)%P%R(J+1) - U2*ds2
173  ELSEIF (curbrw%R(i).EQ.finbrw%R(i)) THEN
174  END2X = DKAXCR%ADR(I)%P%R(NPONBR%I(I))+ U1*
175  & (pondsb%ADR(i)%P%R(nponbr%I(i)-1)/10.d0)
176  END2Y = DKAYCR%ADR(I)%P%R(NPONBR%I(I))+ U2*
177  & (pondsb%ADR(i)%P%R(nponbr%I(i)-1)/10.d0)
178  ELSE
179  WRITE(lu,400) i
180  stop
181  ENDIF
182  EXIT secondend
183  ENDIF
184  END DO secondend
185 !
186  EXIT firstend
187  ENDIF
188  ENDDO firstend
189 400 FORMAT(1x, 'BREACH: BREACH WIDTH CANNOT BE GREATER ',
190  & 'THAN FINAL WIDTH FOR BREACH: ', 1i6)
191 ! ALLOCATION OF LOCAL VARIABLE FOR BREACH DEFINITION
192  istat = 0
193  vecz=end2-end1+3
194  ALLOCATE(xl(vecz), stat=istat)
195  IF(istat.NE.0) THEN
196  WRITE(lu,200) istat
197  stop
198  ENDIF
199  ALLOCATE(yl(vecz), stat=istat)
200  IF(istat.NE.0) THEN
201  WRITE(lu,200) istat
202  stop
203  ENDIF
204 !
205 200 FORMAT(1x,'ERROR DURING ALLOCATION OF VECTOR: ',
206  & 'ERROR CODE: ',1i6)
207 !
208  xl(1)=end1x
209  yl(1)=end1y
210  xl(vecz)=end2x
211  yl(vecz)=end2y
212  DO j=2,vecz-1
213  xl(j)=dkaxcr%ADR(i)%P%R(end1+j-2)
214  yl(j)=dkaycr%ADR(i)%P%R(end1+j-2)
215  ENDDO
216 ! SEARCH MESH POINTS INSIDE THE BREACH DOMAIN
217  istat = 0
218  ALLOCATE(xp(2*vecz), stat=istat)
219  IF(istat.NE.0) THEN
220  WRITE(lu,200) istat
221  stop
222  ENDIF
223  ALLOCATE(yp(2*vecz), stat=istat)
224  IF(istat.NE.0) THEN
225  WRITE(lu,200) istat
226  stop
227  ENDIF
228 !
229  x1 = xl(1)
230  y1 = yl(1)
231  x2 = xl(2)
232  y2 = yl(2)
233  dx = x2 - x1
234  dy = y2 - y1
235  dels=sqrt(dx*dx+dy*dy)
236  IF(dels.GE.0.d0) THEN
237  u1 = dx/dels
238  u2 = dy/dels
239  ELSE
240  WRITE(lu,*) 'PROBLEM IN DEFINITION OF BREACH :',i
241  CALL plante(1)
242  stop
243  ENDIF
244  v1 = -u2
245  v2 = u1
246  xp(1) = x1 + v1*polwdt%R(i)/2.d0
247  yp(1) = y1 + v2*polwdt%R(i)/2.d0
248  xp(2*vecz) = x1 - v1*polwdt%R(i)/2.d0
249  yp(2*vecz) = y1 - v2*polwdt%R(i)/2.d0
250 !
251  DO m = 2,vecz
252  x2 = xl(m)
253  y2 = yl(m)
254  dx = x2 - x1
255  dy = y2 - y1
256  dels=sqrt(dx*dx+dy*dy)
257  IF(dels.GE.0.d0) THEN
258  u1 = dx/dels
259  u2 = dy/dels
260  ELSE
261  WRITE(lu,*) 'PROBLEM IN DEFINITION OF BREACH :',i
262  CALL plante(1)
263  stop
264  ENDIF
265  v1 = -u2
266  v2 = u1
267  xp(m) = x2 + v1*polwdt%R(i)/2.d0
268  yp(m) = y2 + v2*polwdt%R(i)/2.d0
269  xp(2*vecz-m+1) = x2 - v1*polwdt%R(i)/2.d0
270  yp(2*vecz-m+1) = y2 - v2*polwdt%R(i)/2.d0
271  x1=x2
272  y1=y2
273  ENDDO
274  curnbr = 0
275  DO m = 1, npoin
276  IF(inpoly(mesh%X%R(m), mesh%Y%R(m), xp, yp, 2*vecz)) THEN
277  curnbr = curnbr +1
278  tempnd(curnbr) = m
279  ENDIF
280  ENDDO
281 !
282  IF((optnbr%I(i).EQ.2).AND.(tdecbr%R(i).LT.0.d0)) THEN
283  zc = 0.d0
284  n = 0
285  DO j = 1, curnbr
286  k = tempnd(j)
287  IF(h%R(k).GT.0.d0) THEN
288  zw = zf%R(k)+h%R(k)
289  zc = zc + zw
290  n = n + 1
291  ENDIF
292  ENDDO
293  IF(ncsize.GT.1) THEN
294  n = p_sum(n)
295  zc = p_sum(zc)
296  ENDIF
297  IF(n.GT.1) zc = zc/n
298  IF(zc.GT.zdecbr%R(i)) THEN
299  WRITE(lu,20) i, at
300  tdecbr%R(i) = at
301  ENDIF
302  ENDIF
303  IF((optnbr%I(i).EQ.3).AND.(tdecbr%R(i).LT.0.d0)) THEN
304  IF(numpsd%I(i).GT.0) THEN
305  zw = zf%R(numpsd%I(i)) + h%R(numpsd%I(i))
306  ELSE
307  zw = 0.d0
308  ENDIF
309 ! CASE WHERE ONE OF THE ENDS IS NOT IN THE SUB-DOMAIN
310  IF(ncsize.GT.1) THEN
311  zw = p_max(zw)+p_min(zw)
312  ENDIF
313  IF(zw.GT.zdecbr%R(i)) THEN
314  WRITE(lu,20) i, at
315  tdecbr%R(i) = at
316  ENDIF
317  ENDIF
318 !
319 ! DURATION FOR THE VERTICAL EROSION IS TAKEN AS A TENTH
320 ! OF THE TOTAL BREACH DURATION
321  at1 = tdecbr%R(i)
322  at2 = at1 + durbr%R(i)
323  at3= at1 + durbr%R(i)/10.d0
324  IF(at1.GT.0.d0) THEN
325  IF(at.GT.at1) THEN
326  IF(at.GT.at3) THEN
327  zb = zfinbr%R(i)
328  ELSE
329  zb = zcrbr%R(i)+(zfinbr%R(i)-zcrbr%R(i))/(at3-at1)
330  & *(at-at1)
331  ENDIF
332  IF(at.GT.at2) THEN
333  curbrw%R(i)=finbrw%R(i)
334  ELSE
335  curbrw%R(i)=inibrw%R(i)+(finbrw%R(i)-inibrw%R(i))
336  & /(at2-at1)*(at-at1)
337  ENDIF
338  DO j = 1, curnbr
339  k = tempnd(j)
340  zf%R(k)=min(zf%R(k), zb)
341  ENDDO
342  ENDIF
343  ENDIF
344  DEALLOCATE(xl)
345  DEALLOCATE(yl)
346  DEALLOCATE(xp)
347  DEALLOCATE(yp)
348  ENDIF
349  ENDDO
350 !
351 !-----------------------------------------------------------------------
352 ! MESSAGES
353 20 FORMAT(1x,'CREATION OF BREACH : ',i4,/,1x,
354  & 'AT TIME : ',g16.7)
355 !-----------------------------------------------------------------------
356 !
357  RETURN
358  END
359 
type(bief_obj), target h
type(bief_obj), target zfinbr
logical function inpoly(X, Y, XSOM, YSOM, NSOM)
Definition: inpoly.f:7
subroutine breach
Definition: breach.f:4
type(bief_obj), target zdecbr
type(bief_obj), target durbr
type(bief_obj), target optnbr
type(bief_obj), target zf
type(bief_obj), target optero
subroutine lecbreach(IFIC)
Definition: lecbreach.f:7
type(bief_obj), target zcrbr
double precision, target at
type(bief_obj), target nbndbr
type(bief_obj), target tdecbr
type(bief_obj), target indbr
type(bief_obj), target numpsd
type(bief_file), dimension(maxlu_t2d), target t2d_files
Definition: bief.f:3