The TELEMAC-MASCARET system  trunk
mt02pp_star.f
Go to the documentation of this file.
1 ! **********************
2  SUBROUTINE mt02pp_star
3 ! **********************
4 !
5  &(t,xm,xmul,sf,sg,sh,f,g,h,x,y,z,surfac,ikle,nelem,nelmax,inchyd,
6  & formul,nplan)
7 !
8 !***********************************************************************
9 ! BIEF V6P3 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE DIFFUSION MATRIX AS IN MT02PP BUT HERE
13 !+ WITH A PRIOR DECOMPOSITION OF GRADIENTS ALOG PLANES
14 !+ AND ON THE VERTICAL, THIS GIVES 4 TERMS WHICH ARE
15 !+ SUCCESSIVELY COMPUTED HERE, AND MAY BE COMPUTED
16 !+ SEPARATELY. AS THERE ARE APPROXIMATIONS WHICH ARE
17 !+ DIFFERENT, THE RESULT IS SLIGHTLY DIFFERENT FROM MT02PP
18 !+ THIS IS USED E.G. FOR CORRECTING HORIZONTAL AND VERTICAL
19 !+ TO GET DIVERGENCE FREE FLUXES.
20 !+
21 !+ THE FUNCTION DIFFUSION COEFFICIENT IS HERE A P1
22 !+ DIAGONAL TENSOR.
23 !+
24 !+ CASE OF THE PRISM.
25 !code
26 !+ STORAGE CONVENTION FOR EXTRA-DIAGONAL TERMS:
27 !+
28 !+ XM(IELEM, 1) ----> M(1,2) = M(2,1)
29 !+ XM(IELEM, 2) ----> M(1,3) = M(3,1)
30 !+ XM(IELEM, 3) ----> M(1,4) = M(4,1)
31 !+ XM(IELEM, 4) ----> M(1,5) = M(5,1)
32 !+ XM(IELEM, 5) ----> M(1,6) = M(6,1)
33 !+ XM(IELEM, 6) ----> M(2,3) = M(3,2)
34 !+ XM(IELEM, 7) ----> M(2,4) = M(4,2)
35 !+ XM(IELEM, 8) ----> M(2,5) = M(5,2)
36 !+ XM(IELEM, 9) ----> M(2,6) = M(6,2)
37 !+ XM(IELEM,10) ----> M(3,4) = M(4,3)
38 !+ XM(IELEM,11) ----> M(3,5) = M(5,3)
39 !+ XM(IELEM,12) ----> M(3,6) = M(6,3)
40 !+ XM(IELEM,13) ----> M(4,5) = M(5,4)
41 !+ XM(IELEM,14) ----> M(4,6) = M(6,4)
42 !+ XM(IELEM,15) ----> M(5,6) = M(6,5)
43 !
44 !warning THE JACOBIAN MUST BE POSITIVE
45 !
46 !history
47 !+ 22/08/07
48 !+
49 !+ CAN OPTIONALLY TREAT (USES DIMDISC) A P0 VERTICAL
50 !
51 !history
52 !+ 06/02/07
53 !+
54 !+ IF FORMUL(14:16)='MON' :
55 !
56 !history JMH
57 !+ 15/03/2010
58 !+
59 !+ PARAMETER EPSILON ADDED
60 !
61 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
62 !+ 20/05/2010
63 !+ V6P0
64 !+
65 !
66 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
67 !+ 13/07/2010
68 !+ V6P0
69 !+ Translation of French comments within the FORTRAN sources into
70 !+ English comments
71 !
72 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
73 !+ 21/08/2010
74 !+ V6P0
75 !+ Creation of DOXYGEN tags for automated documentation and
76 !+ cross-referencing of the FORTRAN sources
77 !
78 !history U.H.MErkel
79 !+ 18/07/2012
80 !+ V6P2
81 !+ Replaced EPSILON with CHOUIA due to nag compiler problems
82 !
83 !history J-M HERVOUET (LNHE)
84 !+ 11/01/2013
85 !+ V6P3
86 !+ XEL and YEL sent instead of XPT and YPT.
87 !
88 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
89 !| F |-->| FUNCTION USED IN THE FORMULA
90 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
91 !| G |-->| FUNCTION USED IN THE FORMULA
92 !| H |-->| FUNCTION USED IN THE FORMULA
93 !| IKLE |-->| CONNECTIVITY TABLE.
94 !| INCHYD |-->| IF YES, TREATS HYDROSTATIC INCONSISTENCIES
95 !| NELEM |-->| NUMBER OF ELEMENTS
96 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
97 !| NPLAN |-->| NUMBER OF PLANES IN THE MESH OF PRISMS
98 !| SF |-->| STRUCTURE OF FUNCTIONS F
99 !| SG |-->| STRUCTURE OF FUNCTIONS G
100 !| SH |-->| STRUCTURE OF FUNCTIONS H
101 !| SURFAC |-->| AREA OF 2D ELEMENTS
102 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
103 !| X |-->| ABSCISSAE OF POINTS, PER ELEMENT
104 !| Y |-->| ORDINATES OF POINTS, PER ELEMENT
105 !| Z |-->| ELEVATIONS OF POINTS, PER POINT
106 !| XM |<->| OFF-DIAGONAL TERMS
107 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
108 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
109 !
110  USE bief, ex_mt02pp_star => mt02pp_star
111 !
112  USE declarations_special
113  IMPLICIT NONE
114 !
115 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
116 !
117  INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN
118  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
119 !
120 ! 15 OR 30
121  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,*)
122  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
123 !
124  DOUBLE PRECISION, INTENT(IN) :: XMUL
125  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),H(*)
126 !
127 ! STRUCTURES DE F,G,H
128 !
129  TYPE(bief_obj), INTENT(IN) :: SF,SG,SH
130 !
131  DOUBLE PRECISION, INTENT(IN) :: X(nelmax,6),Y(nelmax,6),Z(*)
132 !
133  LOGICAL, INTENT(IN) :: INCHYD
134  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
135 !
136 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
137 !
138 ! DECLARATIONS SPECIFIQUES A CE SOUS-PROGRAMME
139 !
140  DOUBLE PRECISION H1,H2,H3
141  DOUBLE PRECISION NF1,NF2,NF3,NG1,NG2,NG3,NH1,NH2,NH3
142  DOUBLE PRECISION XS06,XS24,NPX,NPY,NPXL,NPXU,XS03
143  DOUBLE PRECISION NPYL,NPYU,GX(3),GY(3),NUXMOY,NUYMOY
144 !
145  DOUBLE PRECISION X2,X3,Y2,Y3,SOMVX,X2X3,X2AUX,X3AUX,NPX2,NPY2
146  DOUBLE PRECISION SOMVY,Y2Y3,Y2AUX,Y3AUX,AUX
147 !
148  INTEGER I1,I2,I3,I4,I5,I6,IELEM,II4,II5,II6,NPOU0
149 !
150  DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-4
151 !
152 !***********************************************************************
153 !
154  xs03 =xmul/3.d0
155  xs06 =xmul/6.d0
156  xs24 =xmul/24.d0
157 !
158  IF(sf%ELM.NE.41) THEN
159  WRITE(lu,1001) sf%ELM
160 1001 FORMAT(1x,'MT02PP_STAR (BIEF) : TYPE OF F NOT IMPLEMENTED: ',i6)
161  CALL plante(1)
162  stop
163  ENDIF
164  IF(sg%ELM.NE.41) THEN
165  WRITE(lu,2001) sg%ELM
166 2001 FORMAT(1x,'MT02PP_STAR (BIEF) : TYPE OF G NOT IMPLEMENTED: ',i6)
167  CALL plante(1)
168  stop
169  ENDIF
170  IF(sh%ELM.NE.41) THEN
171  WRITE(lu,3001) sh%ELM
172 3001 FORMAT(1x,'MT02PP_STAR (BIEF) : TYPE OF H NOT IMPLEMENTED: ',i6)
173  CALL plante(1)
174  stop
175  ENDIF
176 !
177 ! POUR LE TRAITEMENT DE LA VISCOSITE VERTICALE P0 SUR LA VERTICALE
178 ! VOIR VISCLM DE TELEMAC-3D
179 !
180  IF(sh%DIMDISC.EQ.0) THEN
181 ! VISCOSITE VERTICALE P1
182  npou0=sh%DIM1/nplan
183  ELSEIF(sh%DIMDISC.EQ.4111) THEN
184 ! VISCOSITE VERTICALE P0 SUR LA VERTICALE (VOIR II4,5,6)
185  npou0=0
186  ELSE
187  WRITE(lu,4001) sh%DIMDISC
188 4001 FORMAT(1x,'MT02PP_STAR (BIEF): DIMDISC NOT IMPLEMENTED: ',i6)
189  CALL plante(1)
190  stop
191  ENDIF
192 !
193  IF(formul(11:13).EQ.'MON') THEN
194 !
195  IF(xmul.LT.0.d0) THEN
196  WRITE(lu,4003) xmul
197 4003 FORMAT(1x,'MT02PP (BIEF): NEGATIVE XMUL EXCLUDED: ',g16.7)
198  CALL plante(1)
199  stop
200  ENDIF
201 !
202 !-----------------------------------------------------------------------
203 !
204 ! VERSION WITH MONOTONICITY
205 !
206 !-----------------------------------------------------------------------
207 !
208  DO ielem=1,nelem
209 !
210  i1=ikle(ielem,1)
211  i2=ikle(ielem,2)
212  i3=ikle(ielem,3)
213  i4=ikle(ielem,4)
214  i5=ikle(ielem,5)
215  i6=ikle(ielem,6)
216 ! SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
217  ii4=i1+npou0
218  ii5=i2+npou0
219  ii6=i3+npou0
220 !
221  h1=z(i4)-z(i1)
222  h2=z(i5)-z(i2)
223  h3=z(i6)-z(i3)
224 !
225  IF((inchyd.AND.max(z(i1),z(i2),z(i3)).GT.
226  & min(z(i4),z(i5),z(i6))) .OR.
227  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
228  nh1=0.d0
229  nh2=0.d0
230  nh3=0.d0
231  ELSE
232 ! SUIVANT LES CAS (II4=I1 OU I4, ETC.)
233 ! VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
234 ! VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
235  nh1=0.5d0*(h(i1)+h(ii4))/h1
236  nh2=0.5d0*(h(i2)+h(ii5))/h2
237  nh3=0.5d0*(h(i3)+h(ii6))/h3
238  ENDIF
239 !
240 ! DIFFUSION ALONG PLANES
241 !
242 ! X2 = X(I2)-X(I1)
243 ! X3 = X(I3)-X(I1)
244 ! Y2 = Y(I2)-Y(I1)
245 ! Y3 = Y(I3)-Y(I1)
246  x2 = x(ielem,2)
247  x3 = x(ielem,3)
248  y2 = y(ielem,2)
249  y3 = y(ielem,3)
250 !
251  x2x3 = x2*x3
252  y2y3 = y2*y3
253  x2aux = x2x3-x2**2
254  x3aux = x3**2-x2x3
255  y2aux = y2y3-y2**2
256  y3aux = y3**2-y2y3
257 !
258  aux = xs24 / surfac(ielem)
259 !
260 ! 2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
261 ! LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
262  somvx = ( f(i1)*h1+f(i2)*h2+f(i3)*h3 ) * aux
263  somvy = ( g(i1)*h1+g(i2)*h2+g(i3)*h3 ) * aux
264 !
265 ! OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
266 ! WITH MONOTONICITY ENSURED
267 !
268  xm(ielem, 1) = min(0.d0,- somvy * x3aux - somvx * y3aux)
269  xm(ielem, 2) = min(0.d0, somvy * x2aux + somvx * y2aux)
270  xm(ielem, 6) = min(0.d0,- somvy * x2x3 - somvx * y2y3 )
271 !
272 ! 2D DIFFUSION, UPPER LEVEL
273  somvx = ( f(i4)*h1+f(i5)*h2+f(i6)*h3 ) * aux
274  somvy = ( g(i4)*h1+g(i5)*h2+g(i6)*h3 ) * aux
275 !
276 ! OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
277 ! WITH MONOTONICITY ENSURED
278 !
279  xm(ielem,13) = min(0.d0,- somvy * x3aux - somvx * y3aux)
280  xm(ielem,14) = min(0.d0, somvy * x2aux + somvx * y2aux)
281  xm(ielem,15) = min(0.d0,- somvy * x2x3 - somvx * y2y3 )
282 !
283 ! HERE TERMS 2 AND 3 HAVE BEEN REMOVED FOR MONOTONICITY
284 ! CROSSED TERMS CANCELLED
285 !
286  xm(ielem,04)=0.d0
287  xm(ielem,05)=0.d0
288  xm(ielem,07)=0.d0
289  xm(ielem,09)=0.d0
290  xm(ielem,10)=0.d0
291  xm(ielem,11)=0.d0
292 !
293 ! TERM 4
294 !
295 ! HERE SOME HORIZONTAL TERMS HAVE BEEN REMOVED BECAUSE
296 ! HORIZONTAL FLUXES THROUGH BOTTOM AND TOP OF PRISM
297 ! ARE NEGLECTED FOR MONOTONICITY
298 ! SEE THE WHOLE TERM IN THE OPTION WITHOUT MONOTONICITY
299 !
300  xm(ielem,03)=-nh1*surfac(ielem)*xs03
301  xm(ielem,08)=-nh2*surfac(ielem)*xs03
302  xm(ielem,12)=-nh3*surfac(ielem)*xs03
303 !
304 !-----------------------------------------------------------------------
305 !
306  ENDDO
307 !
308 !-----------------------------------------------------------------------
309 !
310  ELSEIF(formul(10:13).EQ.'1234') THEN
311 !
312 !-----------------------------------------------------------------------
313 !
314 ! VERSION WITHOUT MONOTONICITY AND ALL TERMS
315 !
316 !-----------------------------------------------------------------------
317 !
318  DO ielem=1,nelem
319 !
320  i1=ikle(ielem,1)
321  i2=ikle(ielem,2)
322  i3=ikle(ielem,3)
323  i4=ikle(ielem,4)
324  i5=ikle(ielem,5)
325  i6=ikle(ielem,6)
326 ! SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
327  ii4=i1+npou0
328  ii5=i2+npou0
329  ii6=i3+npou0
330 !
331  h1=z(i4)-z(i1)
332  h2=z(i5)-z(i2)
333  h3=z(i6)-z(i3)
334 !
335  IF((inchyd.AND.max(z(i1),z(i2),z(i3)).GT.
336  & min(z(i4),z(i5),z(i6))) .OR.
337  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
338  nf1=0.d0
339  nf2=0.d0
340  nf3=0.d0
341  ng1=0.d0
342  ng2=0.d0
343  ng3=0.d0
344  nh1=0.d0
345  nh2=0.d0
346  nh3=0.d0
347  aux=0.d0
348  nuxmoy=0.d0
349  nuymoy=0.d0
350  ELSE
351 ! SUIVANT LES CAS (II4=I1 OU I4, ETC.)
352 ! VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
353 ! VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
354  nf1=0.5d0*(f(i1)+f(ii4))/h1
355  nf2=0.5d0*(f(i2)+f(ii5))/h2
356  nf3=0.5d0*(f(i3)+f(ii6))/h3
357  ng1=0.5d0*(g(i1)+g(ii4))/h1
358  ng2=0.5d0*(g(i2)+g(ii5))/h2
359  ng3=0.5d0*(g(i3)+g(ii6))/h3
360  nh1=0.5d0*(h(i1)+h(ii4))/h1
361  nh2=0.5d0*(h(i2)+h(ii5))/h2
362  nh3=0.5d0*(h(i3)+h(ii6))/h3
363  aux = xs24 / surfac(ielem)
364 ! TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
365  nuxmoy=(f(i1)+f(i2)+f(i3)+f(i4)+f(i5)+f(i6))/6.d0
366  nuymoy=(g(i1)+g(i2)+g(i3)+g(i4)+g(i5)+g(i6))/6.d0
367  ENDIF
368 !
369 ! X2 = X(I2)-X(I1)
370 ! X3 = X(I3)-X(I1)
371 ! Y2 = Y(I2)-Y(I1)
372 ! Y3 = Y(I3)-Y(I1)
373  x2 = x(ielem,2)
374  x3 = x(ielem,3)
375  y2 = y(ielem,2)
376  y3 = y(ielem,3)
377 !
378  x2x3 = x2*x3
379  y2y3 = y2*y3
380  x2aux = x2x3-x2**2
381  x3aux = x3**2-x2x3
382  y2aux = y2y3-y2**2
383  y3aux = y3**2-y2y3
384 !
385 ! 2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
386 ! LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
387  somvx = ( f(i1)*h1+f(i2)*h2+f(i3)*h3 ) * aux
388  somvy = ( g(i1)*h1+g(i2)*h2+g(i3)*h3 ) * aux
389 !
390 ! OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
391 ! WITH MONOTONICITY ENSURED
392 !
393 ! TERM 1
394 !
395  xm(ielem, 1) = - somvy * x3aux - somvx * y3aux
396  xm(ielem, 2) = somvy * x2aux + somvx * y2aux
397  xm(ielem, 6) = - somvy * x2x3 - somvx * y2y3
398 !
399 ! 2D DIFFUSION, UPPER LEVEL
400  somvx = ( f(i4)*h1+f(i5)*h2+f(i6)*h3 ) * aux
401  somvy = ( g(i4)*h1+g(i5)*h2+g(i6)*h3 ) * aux
402 !
403 ! OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
404 ! WITH MONOTONICITY ENSURED
405 !
406 ! TERM 1
407 !
408  xm(ielem,13) = - somvy * x3aux - somvx * y3aux
409  xm(ielem,14) = somvy * x2aux + somvx * y2aux
410  xm(ielem,15) = - somvy * x2x3 - somvx * y2y3
411 !
412 ! AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
413 ! ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
414 ! NPX=-DZ/DX NPY=-DZ/DY
415 !
416  npxl=-0.5d0*(y2*(z(i1)-z(i3))+y3*(z(i2)-z(i1)))
417  npxu=-0.5d0*(y2*(z(i4)-z(i6))+y3*(z(i5)-z(i4)))
418  npyl=-0.5d0*(x2*(z(i3)-z(i1))+x3*(z(i1)-z(i2)))
419  npyu=-0.5d0*(x2*(z(i6)-z(i4))+x3*(z(i4)-z(i5)))
420  npx=0.5d0*(npxl+npxu)/surfac(ielem)
421  npy=0.5d0*(npyl+npyu)/surfac(ielem)
422  npx2=npx**2
423  npy2=npy**2
424 !
425 ! TERMS 2 AND 3
426 !
427  npx=npx*nuxmoy
428  npy=npy*nuymoy
429 !
430 ! 2D GRADIENT MATRIX GX(3,3) AND GY(3,3)
431 ! BUT GX(1,J)=GX(2,J)=GX(3,J)
432 ! AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
433 !
434  gx(2) = y3*xs06
435  gx(3) = - y2*xs06
436  gx(1) = - gx(2) - gx(3)
437 !
438  gy(2) = - x3 * xs06
439  gy(3) = x2 * xs06
440  gy(1) = - gy(2) - gy(3)
441 !
442 ! OFF-DIAGONAL TERMS, SOME (UP AND DOWN) ALREADY INITIALISED,
443 ! SOME (CROSSED TERMS) NOT
444 !
445 ! TERMS 2
446 !
447 ! TERM 1-2 (01)
448  xm(ielem,01)=xm(ielem,01)+0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
449 ! TERM 1-3 (02)
450  xm(ielem,02)=xm(ielem,02)+0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
451 ! TERM 1-4 (03)
452  xm(ielem,03)= +0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
453 ! TERM 1-5 (04)
454  xm(ielem,04)= +0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
455 ! TERM 1-6 (05)
456  xm(ielem,05)= +0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
457 ! TERM 2-3 (06)
458  xm(ielem,06)=xm(ielem,06)+0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
459 ! TERM 2-4 (07)
460  xm(ielem,07)= +0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
461 ! TERM 2-5 (08)
462  xm(ielem,08)= +0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
463 ! TERM 2-6 (09)
464  xm(ielem,09)= +0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
465 ! TERM 3-4 (10)
466  xm(ielem,10)= +0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
467 ! TERM 3-5 (11)
468  xm(ielem,11)= +0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
469 ! TERM 3-6 (12)
470  xm(ielem,12)= +0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
471 ! TERME 4-5 (13)
472  xm(ielem,13)=xm(ielem,13)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
473 ! TERM 4-6 (14)
474  xm(ielem,14)=xm(ielem,14)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
475 ! TERM 5-6 (15)
476  xm(ielem,15)=xm(ielem,15)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
477 !
478 ! TERMS 3
479 !
480 ! TERM 1-2 (01)
481  xm(ielem,01)=xm(ielem,01)+0.5d0*( - ( npx*gx(1)+npy*gy(1)) )
482 ! TERM 1-3 (02)
483  xm(ielem,02)=xm(ielem,02)+0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
484 ! TERM 1-4 (03)
485  xm(ielem,03)=xm(ielem,03)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
486 ! TERM 1-5 (04)
487  xm(ielem,04)=xm(ielem,04)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
488 ! TERM 1-6 (05)
489  xm(ielem,05)=xm(ielem,05)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
490 ! TERM 2-3 (06)
491  xm(ielem,06)=xm(ielem,06)+0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
492 ! TERM 2-4 (07)
493  xm(ielem,07)=xm(ielem,07)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
494 ! TERM 2-5 (08)
495  xm(ielem,08)=xm(ielem,08)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
496 ! TERM 2-6 (09)
497  xm(ielem,09)=xm(ielem,09)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
498 ! TERM 3-4 (10)
499  xm(ielem,10)=xm(ielem,10)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
500 ! TERM 3-5 (11)
501  xm(ielem,11)=xm(ielem,11)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
502 ! TERM 3-6 (12)
503  xm(ielem,12)=xm(ielem,12)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
504 ! TERME 4-5 (13)
505  xm(ielem,13)=xm(ielem,13)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
506 ! TERM 4-6 (14)
507  xm(ielem,14)=xm(ielem,14)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
508 ! TERM 5-6 (15)
509  xm(ielem,15)=xm(ielem,15)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
510 !
511 ! TERM 4
512 !
513  xm(ielem,03)=xm(ielem,03)
514  & -(npx2*nf1+npy2*ng1+nh1)*surfac(ielem)*xs03
515  xm(ielem,08)=xm(ielem,08)
516  & -(npx2*nf2+npy2*ng2+nh2)*surfac(ielem)*xs03
517  xm(ielem,12)=xm(ielem,12)
518  & -(npx2*nf3+npy2*ng3+nh3)*surfac(ielem)*xs03
519 !
520 !-----------------------------------------------------------------------
521 !
522  ENDDO
523 !
524 !-----------------------------------------------------------------------
525 !
526  ELSEIF(formul(10:13).EQ.'1 3 ') THEN
527 !
528 !-----------------------------------------------------------------------
529 !
530 ! VERSION WITHOUT MONOTONICITY AND ONLY TERMS 1 AND 3
531 !
532 !-----------------------------------------------------------------------
533 !
534  DO ielem=1,nelem
535 !
536  i1=ikle(ielem,1)
537  i2=ikle(ielem,2)
538  i3=ikle(ielem,3)
539  i4=ikle(ielem,4)
540  i5=ikle(ielem,5)
541  i6=ikle(ielem,6)
542 ! SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
543  ii4=i1+npou0
544  ii5=i2+npou0
545  ii6=i3+npou0
546 !
547  h1=z(i4)-z(i1)
548  h2=z(i5)-z(i2)
549  h3=z(i6)-z(i3)
550 !
551  IF((inchyd.AND.max(z(i1),z(i2),z(i3)).GT.
552  & min(z(i4),z(i5),z(i6))) .OR.
553  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
554  nf1=0.d0
555  nf2=0.d0
556  nf3=0.d0
557  ng1=0.d0
558  ng2=0.d0
559  ng3=0.d0
560  nh1=0.d0
561  nh2=0.d0
562  nh3=0.d0
563  aux=0.d0
564  nuxmoy=0.d0
565  nuymoy=0.d0
566  ELSE
567 ! SUIVANT LES CAS (II4=I1 OU I4, ETC.)
568 ! VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
569 ! VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
570  nf1=0.5d0*(f(i1)+f(ii4))/h1
571  nf2=0.5d0*(f(i2)+f(ii5))/h2
572  nf3=0.5d0*(f(i3)+f(ii6))/h3
573  ng1=0.5d0*(g(i1)+g(ii4))/h1
574  ng2=0.5d0*(g(i2)+g(ii5))/h2
575  ng3=0.5d0*(g(i3)+g(ii6))/h3
576  nh1=0.5d0*(h(i1)+h(ii4))/h1
577  nh2=0.5d0*(h(i2)+h(ii5))/h2
578  nh3=0.5d0*(h(i3)+h(ii6))/h3
579  aux = xs24 / surfac(ielem)
580 ! TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
581  nuxmoy=(f(i1)+f(i2)+f(i3)+f(i4)+f(i5)+f(i6))/6.d0
582  nuymoy=(g(i1)+g(i2)+g(i3)+g(i4)+g(i5)+g(i6))/6.d0
583  ENDIF
584 !
585 ! X2 = X(I2)-X(I1)
586 ! X3 = X(I3)-X(I1)
587 ! Y2 = Y(I2)-Y(I1)
588 ! Y3 = Y(I3)-Y(I1)
589  x2 = x(ielem,2)
590  x3 = x(ielem,3)
591  y2 = y(ielem,2)
592  y3 = y(ielem,3)
593 !
594  x2x3 = x2*x3
595  y2y3 = y2*y3
596  x2aux = x2x3-x2**2
597  x3aux = x3**2-x2x3
598  y2aux = y2y3-y2**2
599  y3aux = y3**2-y2y3
600 !
601 ! 2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
602 ! LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
603  somvx = ( f(i1)*h1+f(i2)*h2+f(i3)*h3 ) * aux
604  somvy = ( g(i1)*h1+g(i2)*h2+g(i3)*h3 ) * aux
605 !
606 ! OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
607 ! WITH MONOTONICITY ENSURED
608 !
609 ! TERM 1
610 !
611  xm(ielem, 1) = - somvy * x3aux - somvx * y3aux
612  xm(ielem, 2) = somvy * x2aux + somvx * y2aux
613  xm(ielem, 6) = - somvy * x2x3 - somvx * y2y3
614  xm(ielem,16) = xm(ielem, 1)
615  xm(ielem,17) = xm(ielem, 2)
616  xm(ielem,21) = xm(ielem, 6)
617 !
618 ! 2D DIFFUSION, UPPER LEVEL
619  somvx = ( f(i4)*h1+f(i5)*h2+f(i6)*h3 ) * aux
620  somvy = ( g(i4)*h1+g(i5)*h2+g(i6)*h3 ) * aux
621 !
622 ! OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
623 ! WITH MONOTONICITY ENSURED
624 !
625 ! TERM 1
626 !
627  xm(ielem,13) = - somvy * x3aux - somvx * y3aux
628  xm(ielem,14) = somvy * x2aux + somvx * y2aux
629  xm(ielem,15) = - somvy * x2x3 - somvx * y2y3
630  xm(ielem,28) = xm(ielem,13)
631  xm(ielem,29) = xm(ielem,14)
632  xm(ielem,30) = xm(ielem,15)
633 !
634 ! AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
635 ! ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
636 ! NPX=-DZ/DX NPY=-DZ/DY
637 !
638  npxl=-0.5d0*(y2*(z(i1)-z(i3))+y3*(z(i2)-z(i1)))
639  npxu=-0.5d0*(y2*(z(i4)-z(i6))+y3*(z(i5)-z(i4)))
640  npyl=-0.5d0*(x2*(z(i3)-z(i1))+x3*(z(i1)-z(i2)))
641  npyu=-0.5d0*(x2*(z(i6)-z(i4))+x3*(z(i4)-z(i5)))
642  npx=0.5d0*(npxl+npxu)/surfac(ielem)
643  npy=0.5d0*(npyl+npyu)/surfac(ielem)
644  npx2=npx**2
645  npy2=npy**2
646 !
647 ! TERM 3
648 !
649  npx=npx*nuxmoy
650  npy=npy*nuymoy
651 !
652 ! 2D GRADIENT MATRIX GX(3,3) AND GY(3,3)
653 ! BUT GX(1,J)=GX(2,J)=GX(3,J)
654 ! AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
655 !
656  gx(2) = y3*xs06
657  gx(3) = - y2*xs06
658  gx(1) = - gx(2) - gx(3)
659 !
660  gy(2) = - x3 * xs06
661  gy(3) = x2 * xs06
662  gy(1) = - gy(2) - gy(3)
663 !
664 ! OFF-DIAGONAL TERMS, SOME (UP AND DOWN) ALREADY INITIALISED,
665 ! SOME (CROSSED TERMS) NOT
666 !
667 ! TERMS 3
668 !
669 ! TERM 1-2 (01)
670  xm(ielem,01)=xm(ielem,01)+0.5d0*( - ( npx*gx(1)+npy*gy(1)) )
671 ! TERM 1-3 (02)
672  xm(ielem,02)=xm(ielem,02)+0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
673 ! TERM 1-4 (03)
674  xm(ielem,03)= +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
675 ! TERM 1-5 (04)
676  xm(ielem,04)= +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
677 ! TERM 1-6 (05)
678  xm(ielem,05)= +0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
679 ! TERM 2-3 (06)
680  xm(ielem,06)=xm(ielem,06)+0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
681 ! TERM 2-4 (07)
682  xm(ielem,07)= +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
683 ! TERM 2-5 (08)
684  xm(ielem,08)= +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
685 ! TERM 2-6 (09)
686  xm(ielem,09)= +0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
687 ! TERM 3-4 (10)
688  xm(ielem,10)= +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
689 ! TERM 3-5 (11)
690  xm(ielem,11)= +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
691 ! TERM 3-6 (12)
692  xm(ielem,12)= +0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
693 ! TERME 4-5 (13)
694  xm(ielem,13)=xm(ielem,13)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
695 ! TERM 4-6 (14)
696  xm(ielem,14)=xm(ielem,14)+0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
697 ! TERM 5-6 (15)
698  xm(ielem,15)=xm(ielem,15)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
699 !
700 ! TERM 2-1 (16)
701  xm(ielem,16)=xm(ielem,16)+0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
702 ! TERM 3-1 (17)
703  xm(ielem,17)=xm(ielem,17)+0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
704 ! TERM 4-1 (18)
705  xm(ielem,18)= 0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
706 ! TERM 5-1 (19)
707  xm(ielem,19)= 0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
708 ! TERM 6-1 (20)
709  xm(ielem,20)= 0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
710 ! TERM 3-2 (21)
711  xm(ielem,21)=xm(ielem,21)+0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
712 ! TERM 4-2 (22)
713  xm(ielem,22)= 0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
714 ! TERM 5-2 (23)
715  xm(ielem,23)= 0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
716 ! TERM 6-2 (24)
717  xm(ielem,24)= 0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
718 ! TERM 4-3 (25)
719  xm(ielem,25)= 0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
720 ! TERM 5-3 (26)
721  xm(ielem,26)= 0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
722 ! TERM 5-4 (28)
723  xm(ielem,28)=xm(ielem,28)+0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
724 ! TERM 6-3 (27)
725  xm(ielem,27)= 0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
726 ! TERM 6-4 (29)
727  xm(ielem,29)=xm(ielem,29)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
728 ! TERM 6-5 (30)
729  xm(ielem,30)=xm(ielem,30)+0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
730 !
731 !-----------------------------------------------------------------------
732 !
733  ENDDO
734 !
735 !-----------------------------------------------------------------------
736 !
737  ELSEIF(formul(10:13).EQ.' 2 4') THEN
738 !
739 !-----------------------------------------------------------------------
740 !
741 ! VERSION WITHOUT MONOTONICITY AND ONLY TERMS 2 AND 4
742 !
743 !-----------------------------------------------------------------------
744 !
745  DO ielem=1,nelem
746 !
747  i1=ikle(ielem,1)
748  i2=ikle(ielem,2)
749  i3=ikle(ielem,3)
750  i4=ikle(ielem,4)
751  i5=ikle(ielem,5)
752  i6=ikle(ielem,6)
753 ! SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
754  ii4=i1+npou0
755  ii5=i2+npou0
756  ii6=i3+npou0
757 !
758  h1=z(i4)-z(i1)
759  h2=z(i5)-z(i2)
760  h3=z(i6)-z(i3)
761 !
762  IF((inchyd.AND.max(z(i1),z(i2),z(i3)).GT.
763  & min(z(i4),z(i5),z(i6))) .OR.
764  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
765  nf1=0.d0
766  nf2=0.d0
767  nf3=0.d0
768  ng1=0.d0
769  ng2=0.d0
770  ng3=0.d0
771  nh1=0.d0
772  nh2=0.d0
773  nh3=0.d0
774  nuxmoy=0.d0
775  nuymoy=0.d0
776  ELSE
777 ! SUIVANT LES CAS (II4=I1 OU I4, ETC.)
778 ! VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
779 ! VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
780  nf1=0.5d0*(f(i1)+f(ii4))/h1
781  nf2=0.5d0*(f(i2)+f(ii5))/h2
782  nf3=0.5d0*(f(i3)+f(ii6))/h3
783  ng1=0.5d0*(g(i1)+g(ii4))/h1
784  ng2=0.5d0*(g(i2)+g(ii5))/h2
785  ng3=0.5d0*(g(i3)+g(ii6))/h3
786  nh1=0.5d0*(h(i1)+h(ii4))/h1
787  nh2=0.5d0*(h(i2)+h(ii5))/h2
788  nh3=0.5d0*(h(i3)+h(ii6))/h3
789 ! TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
790  nuxmoy=(f(i1)+f(i2)+f(i3)+f(i4)+f(i5)+f(i6))/6.d0
791  nuymoy=(g(i1)+g(i2)+g(i3)+g(i4)+g(i5)+g(i6))/6.d0
792  ENDIF
793 !
794 ! X2 = X(I2)-X(I1)
795 ! X3 = X(I3)-X(I1)
796 ! Y2 = Y(I2)-Y(I1)
797 ! Y3 = Y(I3)-Y(I1)
798  x2 = x(ielem,2)
799  x3 = x(ielem,3)
800  y2 = y(ielem,2)
801  y3 = y(ielem,3)
802 !
803 ! AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
804 ! ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
805 ! NPX=-DZ/DX NPY=-DZ/DY
806 !
807  npxl=-0.5d0*(y2*(z(i1)-z(i3))+y3*(z(i2)-z(i1)))
808  npxu=-0.5d0*(y2*(z(i4)-z(i6))+y3*(z(i5)-z(i4)))
809  npyl=-0.5d0*(x2*(z(i3)-z(i1))+x3*(z(i1)-z(i2)))
810  npyu=-0.5d0*(x2*(z(i6)-z(i4))+x3*(z(i4)-z(i5)))
811  npx=0.5d0*(npxl+npxu)/surfac(ielem)
812  npy=0.5d0*(npyl+npyu)/surfac(ielem)
813  npx2=npx**2
814  npy2=npy**2
815 !
816 ! TERMS 2
817 !
818  npx=npx*nuxmoy
819  npy=npy*nuymoy
820 !
821 ! 2D GRADIENT MATRIX GX(3,3) AND GY(3,3)
822 ! BUT GX(1,J)=GX(2,J)=GX(3,J)
823 ! AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
824 !
825  gx(2) = y3*xs06
826  gx(3) = - y2*xs06
827  gx(1) = - gx(2) - gx(3)
828 !
829  gy(2) = - x3 * xs06
830  gy(3) = x2 * xs06
831  gy(1) = - gy(2) - gy(3)
832 !
833 ! OFF-DIAGONAL TERMS
834 !
835 ! TERMS 2
836 !
837 ! TERM 1-2 (01)
838  xm(ielem,01)=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
839 ! TERM 1-3 (02)
840  xm(ielem,02)=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
841 ! TERM 1-4 (03)
842  xm(ielem,03)=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
843 ! TERM 1-5 (04)
844  xm(ielem,04)=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
845 ! TERM 1-6 (05)
846  xm(ielem,05)=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
847 ! TERM 2-3 (06)
848  xm(ielem,06)=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
849 ! TERM 2-4 (07)
850  xm(ielem,07)=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
851 ! TERM 2-5 (08)
852  xm(ielem,08)=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
853 ! TERM 2-6 (09)
854  xm(ielem,09)=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
855 ! TERM 3-4 (10)
856  xm(ielem,10)=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
857 ! TERM 3-5 (11)
858  xm(ielem,11)=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
859 ! TERM 3-6 (12)
860  xm(ielem,12)=0.5d0*( - ( npx*gx(3)+npy*gy(3) ) )
861 ! TERM 4-5 (13)
862  xm(ielem,13)=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
863 ! TERM 4-6 (14)
864  xm(ielem,14)=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
865 ! TERM 5-6 (15)
866  xm(ielem,15)=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
867 !
868 ! TERM 2-1 (16)
869  xm(ielem,16)=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
870 ! TERM 3-1 (17)
871  xm(ielem,17)=0.5d0*( - ( npx*gx(1)+npy*gy(1) ) )
872 ! TERM 4-1 (18)
873  xm(ielem,18)=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
874 ! TERM 5-1 (19)
875  xm(ielem,19)=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
876 ! TERM 6-1 (20)
877  xm(ielem,20)=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
878 ! TERM 3-2 (21)
879  xm(ielem,21)=0.5d0*( - ( npx*gx(2)+npy*gy(2) ) )
880 ! TERM 4-2 (22)
881  xm(ielem,22)=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
882 ! TERM 5-2 (23)
883  xm(ielem,23)=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
884 ! TERM 6-2 (24)
885  xm(ielem,24)=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
886 ! TERM 4-3 (25)
887  xm(ielem,25)=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
888 ! TERM 5-3 (26)
889  xm(ielem,26)=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
890 ! TERM 6-3 (27)
891  xm(ielem,27)=0.5d0*( + ( npx*gx(3)+npy*gy(3) ) )
892 ! TERM 5-4 (28)
893  xm(ielem,28)=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
894 ! TERM 6-4 (29)
895  xm(ielem,29)=0.5d0*( + ( npx*gx(1)+npy*gy(1) ) )
896 ! TERM 6-5 (30)
897  xm(ielem,30)=0.5d0*( + ( npx*gx(2)+npy*gy(2) ) )
898 !
899 ! TERM 4
900 !
901  xm(ielem,03)=xm(ielem,03)
902  & -(npx2*nf1+npy2*ng1+nh1)*surfac(ielem)*xs03
903  xm(ielem,08)=xm(ielem,08)
904  & -(npx2*nf2+npy2*ng2+nh2)*surfac(ielem)*xs03
905  xm(ielem,12)=xm(ielem,12)
906  & -(npx2*nf3+npy2*ng3+nh3)*surfac(ielem)*xs03
907  xm(ielem,18)=xm(ielem,18)
908  & -(npx2*nf1+npy2*ng1+nh1)*surfac(ielem)*xs03
909  xm(ielem,23)=xm(ielem,23)
910  & -(npx2*nf2+npy2*ng2+nh2)*surfac(ielem)*xs03
911  xm(ielem,27)=xm(ielem,27)
912  & -(npx2*nf3+npy2*ng3+nh3)*surfac(ielem)*xs03
913 !
914 !-----------------------------------------------------------------------
915 !
916  ENDDO
917 !
918 !-----------------------------------------------------------------------
919 !
920  ELSE
921  WRITE(lu,*) 'MT02PP_STAR (BIEF): UNKNOWN FORMULA'
922  WRITE(lu,*) formul
923  CALL plante(1)
924  stop
925  ENDIF
926 !
927 !-----------------------------------------------------------------------
928 !
929 ! DIAGONAL TERMS OBTAINED BY THE FACT THAT THE SUM OF TERMS IN A
930 ! LINE IS 0
931 !
932 !-----------------------------------------------------------------------
933 !
934 ! IN SYMMETRIC MODE
935 !
936  IF(formul(11:13).EQ.'MON'.OR.formul(10:13).EQ.'1234') THEN
937 !
938  DO ielem=1,nelem
939 !
940  t(ielem,1)= -xm(ielem,01)
941  & -xm(ielem,02)
942  & -xm(ielem,03)
943  & -xm(ielem,04)
944  & -xm(ielem,05)
945  t(ielem,2)= -xm(ielem,01)
946  & -xm(ielem,06)
947  & -xm(ielem,07)
948  & -xm(ielem,08)
949  & -xm(ielem,09)
950  t(ielem,3)= -xm(ielem,02)
951  & -xm(ielem,06)
952  & -xm(ielem,10)
953  & -xm(ielem,11)
954  & -xm(ielem,12)
955  t(ielem,4)= -xm(ielem,03)
956  & -xm(ielem,07)
957  & -xm(ielem,10)
958  & -xm(ielem,13)
959  & -xm(ielem,14)
960  t(ielem,5)= -xm(ielem,04)
961  & -xm(ielem,08)
962  & -xm(ielem,11)
963  & -xm(ielem,13)
964  & -xm(ielem,15)
965  t(ielem,6)= -xm(ielem,05)
966  & -xm(ielem,09)
967  & -xm(ielem,12)
968  & -xm(ielem,14)
969  & -xm(ielem,15)
970 !
971  ENDDO
972 !
973 ! IN NON SYMMETRIC MODE
974 !
975  ELSEIF(formul(10:13).EQ.'1 3 '.OR.formul(10:13).EQ.' 2 4') THEN
976 !
977  DO ielem=1,nelem
978 !
979  t(ielem,1)= -xm(ielem,01)
980  & -xm(ielem,02)
981  & -xm(ielem,03)
982  & -xm(ielem,04)
983  & -xm(ielem,05)
984  t(ielem,2)= -xm(ielem,16)
985  & -xm(ielem,06)
986  & -xm(ielem,07)
987  & -xm(ielem,08)
988  & -xm(ielem,09)
989  t(ielem,3)= -xm(ielem,17)
990  & -xm(ielem,21)
991  & -xm(ielem,10)
992  & -xm(ielem,11)
993  & -xm(ielem,12)
994  t(ielem,4)= -xm(ielem,18)
995  & -xm(ielem,22)
996  & -xm(ielem,25)
997  & -xm(ielem,13)
998  & -xm(ielem,14)
999  t(ielem,5)= -xm(ielem,19)
1000  & -xm(ielem,23)
1001  & -xm(ielem,26)
1002  & -xm(ielem,28)
1003  & -xm(ielem,15)
1004  t(ielem,6)= -xm(ielem,20)
1005  & -xm(ielem,24)
1006  & -xm(ielem,27)
1007  & -xm(ielem,29)
1008  & -xm(ielem,30)
1009 !
1010  ENDDO
1011 !
1012  ELSE
1013  WRITE(lu,*) 'MT02PP_STAR (BIEF): UNKNOWN FORMULA'
1014  WRITE(lu,*) formul
1015  CALL plante(1)
1016  stop
1017  ENDIF
1018 !
1019 !-----------------------------------------------------------------------
1020 !
1021  RETURN
1022  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt02pp_star
subroutine mt02pp_star(T, XM, XMUL, SF, SG, SH, F, G, H, X, Y, Z, SURFAC, IKLE, NELEM, NELMAX, INCHYD, FORMUL, NPLAN)
Definition: mt02pp_star.f:8
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:46 for The TELEMAC-MASCARET system by doxygen 1.8.13