The TELEMAC-MASCARET system  trunk
flusec_sisyphe.f
Go to the documentation of this file.
1 ! *************************
2  SUBROUTINE flusec_sisyphe
3 ! *************************
4 !
5  &(u,v,h,qsxc,qsyc,charr,qsxs,qsys,susp,
6  & ikle,nelmax,nelem,x,y,dt,ncp,ctrlsc,info,tps)
7 !
8 !***********************************************************************
9 ! SISYPHE V7P0 21/07/2011
10 !***********************************************************************
11 !
12 !brief COMPUTES FLUXES THROUGH CONTROL SECTIONS
13 !+ AND ADDS THEM UP TO OBTAIN OSCILLATING VOLUMES.
14 !+
15 !+ MESHES OF DIMENSION 2 AND CONSIDERED WATER DEPTH.
16 !
17 !history J-M HERVOUET (LNHE)
18 !+ 27/12/2006
19 !+ V5P7
20 !+
21 !
22 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
23 !+ 13/07/2010
24 !+ V6P0
25 !+ Translation of French comments within the FORTRAN sources into
26 !+ English comments
27 !
28 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
29 !+ 21/08/2010
30 !+ V6P0
31 !+ Creation of DOXYGEN tags for automated documentation and
32 !+ cross-referencing of the FORTRAN sources
33 !
34 !history J-M HERVOUET (EDF R&D, LNHE)
35 !+ 02/01/2014
36 !+ V7P0
37 !+ KNOGL removed.
38 !
39 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
40 !| CHARR |-->| LOGICAL, BEDLOAD
41 !| CTRLSC |-->| CONTROL SECTION
42 !| DT |-->| TIME STEP
43 !| H |-->| WATER DEPTH
44 !| IKLE |-->| CONNECTIVITY TABLE
45 !| INFO |-->| IF YES, PRINT
46 !| NCP |-->| TWO TIMES THE NUMBER OF CONTROL SECTIONS
47 !| NELEM |-->| NUMBER OF ELEMENTS
48 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
49 !| QSXC |<->| BEDLOAD TRANSPORT RATE X-DIRECTION
50 !| QSXS |<->| SUSPENSION TRANSPORT RATE X-DIRECTION
51 !| QSYC |<->| BEDLOAD TRANSPORT RATE Y-DIRECTION
52 !| QSYS |<->| SUSPENSION TRANSPORT RATE Y-DIRECTION
53 !| SUSP |-->| LOGICAL, SUSPENSION
54 !| TPS |-->| TEMPS
55 !| U,V |-->| VELOCITY FIELD COMPONENTS
56 !| X,Y |-->| NODES COORDINATES
57 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
58 !
59  USE bief
60  USE declarations_sisyphe, ONLY: chain,mesh,old_method_flusec,
61  & liste,deja_flusec,nseg,volneg,
62  & volpos,flx,flxs,volnegs,volposs,
63  & flxc,volnegc,volposc
64  USE interface_sisyphe, ex_flusec_sisyphe => flusec_sisyphe
65 !
66  USE declarations_special
67  IMPLICIT NONE
68 !
69 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
70 !
71  INTEGER, INTENT(IN) :: NELMAX,NELEM,NCP
72  INTEGER, INTENT(IN) :: IKLE(nelmax,*)
73  INTEGER, INTENT(IN) :: CTRLSC(ncp)
74  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),TPS,DT
75  LOGICAL, INTENT(IN) :: INFO,SUSP,CHARR
76  TYPE(bief_obj), INTENT(IN) :: U,V,H,QSXC,QSYC,QSXS,QSYS
77 !
78 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
79 !
80  INTEGER ERR
81  INTEGER, PARAMETER :: NSEMAX = 100
82 !
83  INTEGER IELEM,I1,I2,I3,ELBEST,IGBEST,ILBEST
84  INTEGER ILPREC,ISEG,ISEC,NSEC,PT,DEP,ARR
85 !
86  DOUBLE PRECISION DIST,DIST1,DIST2,DIST3
87  DOUBLE PRECISION H1,H2,X1,Y1,X2,Y2,UN1,UN2,NX,NY,SUR6
88 !
89 !
90 !-----------------------------------------------------------------------
91 !
92  sur6 = 1.d0/6.d0
93  nsec = ncp/2
94 !
95 ! LOOKS FOR WAYS THAT JOIN THE POINT COUPLES:
96 !
97  IF(.NOT.deja_flusec) THEN
98 !
99 ! DYNAMICALLY ALLOCATES FLX, VOLNEG, VOLPOS, ETC.
100 !
101  ALLOCATE(flx(nsec) ,stat=err)
102  ALLOCATE(volneg(nsec) ,stat=err)
103  ALLOCATE(volpos(nsec) ,stat=err)
104  ALLOCATE(nseg(ncp) ,stat=err)
105  ALLOCATE(liste(ncp,nsemax,2) ,stat=err)
106 ! S FOR SUSPENSION, C FOR BEDLOAD
107  ALLOCATE(flxs(nsec) ,stat=err)
108  ALLOCATE(volnegs(nsec) ,stat=err)
109  ALLOCATE(volposs(nsec) ,stat=err)
110  ALLOCATE(flxc(nsec) ,stat=err)
111  ALLOCATE(volnegc(nsec) ,stat=err)
112  ALLOCATE(volposc(nsec) ,stat=err)
113 !
114  IF(err.NE.0) THEN
115  WRITE(lu,200) err
116 200 FORMAT(1x,'FLUSEC: ERROR DURING ALLOCATION OF MEMORY: ',/,1x,
117  & 'ERROR CODE: ',1i6)
118  ENDIF
119 !
120  IF (.NOT.ALLOCATED(chain)) old_method_flusec=.true.
121 !
122  DO isec=1,nsec
123  flx(isec)=0.0d0
124  volneg(isec) =0.d0
125  volpos(isec) =0.d0
126  volnegs(isec)=0.d0
127  volposs(isec)=0.d0
128  volnegc(isec)=0.d0
129  volposc(isec)=0.d0
130  ENDDO
131 !
132  DO isec =1,nsec
133 !
134 !!JAJ #### IN THE SERIAL CASE, OR "CLASSICAL" IN PARALLEL,
135 ! FOLLOW THE ALGORITHM OF FINDING SEGMENT CHAINS
136 !
137 ! NOTE: IF YOU CHANGE THE ALGORITHM, CHANGE IT IN PARTEL AS WELL
138 !
139  IF (ncsize.LE.1 .OR. old_method_flusec) THEN
140 !
141  dep = ctrlsc(1+2*(isec-1))
142  arr = ctrlsc(2+2*(isec-1))
143  IF(ncsize.GT.1) THEN
144  dep=global_to_local_point(dep,mesh)
145  arr=global_to_local_point(arr,mesh)
146  IF(dep.EQ.0.AND.arr.EQ.0) THEN
147  nseg(isec)=0
148  EXIT
149  ENDIF
150  IF((dep.EQ.0.AND.arr.NE.0).OR.(dep.NE.0.AND.arr.EQ.0)) THEN
151  nseg(isec)=-1
152  EXIT
153  ENDIF
154  ENDIF
155  pt = dep
156  iseg = 0
157  dist=(x(dep)-x(arr))**2+(y(dep)-y(arr))**2
158 10 CONTINUE
159 !
160  DO ielem =1,nelem
161 !
162  i1 = ikle(ielem,1)
163  i2 = ikle(ielem,2)
164  i3 = ikle(ielem,3)
165 ! IF THE ELEMENT CONTAINS THE CURRENT POINT:
166  IF(pt.EQ.i1.OR.pt.EQ.i2.OR.pt.EQ.i3) THEN
167  dist1 = (x(i1)-x(arr))**2 + (y(i1)-y(arr))**2
168  dist2 = (x(i2)-x(arr))**2 + (y(i2)-y(arr))**2
169  dist3 = (x(i3)-x(arr))**2 + (y(i3)-y(arr))**2
170  IF(dist1.LT.dist) THEN
171  dist = dist1
172  elbest = ielem
173  igbest = i1
174  ilbest = 1
175  IF(i1.EQ.pt) ilprec = 1
176  IF(i2.EQ.pt) ilprec = 2
177  IF(i3.EQ.pt) ilprec = 3
178  ENDIF
179  IF(dist2.LT.dist) THEN
180  dist = dist2
181  elbest = ielem
182  igbest = i2
183  ilbest = 2
184  IF(i1.EQ.pt) ilprec = 1
185  IF(i2.EQ.pt) ilprec = 2
186  IF(i3.EQ.pt) ilprec = 3
187  ENDIF
188  IF(dist3.LT.dist) THEN
189  dist = dist3
190  elbest = ielem
191  igbest = i3
192  ilbest = 3
193  IF(i1.EQ.pt) ilprec = 1
194  IF(i2.EQ.pt) ilprec = 2
195  IF(i3.EQ.pt) ilprec = 3
196  ENDIF
197  ENDIF
198 !
199  ENDDO !IELEM
200 !
201  IF(igbest.EQ.pt) THEN
202  WRITE(lu,33)
203 33 FORMAT(1x,'FLUSEC : ALGORITHM FAILED')
204  CALL plante(1)
205  stop
206  ELSE
207  pt = igbest
208  iseg = iseg + 1
209  IF(iseg.GT.nsemax) THEN
210  WRITE(lu,*) 'FLUSEC: TOO MANY SEGMENTS IN A '
211  WRITE(lu,*) ' SECTION. INCREASE NSEMAX'
212  CALL plante(1)
213  stop
214  ENDIF
215  liste(isec,iseg,1) = ikle(elbest,ilprec)
216  liste(isec,iseg,2) = ikle(elbest,ilbest)
217  IF(igbest.NE.arr) GO TO 10
218  ENDIF
219 !
220  nseg(isec) = iseg
221 !
222 !JAJ #### THIS PART TO BE DONE IN THE PARALLEL CASE; FILL LISTE
223 ! WITH READY SEGMENT CHAINS PROVIDED BY PARTEL: SEE READ_SECTIONS
224 ! NOTE: FUTURE OPTIMISATION - USE CHAIN STRUCTURE IN THE WHOLE ROUTINE
225 !
226  ELSE
227 !
228  nseg(isec) = chain(isec)%NSEG
229  liste(isec,:,:)=0
230  DO iseg=1,nseg(isec)
231  liste(isec,iseg,1) = chain(isec)%LISTE(iseg,1)
232  liste(isec,iseg,2) = chain(isec)%LISTE(iseg,2)
233  END DO
234  WRITE(lu,*) 'CHAIN@SISYPHE -> LISTE@SISYPHE:'
235  WRITE(lu,*) 'ISEC,NSEG(ISEC): ',isec,nseg(isec)
236  DO iseg=1,nseg(isec)
237  WRITE(lu,*) liste(isec,iseg,:)
238  END DO
239  ENDIF
240 !
241  ENDDO
242 !
243 ! IF(.NOT.DEJA_FLUSEC) THEN
244  ENDIF
245 !
246 !-----------------------------------------------------------------------
247 !
248  deja_flusec = .true.
249 !
250 !-----------------------------------------------------------------------
251 !
252  DO isec = 1 , nsec
253 !
254  flx(isec) = 0.d0
255  flxs(isec) = 0.d0
256  flxc(isec) = 0.d0
257 !
258  IF(nseg(isec).GE.1) THEN
259 !
260 ! COMPUTES THE FLUX DIRECTLY, REGARDLESS OF THE WEAK FORM
261 ! OF THE IMPERMEABILITY CONDITION
262 !
263  DO iseg = 1 , nseg(isec)
264  i1 = liste(isec,iseg,1)
265  i2 = liste(isec,iseg,2)
266  x1 = x(i1)
267  x2 = x(i2)
268  y1 = y(i1)
269  y2 = y(i2)
270  h1 = h%R(i1)
271  h2 = h%R(i2)
272  nx = y1-y2
273  ny = x2-x1
274  un1= u%R(i1)*nx + v%R(i1)*ny
275  un2= u%R(i2)*nx + v%R(i2)*ny
276  flx(isec) = flx(isec) + ((h1+h2)*(un1+un2)+h2*un2+h1*un1)*sur6
277  IF(susp) THEN
278  un1= qsxs%R(i1)*nx + qsys%R(i1)*ny
279  un2= qsxs%R(i2)*nx + qsys%R(i2)*ny
280  flxs(isec) = flxs(isec) + 0.5d0*(un1+un2)
281  ENDIF
282  IF(charr) THEN
283  un1= qsxc%R(i1)*nx + qsyc%R(i1)*ny
284  un2= qsxc%R(i2)*nx + qsyc%R(i2)*ny
285  flxc(isec) = flxc(isec) + 0.5d0*(un1+un2)
286  ENDIF
287  ENDDO
288 !
289  IF(flx(isec).GT.0.d0) THEN
290  volpos(isec) = volpos(isec) + flx(isec)*dt
291  ELSE
292  volneg(isec) = volneg(isec) + flx(isec)*dt
293  ENDIF
294 !
295  IF(susp) THEN
296  IF(flxs(isec).GT.0.d0) THEN
297  volposs(isec) = volposs(isec) + flxs(isec)*dt
298  ELSE
299  volnegs(isec) = volnegs(isec) + flxs(isec)*dt
300  ENDIF
301  ENDIF
302 !
303  IF(charr) THEN
304  IF(flxc(isec).GT.0.d0) THEN
305  volposc(isec) = volposc(isec) + flxc(isec)*dt
306  ELSE
307  volnegc(isec) = volnegc(isec) + flxc(isec)*dt
308  ENDIF
309  ENDIF
310 !
311 ! IF(NSEG(ISEC).GT.1)...
312  ENDIF
313 !
314  ENDDO
315 !
316 !-----------------------------------------------------------------------
317 !
318 ! PRINTS THE RESULTS / !JAJ HERE ALLREDUCES FOR VALUES
319 !
320  CALL fluxpr_sisyphe(nsec,ctrlsc,flx,volneg,volpos,info,tps,
321  & nseg,ncsize,
322  & flxs,volnegs,volposs,susp,
323  & flxc,volnegc,volposc,charr)
324 !
325 !-----------------------------------------------------------------------
326 !
327 ! WRITE(LU,*) '-> LEAVING FLUSEC_SISYPHE'
328  RETURN
329  END
fluxpr_sisyphe
subroutine fluxpr_sisyphe(NSEC, CTRLSC, FLX, VOLNEG, VOLPOS, INFO, TPS, NSEG, NCSIZE, FLXS, VOLNEGS, VOLPOSS, SUSP, FLXC, VOLNEGC, VOLPOSC, CHARR)
Definition: fluxpr_sisyphe.f:8
declarations_special
Definition: declarations_special.F:3
declarations_sisyphe::old_method_flusec
logical old_method_flusec
Definition: declarations_sisyphe.f:1469
flusec_sisyphe
subroutine flusec_sisyphe(U, V, H, QSXC, QSYC, CHARR, QSXS, QSYS, SUSP, IKLE, NELMAX, NELEM, X, Y, DT, NCP, CTRLSC, INFO, TPS)
Definition: flusec_sisyphe.f:8
declarations_special::lu
integer lu
Definition: declarations_special.F:45
global_to_local_point
integer function global_to_local_point(IPOIN, MESH)
Definition: global_to_local_point.f:7
declarations_sisyphe
Definition: declarations_sisyphe.f:3
interface_sisyphe
Definition: interface_sisyphe.f:3
declarations_sisyphe::chain
type(chain_type), dimension(:), allocatable chain
Definition: declarations_sisyphe.f:1436
declarations_sisyphe::mesh
type(bief_mesh), target mesh
Definition: declarations_sisyphe.f:1354
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:34:51 for The TELEMAC-MASCARET system by doxygen 1.8.13