The TELEMAC-MASCARET system  trunk
flusec_gaia.f
Go to the documentation of this file.
1 ! *************************
2  SUBROUTINE flusec_gaia
3 ! *************************
4 !
5  &(u,v,h,qsxc,qsyc,charr,susp,
6  & ikle,nelmax,nelem,x,y,dt,ncp,ctrlsc,info,tps)
7 !
8 !***********************************************************************
9 ! GAIA V7P0 21/07/2011
10 !***********************************************************************
11 !
16 !
17 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
35 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36 !
37  USE bief
39  & liste,deja_flusec,nseg,volneg,
40  & volpos,flx,flxs,volnegs,volposs,
41  & flxc,volnegc,volposc
42  USE interface_gaia, ex_flusec_gaia => flusec_gaia
43 !
45  IMPLICIT NONE
46 !
47 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
48 !
49  INTEGER, INTENT(IN) :: NELMAX,NELEM,NCP
50  INTEGER, INTENT(IN) :: IKLE(nelmax,*)
51  INTEGER, INTENT(IN) :: CTRLSC(ncp)
52  DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),TPS,DT
53  LOGICAL, INTENT(IN) :: INFO,SUSP,CHARR
54  TYPE(bief_obj), INTENT(IN) :: U,V,H,QSXC,QSYC
55 !
56 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
57 !
58  INTEGER ERR
59  INTEGER, PARAMETER :: NSEMAX = 100
60 !
61  INTEGER IELEM,I1,I2,I3,ELBEST,IGBEST,ILBEST
62  INTEGER ILPREC,ISEG,ISEC,NSEC,PT,DEP,ARR
63 !
64  DOUBLE PRECISION DIST,DIST1,DIST2,DIST3
65  DOUBLE PRECISION H1,H2,X1,Y1,X2,Y2,UN1,UN2,NX,NY,SUR6
66 !
67 !
68 !-----------------------------------------------------------------------
69 !
70  sur6 = 1.d0/6.d0
71  nsec = ncp/2
72 !
73 ! LOOKS FOR WAYS THAT JOIN THE POINT COUPLES:
74 !
75  IF(.NOT.deja_flusec) THEN
76 !
77 ! DYNAMICALLY ALLOCATES FLX, VOLNEG, VOLPOS, ETC.
78 !
79  ALLOCATE(flx(nsec) ,stat=err)
80  ALLOCATE(volneg(nsec) ,stat=err)
81  ALLOCATE(volpos(nsec) ,stat=err)
82  ALLOCATE(nseg(ncp) ,stat=err)
83  ALLOCATE(liste(ncp,nsemax,2) ,stat=err)
84 ! S FOR SUSPENSION, C FOR BEDLOAD
85  ALLOCATE(flxs(nsec) ,stat=err)
86  ALLOCATE(volnegs(nsec) ,stat=err)
87  ALLOCATE(volposs(nsec) ,stat=err)
88  ALLOCATE(flxc(nsec) ,stat=err)
89  ALLOCATE(volnegc(nsec) ,stat=err)
90  ALLOCATE(volposc(nsec) ,stat=err)
91 !
92  IF(err.NE.0) THEN
93  WRITE(lu,200) err
94 200 FORMAT(1x,'FLUSEC: ERROR DURING ALLOCATION OF MEMORY: ',/,1x,
95  & 'ERROR CODE: ',1i6)
96  ENDIF
97 !
98  IF (.NOT.ALLOCATED(chain)) old_method_flusec=.true.
99 !
100  DO isec=1,nsec
101  flx(isec)=0.0d0
102  volneg(isec) =0.d0
103  volpos(isec) =0.d0
104  volnegs(isec)=0.d0
105  volposs(isec)=0.d0
106  volnegc(isec)=0.d0
107  volposc(isec)=0.d0
108  ENDDO
109 !
110  DO isec =1,nsec
111 !
112 !!JAJ #### IN THE SERIAL CASE, OR "CLASSICAL" IN PARALLEL,
113 ! FOLLOW THE ALGORITHM OF FINDING SEGMENT CHAINS
114 !
115 ! NOTE: IF YOU CHANGE THE ALGORITHM, CHANGE IT IN PARTEL AS WELL
116 !
117  IF (ncsize.LE.1 .OR. old_method_flusec) THEN
118 !
119  dep = ctrlsc(1+2*(isec-1))
120  arr = ctrlsc(2+2*(isec-1))
121  IF(ncsize.GT.1) THEN
122  dep=global_to_local_point(dep,mesh)
123  arr=global_to_local_point(arr,mesh)
124  IF(dep.EQ.0.AND.arr.EQ.0) THEN
125  nseg(isec)=0
126  EXIT
127  ENDIF
128  IF((dep.EQ.0.AND.arr.NE.0).OR.(dep.NE.0.AND.arr.EQ.0)) THEN
129  nseg(isec)=-1
130  EXIT
131  ENDIF
132  ENDIF
133  pt = dep
134  iseg = 0
135  dist=(x(dep)-x(arr))**2+(y(dep)-y(arr))**2
136 10 CONTINUE
137 !
138  DO ielem =1,nelem
139 !
140  i1 = ikle(ielem,1)
141  i2 = ikle(ielem,2)
142  i3 = ikle(ielem,3)
143 ! IF THE ELEMENT CONTAINS THE CURRENT POINT:
144  IF(pt.EQ.i1.OR.pt.EQ.i2.OR.pt.EQ.i3) THEN
145  dist1 = (x(i1)-x(arr))**2 + (y(i1)-y(arr))**2
146  dist2 = (x(i2)-x(arr))**2 + (y(i2)-y(arr))**2
147  dist3 = (x(i3)-x(arr))**2 + (y(i3)-y(arr))**2
148  IF(dist1.LT.dist) THEN
149  dist = dist1
150  elbest = ielem
151  igbest = i1
152  ilbest = 1
153  IF(i1.EQ.pt) ilprec = 1
154  IF(i2.EQ.pt) ilprec = 2
155  IF(i3.EQ.pt) ilprec = 3
156  ENDIF
157  IF(dist2.LT.dist) THEN
158  dist = dist2
159  elbest = ielem
160  igbest = i2
161  ilbest = 2
162  IF(i1.EQ.pt) ilprec = 1
163  IF(i2.EQ.pt) ilprec = 2
164  IF(i3.EQ.pt) ilprec = 3
165  ENDIF
166  IF(dist3.LT.dist) THEN
167  dist = dist3
168  elbest = ielem
169  igbest = i3
170  ilbest = 3
171  IF(i1.EQ.pt) ilprec = 1
172  IF(i2.EQ.pt) ilprec = 2
173  IF(i3.EQ.pt) ilprec = 3
174  ENDIF
175  ENDIF
176 !
177  ENDDO !IELEM
178 !
179  IF(igbest.EQ.pt) THEN
180  WRITE(lu,33)
181 33 FORMAT(1x,'FLUSEC : ALGORITHM FAILED')
182  CALL plante(1)
183  stop
184  ELSE
185  pt = igbest
186  iseg = iseg + 1
187  IF(iseg.GT.nsemax) THEN
188  WRITE(lu,*) 'FLUSEC: TOO MANY SEGMENTS IN A '
189  WRITE(lu,*) ' SECTION. INCREASE NSEMAX'
190  CALL plante(1)
191  stop
192  ENDIF
193  liste(isec,iseg,1) = ikle(elbest,ilprec)
194  liste(isec,iseg,2) = ikle(elbest,ilbest)
195  IF(igbest.NE.arr) GO TO 10
196  ENDIF
197 !
198  nseg(isec) = iseg
199 !
200 !JAJ #### THIS PART TO BE DONE IN THE PARALLEL CASE; FILL LISTE
201 ! WITH READY SEGMENT CHAINS PROVIDED BY PARTEL: SEE READ_SECTIONS
202 ! NOTE: FUTURE OPTIMISATION - USE CHAIN STRUCTURE IN THE WHOLE ROUTINE
203 !
204  ELSE
205 !
206  nseg(isec) = chain(isec)%NSEG
207  liste(isec,:,:)=0
208  DO iseg=1,nseg(isec)
209  liste(isec,iseg,1) = chain(isec)%LISTE(iseg,1)
210  liste(isec,iseg,2) = chain(isec)%LISTE(iseg,2)
211  END DO
212  WRITE(lu,*) 'CHAIN@GAIA -> LISTE@GAIA:'
213  WRITE(lu,*) 'ISEC,NSEG(ISEC): ',isec,nseg(isec)
214  DO iseg=1,nseg(isec)
215  WRITE(lu,*) liste(isec,iseg,:)
216  END DO
217  ENDIF
218 !
219  ENDDO
220 !
221 ! IF(.NOT.DEJA_FLUSEC) THEN
222  ENDIF
223 !
224 !-----------------------------------------------------------------------
225 !
226  deja_flusec = .true.
227 !
228 !-----------------------------------------------------------------------
229 !
230  DO isec = 1 , nsec
231 !
232  flx(isec) = 0.d0
233  flxs(isec) = 0.d0
234  flxc(isec) = 0.d0
235 !
236  IF(nseg(isec).GE.1) THEN
237 !
238 ! COMPUTES THE FLUX DIRECTLY, REGARDLESS OF THE WEAK FORM
239 ! OF THE IMPERMEABILITY CONDITION
240 !
241  DO iseg = 1 , nseg(isec)
242  i1 = liste(isec,iseg,1)
243  i2 = liste(isec,iseg,2)
244  x1 = x(i1)
245  x2 = x(i2)
246  y1 = y(i1)
247  y2 = y(i2)
248  h1 = h%R(i1)
249  h2 = h%R(i2)
250  nx = y1-y2
251  ny = x2-x1
252  un1= u%R(i1)*nx + v%R(i1)*ny
253  un2= u%R(i2)*nx + v%R(i2)*ny
254  flx(isec) = flx(isec) + ((h1+h2)*(un1+un2)+h2*un2+h1*un1)*sur6
255  IF(charr) THEN
256  un1= qsxc%R(i1)*nx + qsyc%R(i1)*ny
257  un2= qsxc%R(i2)*nx + qsyc%R(i2)*ny
258  flxc(isec) = flxc(isec) + 0.5d0*(un1+un2)
259  ENDIF
260  ENDDO
261 !
262  IF(flx(isec).GT.0.d0) THEN
263  volpos(isec) = volpos(isec) + flx(isec)*dt
264  ELSE
265  volneg(isec) = volneg(isec) + flx(isec)*dt
266  ENDIF
267 !
268  IF(susp) THEN
269  IF(flxs(isec).GT.0.d0) THEN
270  volposs(isec) = volposs(isec) + flxs(isec)*dt
271  ELSE
272  volnegs(isec) = volnegs(isec) + flxs(isec)*dt
273  ENDIF
274  ENDIF
275 !
276  IF(charr) THEN
277  IF(flxc(isec).GT.0.d0) THEN
278  volposc(isec) = volposc(isec) + flxc(isec)*dt
279  ELSE
280  volnegc(isec) = volnegc(isec) + flxc(isec)*dt
281  ENDIF
282  ENDIF
283 !
284 ! IF(NSEG(ISEC).GT.1)...
285  ENDIF
286 !
287  ENDDO
288 !
289 !-----------------------------------------------------------------------
290 !
291 ! PRINTS THE RESULTS / !JAJ HERE ALLREDUCES FOR VALUES
292 !
293  CALL fluxpr_gaia(nsec,ctrlsc,flx,volneg,volpos,info,tps,
294  & nseg,ncsize,
295  & flxs,volnegs,volposs,susp,
296  & flxc,volnegc,volposc,charr)
297 !
298 !-----------------------------------------------------------------------
299 !
300 ! WRITE(LU,*) '-> LEAVING FLUSEC_GAIA'
301  RETURN
302  END
subroutine flusec_gaia(U, V, H, QSXC, QSYC, CHARR, SUSP, IKLE, NELMAX, NELEM, X, Y, DT, NCP, CTRLSC, INFO, TPS)
Definition: flusec_gaia.f:8
subroutine fluxpr_gaia(NSEC, CTRLSC, FLX, VOLNEG, VOLPOS, INFO, TPS, NSEG, NCSIZE, FLXS, VOLNEGS, VOLPOSS, SUSP, FLXC, VOLNEGC, VOLPOSC, CHARR)
Definition: fluxpr_gaia.f:8
integer function global_to_local_point(IPOIN, MESH)
type(bief_mesh), target mesh
Mesh structure.
type(chain_type), dimension(:), allocatable chain
Definition: bief.f:3