The TELEMAC-MASCARET system  trunk
mt05pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt05pp
3 ! *****************
4 !
5  &(t,xm,xmul,su,sv,sw,u,v,w,f,g,
6  & x,y,z,ikle,nelem,nelmax,sigmag,specad,nplan)
7 !
8 !***********************************************************************
9 ! BIEF V6P3 21/08/2010
10 !***********************************************************************
11 !
12 !brief BUILDS THE MATVGR MATRIX.
13 !code
14 !+ COMPUTES THE COEFFICIENTS OF THE FOLLOWING MATRIX:
15 !+
16 !+
17 !+ / -> --->
18 !+ A(I,J) = XMUL / PSI1(I) * U . GRAD(PSI2(J)) D(OMEGA)
19 !+ /OMEGA
20 !+
21 !+
22 !+ BY ELEMENTARY CELL; THE ELEMENT IS THE P1 TRIANGLE
23 !+
24 !+ J(X,Y): JACOBIAN OF THE ISOPARAMETRIC TRANSFORMATION
25 !
26 !note JMH : 09/11/2004 : ACCORDING TO JMJ PRINCIPLE NOTE, SHOULD
27 !+ BE PERFORMED IN THE TRANSFORMED DOMAIN. VELOCITY W SHOULD HERE
28 !+ BE DELTA(Z)*WSTAR. DELTAZ IS H1, H2 AND H3 HERE.
29 !note IF SPECAD=TRUE, THE VELOCITY FIELD IS U+F*GRAD(G).
30 !
31 !warning THE JACOBIAN MUST BE POSITIVE
32 !
33 !history J-M HERVOUET (LNHE)
34 !+ 26/06/06
35 !+ V5P7
36 !+
37 !
38 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
39 !+ 13/07/2010
40 !+ V6P0
41 !+ Translation of French comments within the FORTRAN sources into
42 !+ English comments
43 !
44 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
45 !+ 21/08/2010
46 !+ V6P0
47 !+ Creation of DOXYGEN tags for automated documentation and
48 !+ cross-referencing of the FORTRAN sources
49 !
50 !history J-M HERVOUET (EDF R&D, LNHE)
51 !+ 11/01/2013
52 !+ V6P3
53 !+ XEL and YEL sent instead of XPT and YPT for X and Y.
54 !
55 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
56 !| F |-->| FUNCTION USED IN THE FORMULA
57 !| G |-->| FUNCTION USED IN THE FORMULA
58 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
59 !| H |-->| FUNCTION USED IN THE FORMULA
60 !| IKLE |-->| CONNECTIVITY TABLE.
61 !| NELEM |-->| NUMBER OF ELEMENTS
62 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
63 !| SF |-->| STRUCTURE OF FUNCTIONS F
64 !| SG |-->| STRUCTURE OF FUNCTIONS G
65 !| SH |-->| STRUCTURE OF FUNCTIONS H
66 !| SIGMAG |-->| IF YES, GENERALISED SIGMA TRANSFORMATION
67 !| SPECAD |-->| IF YES, SPECIAL ADVECTION FIELD (SEE IMPLEMENTATION)
68 !| SU |-->| STRUCTURE OF FUNCTIONS U
69 !| SV |-->| STRUCTURE OF FUNCTIONS V
70 !| SW |-->| STRUCTURE OF FUNCTIONS W
71 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
72 !| U |-->| FUNCTION USED IN THE FORMULA
73 !| V |-->| FUNCTION USED IN THE FORMULA
74 !| W |-->| FUNCTION USED IN THE FORMULA
75 !| X |-->| ABSCISSAE OF POINTS, PER ELEMENT
76 !| Y |-->| ORDINATES OF POINTS, PER ELEMENT
77 !| Z |-->| ELEVATIONS OF POINTS, PER POINTS !!!!
78 !| XM |<->| OFF-DIAGONAL TERMS
79 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
80 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
81 !
82  USE bief, ex_mt05pp => mt05pp
83 !
84  USE declarations_special
85  IMPLICIT NONE
86 !
87 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
88 !
89  INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN
90  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
91 !
92  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,30)
93 !
94  DOUBLE PRECISION, INTENT(IN) :: XMUL
95  DOUBLE PRECISION, INTENT(IN) :: U(*),V(*),W(*)
96  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*)
97 !
98 ! STRUCTURES OF U, V, W
99 !
100  TYPE(bief_obj), INTENT(IN) :: SU,SV,SW
101 !
102  DOUBLE PRECISION, INTENT(IN) :: X(nelmax,6),Y(nelmax,6),Z(*)
103  LOGICAL, INTENT(IN) :: SIGMAG,SPECAD
104 !
105 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
106 !
107 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
108 !
109  DOUBLE PRECISION INT1,X2,X3,Y2,Y3,DEN,DET,G2,G3,GRADGX,GRADGY
110  DOUBLE PRECISION PZ1,PX1,PX2,PX3,PY1,PY2,PY3
111  DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
112  DOUBLE PRECISION Q1,Q2,Q3,H1,H2,H3,HT
113  DOUBLE PRECISION SUS,SUI,SHUI,SHUS
114  DOUBLE PRECISION SHUI1,SHUS1,SHUI2,SHUS2,SHUI3,SHUS3
115  DOUBLE PRECISION HU1S3I,HU1I3S,HU1SI,HU2S3I,HU2I3S,HU2SI
116  DOUBLE PRECISION HU3S3I,HU3I3S,HU3SI
117  DOUBLE PRECISION SVS,SVI,SHVI,SHVS
118  DOUBLE PRECISION SHVI1,SHVS1,SHVI2,SHVS2,SHVI3,SHVS3
119  DOUBLE PRECISION HV1S3I,HV1I3S,HV1SI,HV2S3I,HV2I3S,HV2SI
120  DOUBLE PRECISION HV3S3I,HV3I3S,HV3SI,HH1,HH2,HH3,DZ1,DZ2,DZ3
121 !
122  INTEGER I1,I2,I3,I4,I5,I6,IELEM,IELEM2,IELEM3,NELEM2
123 !
124 !**********************************************************************
125 !
126  IF(su%ELM.NE.41) THEN
127  WRITE(lu,1001) su%ELM
128 1001 FORMAT(1x,'MT05PP (BIEF) : TYPE OF U NOT IMPLEMENTED: ',i6)
129  CALL plante(1)
130  stop
131  ENDIF
132  IF(sv%ELM.NE.41) THEN
133  WRITE(lu,2001) sv%ELM
134 2001 FORMAT(1x,'MT05PP (BIEF) : TYPE OF V NOT IMPLEMENTED: ',i6)
135  CALL plante(1)
136  stop
137  ENDIF
138 ! IN FACT W DEFINED PER LAYER BUT DECLARED 41
139  IF(sw%ELM.NE.41) THEN
140  WRITE(lu,3001) sw%ELM
141 3001 FORMAT(1x,'MT05PP (BIEF) : TYPE OF W NOT IMPLEMENTED: ',i6)
142  CALL plante(1)
143  stop
144  ENDIF
145 !
146 !-----------------------------------------------------------------------
147 !
148 ! WITH NORMAL SIGMA TRANSFORMATION
149 !
150  IF(.NOT.sigmag) THEN
151 !
152 ! LOOP ON THE ELEMENTS
153 !
154 ! STANDARD CASE
155 !
156  IF(.NOT.specad) THEN
157 !
158  DO ielem=1,nelem
159 !
160  i1 = ikle(ielem,1)
161  i2 = ikle(ielem,2)
162  i3 = ikle(ielem,3)
163  i4 = ikle(ielem,4)
164  i5 = ikle(ielem,5)
165  i6 = ikle(ielem,6)
166 !
167 ! X2 = X(I2) - X(I1)
168 ! X3 = X(I3) - X(I1)
169 ! Y2 = Y(I2) - Y(I1)
170 ! Y3 = Y(I3) - Y(I1)
171  x2 = x(ielem,2)
172  x3 = x(ielem,3)
173  y2 = y(ielem,2)
174  y3 = y(ielem,3)
175 !
176  u1 = u(i1)
177  u2 = u(i2)
178  u3 = u(i3)
179  u4 = u(i4)
180  u5 = u(i5)
181  u6 = u(i6)
182  v1 = v(i1)
183  v2 = v(i2)
184  v3 = v(i3)
185  v4 = v(i4)
186  v5 = v(i5)
187  v6 = v(i6)
188  q1 = w(i1)
189  q2 = w(i2)
190  q3 = w(i3)
191 !
192  h1 = z(i4)-z(i1)
193  h2 = z(i5)-z(i2)
194  h3 = z(i6)-z(i3)
195 !
196 ! INTERMEDIATE COMPUTATIONS
197 !
198  den=xmul/1440.d0
199  px2=y3*den
200  px3=-y2*den
201  px1=-px2-px3
202  py2=-x3*den
203  py3=x2*den
204  py1=-py3-py2
205 !
206  ht = h1+h2+h3
207  sus = u6+u5+u4
208  sui = u3+u2+u1
209  shui = h3*u3+h2*u2+h1*u1
210  shus = h3*u6+h2*u5+h1*u4
211  shus1 = (ht+h1)*(sus+u4)+shus+u4*h1
212  shui1 = (ht+h1)*(sui+u1)+shui+u1*h1
213  shus2 = (ht+h2)*(sus+u5)+shus+u5*h2
214  shui2 = (ht+h2)*(sui+u2)+shui+u2*h2
215  shus3 = (ht+h3)*(sus+u6)+shus+u6*h3
216  shui3 = (ht+h3)*(sui+u3)+shui+u3*h3
217  hu1s3i = shus1+3*shui1
218  hu1i3s = shui1+3*shus1
219  hu1si = shus1+shui1
220  hu2s3i = shus2+3*shui2
221  hu2i3s = shui2+3*shus2
222  hu2si = shus2+shui2
223  hu3s3i = shus3+3*shui3
224  hu3i3s = shui3+3*shus3
225  hu3si = shus3+shui3
226 !
227  svs = v6+v5+v4
228  svi = v3+v2+v1
229  shvi = h3*v3+h2*v2+h1*v1
230  shvs = h3*v6+h2*v5+h1*v4
231  shvs1 = (ht+h1)*(svs+v4)+shvs+v4*h1
232  shvi1 = (ht+h1)*(svi+v1)+shvi+v1*h1
233  shvs2 = (ht+h2)*(svs+v5)+shvs+v5*h2
234  shvi2 = (ht+h2)*(svi+v2)+shvi+v2*h2
235  shvs3 = (ht+h3)*(svs+v6)+shvs+v6*h3
236  shvi3 = (ht+h3)*(svi+v3)+shvi+v3*h3
237  hv1s3i = shvs1+3*shvi1
238  hv1i3s = shvi1+3*shvs1
239  hv1si = shvs1+ shvi1
240  hv2s3i = shvs2+3*shvi2
241  hv2i3s = shvi2+3*shvs2
242  hv2si = shvs2+ shvi2
243  hv3s3i = shvs3+3*shvi3
244  hv3i3s = shvi3+3*shvs3
245  hv3si = shvs3+ shvi3
246 !
247 ! OPTION WITH MASS-LUMPING ON VERTICAL VELOCITIES
248 ! COMPATIBILITY WITH TRIDW2 (LUMPING TO COMPUTE DELTA(Z)*WSTAR)
249 !
250  pz1=-den*(x2*y3-y2*x3)*30
251 !
252  t(ielem,1)=pz1*2*q1
253  xm(ielem,3) =-t(ielem,1)
254  xm(ielem,18)= t(ielem,1)
255  t(ielem,4) =-t(ielem,1)
256 !
257  t(ielem,2)=pz1*2*q2
258  xm(ielem,8) =-t(ielem,2)
259  xm(ielem,23)= t(ielem,2)
260  t(ielem,5) =-t(ielem,2)
261 !
262  t(ielem,3)=pz1*2*q3
263  xm(ielem,12)=-t(ielem,3)
264  xm(ielem,27)= t(ielem,3)
265  t(ielem,6) =-t(ielem,3)
266 !
267  xm(ielem,16)=pz1*q1
268  xm(ielem,7) =-xm(ielem,16)
269  xm(ielem,17)= xm(ielem,16)
270  xm(ielem,10)=-xm(ielem,16)
271  xm(ielem,19)= xm(ielem,16)
272  xm(ielem,28)=-xm(ielem,16)
273  xm(ielem,20)= xm(ielem,16)
274  xm(ielem,29)=-xm(ielem,16)
275 !
276  xm(ielem,1)=pz1*q2
277  xm(ielem,4) =-xm(ielem,1)
278  xm(ielem,22)= xm(ielem,1)
279  xm(ielem,13)=-xm(ielem,1)
280  xm(ielem,21)= xm(ielem,1)
281  xm(ielem,11)=-xm(ielem,1)
282  xm(ielem,24)= xm(ielem,1)
283  xm(ielem,30)=-xm(ielem,1)
284 !
285  xm(ielem,2)=pz1*q3
286  xm(ielem,5)=-xm(ielem,2)
287  xm(ielem,25)=xm(ielem,2)
288  xm(ielem,14)=-xm(ielem,2)
289  xm(ielem,26)=xm(ielem,2)
290  xm(ielem,15)=-xm(ielem,2)
291  xm(ielem,6)=xm(ielem,2)
292  xm(ielem,9)=-xm(ielem,2)
293 !
294  t(ielem,1)=t(ielem,1)+px1*hu1s3i+py1*hv1s3i
295  t(ielem,2)=t(ielem,2)+px2*hu2s3i+py2*hv2s3i
296  t(ielem,3)=t(ielem,3)+px3*hu3s3i+py3*hv3s3i
297  t(ielem,4)=t(ielem,4)+px1*hu1i3s+py1*hv1i3s
298  t(ielem,5)=t(ielem,5)+px2*hu2i3s+py2*hv2i3s
299  t(ielem,6)=t(ielem,6)+px3*hu3i3s+py3*hv3i3s
300 !
301  xm(ielem,01)=xm(ielem,01)+px2*hu1s3i+py2*hv1s3i
302  xm(ielem,02)=xm(ielem,02)+px3*hu1s3i+py3*hv1s3i
303  xm(ielem,03)=xm(ielem,03)+px1*hu1si +py1*hv1si
304  xm(ielem,04)=xm(ielem,04)+px2*hu1si +py2*hv1si
305  xm(ielem,05)=xm(ielem,05)+px3*hu1si +py3*hv1si
306  xm(ielem,06)=xm(ielem,06)+px3*hu2s3i+py3*hv2s3i
307  xm(ielem,07)=xm(ielem,07)+px1*hu2si +py1*hv2si
308  xm(ielem,08)=xm(ielem,08)+px2*hu2si +py2*hv2si
309  xm(ielem,09)=xm(ielem,09)+px3*hu2si +py3*hv2si
310  xm(ielem,10)=xm(ielem,10)+px1*hu3si +py1*hv3si
311  xm(ielem,11)=xm(ielem,11)+px2*hu3si +py2*hv3si
312  xm(ielem,12)=xm(ielem,12)+px3*hu3si +py3*hv3si
313  xm(ielem,13)=xm(ielem,13)+px2*hu1i3s+py2*hv1i3s
314  xm(ielem,14)=xm(ielem,14)+px3*hu1i3s+py3*hv1i3s
315  xm(ielem,15)=xm(ielem,15)+px3*hu2i3s+py3*hv2i3s
316  xm(ielem,16)=xm(ielem,16)+px1*hu2s3i+py1*hv2s3i
317  xm(ielem,17)=xm(ielem,17)+px1*hu3s3i+py1*hv3s3i
318  xm(ielem,18)=xm(ielem,18)+px1*hu1si +py1*hv1si
319  xm(ielem,19)=xm(ielem,19)+px1*hu2si +py1*hv2si
320  xm(ielem,20)=xm(ielem,20)+px1*hu3si +py1*hv3si
321  xm(ielem,21)=xm(ielem,21)+px2*hu3s3i+py2*hv3s3i
322  xm(ielem,22)=xm(ielem,22)+px2*hu1si +py2*hv1si
323  xm(ielem,23)=xm(ielem,23)+px2*hu2si +py2*hv2si
324  xm(ielem,24)=xm(ielem,24)+px2*hu3si +py2*hv3si
325  xm(ielem,25)=xm(ielem,25)+px3*hu1si +py3*hv1si
326  xm(ielem,26)=xm(ielem,26)+px3*hu2si +py3*hv2si
327  xm(ielem,27)=xm(ielem,27)+px3*hu3si +py3*hv3si
328  xm(ielem,28)=xm(ielem,28)+px1*hu2i3s+py1*hv2i3s
329  xm(ielem,29)=xm(ielem,29)+px1*hu3i3s+py1*hv3i3s
330  xm(ielem,30)=xm(ielem,30)+px2*hu3i3s+py2*hv3i3s
331 !
332  ENDDO
333 !
334  ELSE
335 !
336 ! WITH SPECIFIC ADVECTING FIELD
337 !
338  nelem2 = nelem/(nplan-1)
339 !
340  DO ielem=1,nelem
341 !
342  i1 = ikle(ielem,1)
343  i2 = ikle(ielem,2)
344  i3 = ikle(ielem,3)
345  i4 = ikle(ielem,4)
346  i5 = ikle(ielem,5)
347  i6 = ikle(ielem,6)
348 !
349 ! X2 = X(I2) - X(I1)
350 ! X3 = X(I3) - X(I1)
351 ! Y2 = Y(I2) - Y(I1)
352 ! Y3 = Y(I3) - Y(I1)
353  x2 = x(ielem,2)
354  x3 = x(ielem,3)
355  y2 = y(ielem,2)
356  y3 = y(ielem,3)
357 !
358  det= x2*y3-x3*y2
359  ielem2 = mod(ielem-1,nelem2) + 1
360 ! G IS PIECE-WISE LINEAR
361 ! IT IS ZCONV IN TELEMAC-3D
362  g2 = g(ielem2+nelem2)-g(ielem2)
363  g3 = g(ielem2+2*nelem2)-g(ielem2)
364  gradgx=(g2*y3-g3*y2)/det
365  gradgy=(x2*g3-x3*g2)/det
366 !
367  u1=u(i1)+f(i1)*gradgx
368  u2=u(i2)+f(i2)*gradgx
369  u3=u(i3)+f(i3)*gradgx
370  u4=u(i4)+f(i4)*gradgx
371  u5=u(i5)+f(i5)*gradgx
372  u6=u(i6)+f(i6)*gradgx
373  v1=v(i1)+f(i1)*gradgy
374  v2=v(i2)+f(i2)*gradgy
375  v3=v(i3)+f(i3)*gradgy
376  v4=v(i4)+f(i4)*gradgy
377  v5=v(i5)+f(i5)*gradgy
378  v6=v(i6)+f(i6)*gradgy
379 !
380  q1 = w(i1)
381  q2 = w(i2)
382  q3 = w(i3)
383 !
384  h1 = z(i4)-z(i1)
385  h2 = z(i5)-z(i2)
386  h3 = z(i6)-z(i3)
387 !
388 ! INTERMEDIATE COMPUTATIONS
389 !
390  den=xmul/1440.d0
391  px2=y3*den
392  px3=-y2*den
393  px1=-px2-px3
394  py2=-x3*den
395  py3=x2*den
396  py1=-py3-py2
397 !
398  ht = h1+h2+h3
399  sus = u6+u5+u4
400  sui = u3+u2+u1
401  shui = h3*u3+h2*u2+h1*u1
402  shus = h3*u6+h2*u5+h1*u4
403  shus1 = (ht+h1)*(sus+u4)+shus+u4*h1
404  shui1 = (ht+h1)*(sui+u1)+shui+u1*h1
405  shus2 = (ht+h2)*(sus+u5)+shus+u5*h2
406  shui2 = (ht+h2)*(sui+u2)+shui+u2*h2
407  shus3 = (ht+h3)*(sus+u6)+shus+u6*h3
408  shui3 = (ht+h3)*(sui+u3)+shui+u3*h3
409  hu1s3i = shus1+3*shui1
410  hu1i3s = shui1+3*shus1
411  hu1si = shus1+shui1
412  hu2s3i = shus2+3*shui2
413  hu2i3s = shui2+3*shus2
414  hu2si = shus2+shui2
415  hu3s3i = shus3+3*shui3
416  hu3i3s = shui3+3*shus3
417  hu3si = shus3+shui3
418 !
419  svs = v6+v5+v4
420  svi = v3+v2+v1
421  shvi = h3*v3+h2*v2+h1*v1
422  shvs = h3*v6+h2*v5+h1*v4
423  shvs1 = (ht+h1)*(svs+v4)+shvs+v4*h1
424  shvi1 = (ht+h1)*(svi+v1)+shvi+v1*h1
425  shvs2 = (ht+h2)*(svs+v5)+shvs+v5*h2
426  shvi2 = (ht+h2)*(svi+v2)+shvi+v2*h2
427  shvs3 = (ht+h3)*(svs+v6)+shvs+v6*h3
428  shvi3 = (ht+h3)*(svi+v3)+shvi+v3*h3
429  hv1s3i = shvs1+3*shvi1
430  hv1i3s = shvi1+3*shvs1
431  hv1si = shvs1+ shvi1
432  hv2s3i = shvs2+3*shvi2
433  hv2i3s = shvi2+3*shvs2
434  hv2si = shvs2+ shvi2
435  hv3s3i = shvs3+3*shvi3
436  hv3i3s = shvi3+3*shvs3
437  hv3si = shvs3+ shvi3
438 !
439 ! OPTION WITH MASS-LUMPING ON VERTICAL VELOCITIES
440 ! COMPATIBILITY WITH TRIDW2 (LUMPING TO COMPUTE DELTA(Z)*WSTAR)
441 !
442  pz1=-den*(x2*y3-y2*x3)*30
443 !
444  int1=pz1*2*q1
445  t(ielem,1)=int1
446  xm(ielem,3)=-int1
447  xm(ielem,18)=int1
448  t(ielem,4)=-int1
449 !
450  int1=pz1*2*q2
451  t(ielem,2)=int1
452  xm(ielem,8)= -int1
453  xm(ielem,23)=int1
454  t(ielem,5)=-int1
455 !
456  int1=pz1*2*q3
457  t(ielem,3)=int1
458  xm(ielem,12)=-int1
459  xm(ielem,27)=int1
460  t(ielem,6)=-int1
461 !
462  int1=pz1*q1
463  xm(ielem,16)=int1
464  xm(ielem,7)=-int1
465  xm(ielem,17)=int1
466  xm(ielem,10)=-int1
467  xm(ielem,19)=int1
468  xm(ielem,28)=-int1
469  xm(ielem,20)=int1
470  xm(ielem,29)=-int1
471 !
472  int1=pz1*q2
473  xm(ielem,1)=int1
474  xm(ielem,4)=-int1
475  xm(ielem,22)=int1
476  xm(ielem,13)=-int1
477  xm(ielem,21)=int1
478  xm(ielem,11)=-int1
479  xm(ielem,24)=int1
480  xm(ielem,30)=-int1
481 !
482  int1=pz1*q3
483  xm(ielem,2)=int1
484  xm(ielem,5)=-int1
485  xm(ielem,25)=int1
486  xm(ielem,14)=-int1
487  xm(ielem,26)=int1
488  xm(ielem,15)=-int1
489  xm(ielem,6)=int1
490  xm(ielem,9)=-int1
491 !
492  t(ielem,1)=t(ielem,1)+px1*hu1s3i+py1*hv1s3i
493  t(ielem,2)=t(ielem,2)+px2*hu2s3i+py2*hv2s3i
494  t(ielem,3)=t(ielem,3)+px3*hu3s3i+py3*hv3s3i
495  t(ielem,4)=t(ielem,4)+px1*hu1i3s+py1*hv1i3s
496  t(ielem,5)=t(ielem,5)+px2*hu2i3s+py2*hv2i3s
497  t(ielem,6)=t(ielem,6)+px3*hu3i3s+py3*hv3i3s
498 !
499  xm(ielem,01)=xm(ielem,01)+px2*hu1s3i+py2*hv1s3i
500  xm(ielem,02)=xm(ielem,02)+px3*hu1s3i+py3*hv1s3i
501  xm(ielem,03)=xm(ielem,03)+px1*hu1si +py1*hv1si
502  xm(ielem,04)=xm(ielem,04)+px2*hu1si +py2*hv1si
503  xm(ielem,05)=xm(ielem,05)+px3*hu1si +py3*hv1si
504  xm(ielem,06)=xm(ielem,06)+px3*hu2s3i+py3*hv2s3i
505  xm(ielem,07)=xm(ielem,07)+px1*hu2si +py1*hv2si
506  xm(ielem,08)=xm(ielem,08)+px2*hu2si +py2*hv2si
507  xm(ielem,09)=xm(ielem,09)+px3*hu2si +py3*hv2si
508  xm(ielem,10)=xm(ielem,10)+px1*hu3si +py1*hv3si
509  xm(ielem,11)=xm(ielem,11)+px2*hu3si +py2*hv3si
510  xm(ielem,12)=xm(ielem,12)+px3*hu3si +py3*hv3si
511  xm(ielem,13)=xm(ielem,13)+px2*hu1i3s+py2*hv1i3s
512  xm(ielem,14)=xm(ielem,14)+px3*hu1i3s+py3*hv1i3s
513  xm(ielem,15)=xm(ielem,15)+px3*hu2i3s+py3*hv2i3s
514  xm(ielem,16)=xm(ielem,16)+px1*hu2s3i+py1*hv2s3i
515  xm(ielem,17)=xm(ielem,17)+px1*hu3s3i+py1*hv3s3i
516  xm(ielem,18)=xm(ielem,18)+px1*hu1si +py1*hv1si
517  xm(ielem,19)=xm(ielem,19)+px1*hu2si +py1*hv2si
518  xm(ielem,20)=xm(ielem,20)+px1*hu3si +py1*hv3si
519  xm(ielem,21)=xm(ielem,21)+px2*hu3s3i+py2*hv3s3i
520  xm(ielem,22)=xm(ielem,22)+px2*hu1si +py2*hv1si
521  xm(ielem,23)=xm(ielem,23)+px2*hu2si +py2*hv2si
522  xm(ielem,24)=xm(ielem,24)+px2*hu3si +py2*hv3si
523  xm(ielem,25)=xm(ielem,25)+px3*hu1si +py3*hv1si
524  xm(ielem,26)=xm(ielem,26)+px3*hu2si +py3*hv2si
525  xm(ielem,27)=xm(ielem,27)+px3*hu3si +py3*hv3si
526  xm(ielem,28)=xm(ielem,28)+px1*hu2i3s+py1*hv2i3s
527  xm(ielem,29)=xm(ielem,29)+px1*hu3i3s+py1*hv3i3s
528  xm(ielem,30)=xm(ielem,30)+px2*hu3i3s+py2*hv3i3s
529 !
530  ENDDO
531 !
532  ENDIF
533 !
534 !-----------------------------------------------------------------------
535 !
536 ! WITH GENERALIZED SIGMA TRANSFORMATION
537 !
538  ELSE
539 !
540  IF(specad) THEN
541  WRITE(lu,202)
542 202 FORMAT(1x,'MT05PP (BIEF) :',/,
543  & 1x,'SIGMAG=T AND SPECAD=T : ',1i6,' NOT IMPLEMENTED',/,
544  & 1x,'REAL NAME OF U : ',a6)
545  CALL plante(1)
546  stop
547  ENDIF
548 !
549  nelem2 = nelem/(nplan-1)
550 !
551  DO ielem=1,nelem
552 !
553  ielem2 = mod(ielem-1,nelem2) + 1
554  ielem3 = nelem - nelem2 + ielem2
555 !
556  i1 = ikle(ielem,1)
557  i2 = ikle(ielem,2)
558  i3 = ikle(ielem,3)
559  i4 = ikle(ielem,4)
560  i5 = ikle(ielem,5)
561  i6 = ikle(ielem,6)
562 !
563 ! COMPUTES DELTA Z
564 !
565  dz1 = z(i4) - z(i1)
566  dz2 = z(i5) - z(i2)
567  dz3 = z(i6) - z(i3)
568 !
569 ! COMPUTES WATER DEPTH
570 !
571  h1 = z(ikle(ielem3,4))-z(ikle(ielem2,1))
572  h2 = z(ikle(ielem3,5))-z(ikle(ielem2,2))
573  h3 = z(ikle(ielem3,6))-z(ikle(ielem2,3))
574 !
575  hh1=max(h1,1.d-6)
576  hh2=max(h2,1.d-6)
577  hh3=max(h3,1.d-6)
578 !
579 ! X2 = X(I2) - X(I1)
580 ! X3 = X(I3) - X(I1)
581 ! Y2 = Y(I2) - Y(I1)
582 ! Y3 = Y(I3) - Y(I1)
583  x2 = x(ielem,2)
584  x3 = x(ielem,3)
585  y2 = y(ielem,2)
586  y3 = y(ielem,3)
587 !
588  u1=dz1*u(i1)/hh1
589  u2=dz2*u(i2)/hh2
590  u3=dz3*u(i3)/hh3
591  u4=dz1*u(i4)/hh1
592  u5=dz2*u(i5)/hh2
593  u6=dz3*u(i6)/hh3
594  v1=dz1*v(i1)/hh1
595  v2=dz2*v(i2)/hh2
596  v3=dz3*v(i3)/hh3
597  v4=dz1*v(i4)/hh1
598  v5=dz2*v(i5)/hh2
599  v6=dz3*v(i6)/hh3
600 !
601  q1 = w(i1)
602  q2 = w(i2)
603  q3 = w(i3)
604 !
605 ! INTERMEDIATE COMPUTATIONS
606 !
607  den=xmul/1440.d0
608  px2=y3*den
609  px3=-y2*den
610  px1=-px2-px3
611  py2=-x3*den
612  py3=x2*den
613  py1=-py3-py2
614 !
615  ht = h1+h2+h3
616  sus = u6+u5+u4
617  sui = u3+u2+u1
618  shui = h3*u3+h2*u2+h1*u1
619  shus = h3*u6+h2*u5+h1*u4
620  shus1 = (ht+h1)*(sus+u4)+shus+u4*h1
621  shui1 = (ht+h1)*(sui+u1)+shui+u1*h1
622  shus2 = (ht+h2)*(sus+u5)+shus+u5*h2
623  shui2 = (ht+h2)*(sui+u2)+shui+u2*h2
624  shus3 = (ht+h3)*(sus+u6)+shus+u6*h3
625  shui3 = (ht+h3)*(sui+u3)+shui+u3*h3
626  hu1s3i = shus1+3*shui1
627  hu1i3s = shui1+3*shus1
628  hu1si = shus1+shui1
629  hu2s3i = shus2+3*shui2
630  hu2i3s = shui2+3*shus2
631  hu2si = shus2+shui2
632  hu3s3i = shus3+3*shui3
633  hu3i3s = shui3+3*shus3
634  hu3si = shus3+shui3
635 !
636  svs = v6+v5+v4
637  svi = v3+v2+v1
638  shvi = h3*v3+h2*v2+h1*v1
639  shvs = h3*v6+h2*v5+h1*v4
640  shvs1 = (ht+h1)*(svs+v4)+shvs+v4*h1
641  shvi1 = (ht+h1)*(svi+v1)+shvi+v1*h1
642  shvs2 = (ht+h2)*(svs+v5)+shvs+v5*h2
643  shvi2 = (ht+h2)*(svi+v2)+shvi+v2*h2
644  shvs3 = (ht+h3)*(svs+v6)+shvs+v6*h3
645  shvi3 = (ht+h3)*(svi+v3)+shvi+v3*h3
646  hv1s3i = shvs1+3*shvi1
647  hv1i3s = shvi1+3*shvs1
648  hv1si = shvs1+ shvi1
649  hv2s3i = shvs2+3*shvi2
650  hv2i3s = shvi2+3*shvs2
651  hv2si = shvs2+ shvi2
652  hv3s3i = shvs3+3*shvi3
653  hv3i3s = shvi3+3*shvs3
654  hv3si = shvs3+ shvi3
655 !
656 ! OPTION WITH MASS-LUMPING ON VERTICAL VELOCITIES
657 ! COMPATIBILITY WITH TRIDW2 (LUMPING TO COMPUTE DELTA(Z)*WSTAR)
658 !
659  pz1=-den*(x2*y3-y2*x3)*30
660 !
661  int1=pz1*2*q1
662  t(ielem,1)=int1
663  xm(ielem,3)=-int1
664  xm(ielem,18)=int1
665  t(ielem,4)=-int1
666 !
667  int1=pz1*2*q2
668  t(ielem,2)=int1
669  xm(ielem,8)= -int1
670  xm(ielem,23)=int1
671  t(ielem,5)=-int1
672 !
673  int1=pz1*2*q3
674  t(ielem,3)=int1
675  xm(ielem,12)=-int1
676  xm(ielem,27)=int1
677  t(ielem,6)=-int1
678 !
679  int1=pz1*q1
680  xm(ielem,16)=int1
681  xm(ielem,7)=-int1
682  xm(ielem,17)=int1
683  xm(ielem,10)=-int1
684  xm(ielem,19)=int1
685  xm(ielem,28)=-int1
686  xm(ielem,20)=int1
687  xm(ielem,29)=-int1
688 !
689  int1=pz1*q2
690  xm(ielem,1)=int1
691  xm(ielem,4)=-int1
692  xm(ielem,22)=int1
693  xm(ielem,13)=-int1
694  xm(ielem,21)=int1
695  xm(ielem,11)=-int1
696  xm(ielem,24)=int1
697  xm(ielem,30)=-int1
698 !
699  int1=pz1*q3
700  xm(ielem,2)=int1
701  xm(ielem,5)=-int1
702  xm(ielem,25)=int1
703  xm(ielem,14)=-int1
704  xm(ielem,26)=int1
705  xm(ielem,15)=-int1
706  xm(ielem,6)=int1
707  xm(ielem,9)=-int1
708 !
709  t(ielem,1)=t(ielem,1)+px1*hu1s3i+py1*hv1s3i
710  t(ielem,2)=t(ielem,2)+px2*hu2s3i+py2*hv2s3i
711  t(ielem,3)=t(ielem,3)+px3*hu3s3i+py3*hv3s3i
712  t(ielem,4)=t(ielem,4)+px1*hu1i3s+py1*hv1i3s
713  t(ielem,5)=t(ielem,5)+px2*hu2i3s+py2*hv2i3s
714  t(ielem,6)=t(ielem,6)+px3*hu3i3s+py3*hv3i3s
715 !
716  xm(ielem,01)=xm(ielem,01)+px2*hu1s3i+py2*hv1s3i
717  xm(ielem,02)=xm(ielem,02)+px3*hu1s3i+py3*hv1s3i
718  xm(ielem,03)=xm(ielem,03)+px1*hu1si +py1*hv1si
719  xm(ielem,04)=xm(ielem,04)+px2*hu1si +py2*hv1si
720  xm(ielem,05)=xm(ielem,05)+px3*hu1si +py3*hv1si
721  xm(ielem,06)=xm(ielem,06)+px3*hu2s3i+py3*hv2s3i
722  xm(ielem,07)=xm(ielem,07)+px1*hu2si +py1*hv2si
723  xm(ielem,08)=xm(ielem,08)+px2*hu2si +py2*hv2si
724  xm(ielem,09)=xm(ielem,09)+px3*hu2si +py3*hv2si
725  xm(ielem,10)=xm(ielem,10)+px1*hu3si +py1*hv3si
726  xm(ielem,11)=xm(ielem,11)+px2*hu3si +py2*hv3si
727  xm(ielem,12)=xm(ielem,12)+px3*hu3si +py3*hv3si
728  xm(ielem,13)=xm(ielem,13)+px2*hu1i3s+py2*hv1i3s
729  xm(ielem,14)=xm(ielem,14)+px3*hu1i3s+py3*hv1i3s
730  xm(ielem,15)=xm(ielem,15)+px3*hu2i3s+py3*hv2i3s
731  xm(ielem,16)=xm(ielem,16)+px1*hu2s3i+py1*hv2s3i
732  xm(ielem,17)=xm(ielem,17)+px1*hu3s3i+py1*hv3s3i
733  xm(ielem,18)=xm(ielem,18)+px1*hu1si +py1*hv1si
734  xm(ielem,19)=xm(ielem,19)+px1*hu2si +py1*hv2si
735  xm(ielem,20)=xm(ielem,20)+px1*hu3si +py1*hv3si
736  xm(ielem,21)=xm(ielem,21)+px2*hu3s3i+py2*hv3s3i
737  xm(ielem,22)=xm(ielem,22)+px2*hu1si +py2*hv1si
738  xm(ielem,23)=xm(ielem,23)+px2*hu2si +py2*hv2si
739  xm(ielem,24)=xm(ielem,24)+px2*hu3si +py2*hv3si
740  xm(ielem,25)=xm(ielem,25)+px3*hu1si +py3*hv1si
741  xm(ielem,26)=xm(ielem,26)+px3*hu2si +py3*hv2si
742  xm(ielem,27)=xm(ielem,27)+px3*hu3si +py3*hv3si
743  xm(ielem,28)=xm(ielem,28)+px1*hu2i3s+py1*hv2i3s
744  xm(ielem,29)=xm(ielem,29)+px1*hu3i3s+py1*hv3i3s
745  xm(ielem,30)=xm(ielem,30)+px2*hu3i3s+py2*hv3i3s
746 !
747  ENDDO
748 !
749  ENDIF
750 !
751 !-----------------------------------------------------------------------
752 !
753  RETURN
754  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt05pp
subroutine mt05pp(T, XM, XMUL, SU, SV, SW, U, V, W, F, G, X, Y, Z, IKLE, NELEM, NELMAX, SIGMAG, SPECAD, NPLAN)
Definition: mt05pp.f:8
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:07 for The TELEMAC-MASCARET system by doxygen 1.8.13