The TELEMAC-MASCARET system  trunk
mt02pp.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE mt02pp
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 V6P2 21/08/2010
10 !***********************************************************************
11 !
12 !brief COMPUTES THE DIFFUSION MATRIX.
13 !+
14 !+ THE FUNCTION DIFFUSION COEFFICIENT IS HERE A P1
15 !+ DIAGONAL TENSOR.
16 !+
17 !+ CASE OF THE PRISM.
18 !code
19 !+ STORAGE CONVENTION FOR EXTRA-DIAGONAL TERMS:
20 !+
21 !+ XM(IELEM, 1) ----> M(1,2) = M(2,1)
22 !+ XM(IELEM, 2) ----> M(1,3) = M(3,1)
23 !+ XM(IELEM, 3) ----> M(1,4) = M(4,1)
24 !+ XM(IELEM, 4) ----> M(1,5) = M(5,1)
25 !+ XM(IELEM, 5) ----> M(1,6) = M(6,1)
26 !+ XM(IELEM, 6) ----> M(2,3) = M(3,2)
27 !+ XM(IELEM, 7) ----> M(2,4) = M(4,2)
28 !+ XM(IELEM, 8) ----> M(2,5) = M(5,2)
29 !+ XM(IELEM, 9) ----> M(2,6) = M(6,2)
30 !+ XM(IELEM,10) ----> M(3,4) = M(4,3)
31 !+ XM(IELEM,11) ----> M(3,5) = M(5,3)
32 !+ XM(IELEM,12) ----> M(3,6) = M(6,3)
33 !+ XM(IELEM,13) ----> M(4,5) = M(5,4)
34 !+ XM(IELEM,14) ----> M(4,6) = M(6,4)
35 !+ XM(IELEM,15) ----> M(5,6) = M(6,5)
36 !
37 !warning THE JACOBIAN MUST BE POSITIVE
38 !warning JMH : THIS IS NOT DONE AS IN 2D, TO BE MODIFIED
39 !
40 !history
41 !+ 22/08/07
42 !+
43 !+ CAN OPTIONALLY TREAT (USES DIMDISC) A P0 VERTICAL
44 !
45 !history
46 !+ 06/02/07
47 !+
48 !+ IF FORMUL(14:16)='MON' :
49 !
50 !history JMH
51 !+ 15/03/2010
52 !+
53 !+ PARAMETER EPSILON ADDED
54 !
55 !history J-M HERVOUET (LNHE) ; F LEPEINTRE (LNH)
56 !+ 20/05/2010
57 !+ V6P0
58 !+
59 !
60 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
61 !+ 13/07/2010
62 !+ V6P0
63 !+ Translation of French comments within the FORTRAN sources into
64 !+ English comments
65 !
66 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
67 !+ 21/08/2010
68 !+ V6P0
69 !+ Creation of DOXYGEN tags for automated documentation and
70 !+ cross-referencing of the FORTRAN sources
71 !
72 !history U.H.MErkel
73 !+ 18/07/2012
74 !+ V6P2
75 !+ Replaced EPSILON with CHOUIA due to nag compiler problems
76 !
77 !history J-M HERVOUET (EDF LAB, LNHE)
78 !+ 11/03/2014
79 !+ V7P0
80 !+ Now SH%TYPR checked to cancel vertical diffusion.
81 !
82 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
83 !| F |-->| FUNCTION USED IN THE FORMULA
84 !| FORMUL |-->| FORMULA DESCRIBING THE RESULTING MATRIX
85 !| G |-->| FUNCTION USED IN THE FORMULA
86 !| H |-->| FUNCTION USED IN THE FORMULA
87 !| IKLE |-->| CONNECTIVITY TABLE.
88 !| INCHYD |---| IF YES, TREATS HYDROSTATIC INCONSISTENCIES
89 !| NELEM |-->| NUMBER OF ELEMENTS
90 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
91 !| NPLAN |-->| NUMBER OF PLANES IN THE MESH OF PRISMS
92 !| SF |-->| STRUCTURE OF FUNCTIONS F
93 !| SG |-->| STRUCTURE OF FUNCTIONS G
94 !| SH |-->| STRUCTURE OF FUNCTIONS H
95 !| SURFAC |-->| AREA OF 2D ELEMENTS
96 !| T |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
97 !| X |-->| ABSCISSAE OF POINTS, PER ELEMENT
98 !| Y |-->| ORDINATES OF POINTS, PER ELEMENT
99 !| Z |-->| ELEVATIONS OF POINTS, PER POINT
100 !| XM |<->| OFF-DIAGONAL TERMS
101 !| XMUL |-->| COEFFICIENT FOR MULTIPLICATION
102 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
103 !
104  USE bief, ex_mt02pp => mt02pp
105 !
106  USE declarations_special
107  IMPLICIT NONE
108 !
109 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
110 !
111  INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN
112  INTEGER, INTENT(IN) :: IKLE(nelmax,6)
113 !
114  DOUBLE PRECISION, INTENT(INOUT) :: T(nelmax,6),XM(nelmax,15)
115  DOUBLE PRECISION, INTENT(IN) :: SURFAC(nelmax)
116 !
117  DOUBLE PRECISION, INTENT(IN) :: XMUL
118  DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),H(*)
119 !
120 ! STRUCTURES OF F, G, H
121 !
122  TYPE(bief_obj), INTENT(IN) :: SF,SG,SH
123 !
124  DOUBLE PRECISION, INTENT(IN) :: X(nelmax,6),Y(nelmax,6),Z(*)
125 !
126  LOGICAL, INTENT(IN) :: INCHYD
127  CHARACTER(LEN=16), INTENT(IN) :: FORMUL
128 !
129 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
130 !
131 ! DECLARATIONS SPECIFIC TO THIS SUBROUTINE
132 !
133  DOUBLE PRECISION H1,H2,H3,RI,RS,R,RRI,RRS,RRRI,RRRS
134  DOUBLE PRECISION D1,D2,D3,D12,D13,D23,D
135  DOUBLE PRECISION SH1,SHH,SNI,SNS,SNHI,SNHS,SNHHI,SNHHS,SNHH
136  DOUBLE PRECISION SNHI1,SNHS1,SNHI2,SNHS2,SNHI3,SNHS3
137  DOUBLE PRECISION HHI12,HHS12,HH12,HHI13,HHS13,HH13
138  DOUBLE PRECISION HHI23,HHS23,HH23
139  DOUBLE PRECISION HRI1,HRS1,HRI2,HRS2,HRI3,HRS3
140  DOUBLE PRECISION RR11,RR22,RR33,RR12,RR13,RR23
141  DOUBLE PRECISION NF1,NF2,NF3,NF4,NF5,NF6
142  DOUBLE PRECISION NG1,NG2,NG3,NG4,NG5,NG6
143  DOUBLE PRECISION NH1,NH2,NH3,XS06,XS2880
144 !
145  INTEGER I1,I2,I3,I4,I5,I6,IELEM,II4,II5,II6,NPOU0
146 !
147  DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-4
148 !
149 !-----------------------------------------------------------------------
150 !
151  xs06=xmul/6.d0
152  xs2880=xmul/2880.d0
153 !
154 !-----------------------------------------------------------------------
155 !
156 ! VERY IMPORTANT !!!!!!!!!
157 !
158 ! SH%TYPR CHECKED TO CANCEL DIFFUSION ALONG Z
159 !
160 ! NOTE: XS06 ONLY USED WITH DIFFUSION ALONG Z !!!!!!!!!
161 !
162  IF(sh%TYPR.EQ.'0') THEN
163  xs06=0.d0
164  ENDIF
165 !
166 ! COULD BE OPTIMISED MORE BUT WOULD REQUEST SPLITTING LOOPS
167 ! INTO HORIZONTAL AND VERTICAL
168 !
169 !-----------------------------------------------------------------------
170 !
171  IF(sf%ELM.NE.41) THEN
172  WRITE(lu,1001) sf%ELM
173 1001 FORMAT(1x,'MT02PP (BIEF) : TYPE OF F NOT IMPLEMENTED: ',i6)
174  CALL plante(1)
175  stop
176  ENDIF
177  IF(sg%ELM.NE.41) THEN
178  WRITE(lu,2001) sg%ELM
179 2001 FORMAT(1x,'MT02PP (BIEF) : TYPE OF G NOT IMPLEMENTED: ',i6)
180  CALL plante(1)
181  stop
182  ENDIF
183  IF(sh%ELM.NE.41) THEN
184  WRITE(lu,3001) sh%ELM
185 3001 FORMAT(1x,'MT02PP (BIEF) : TYPE OF H NOT IMPLEMENTED: ',i6)
186  CALL plante(1)
187  stop
188  ENDIF
189 !
190 ! SEE VISCLM OF TELEMAC-3D
191 ! FOR THE TREATMENT OF P0 VERTICAL VISCOSITY ON THE VERTICAL
192 !
193  IF(sh%DIMDISC.EQ.0) THEN
194 ! P1 VERTICAL VISCOSITY
195  npou0=sh%DIM1/nplan
196  ELSEIF(sh%DIMDISC.EQ.4111) THEN
197 ! P0 VERTICAL VISCOSITY ON THE VERTICAL (SEE II4,5,6)
198  npou0=0
199  ELSE
200  WRITE(lu,4001) sh%DIMDISC
201 4001 FORMAT(1x,'MT02PP (BIEF): DIMDISC OF H NOT IMPLEMENTED: ',i6)
202  CALL plante(1)
203  stop
204  ENDIF
205 !
206 ! VERSION WITH TREATMENT OF HYDROSTATIC INCONSISTENCIES
207 !
208  IF(inchyd) THEN
209 !
210 !-----------------------------------------------------------------------
211 !
212 ! DIFFUSION ALONG X
213 !
214 !-----------------------------------------------------------------------
215 !
216 ! LOOP ON THE ELEMENTS
217 !
218  DO ielem=1,nelem
219 !
220  i1=ikle(ielem,1)
221  i2=ikle(ielem,2)
222  i3=ikle(ielem,3)
223  i4=ikle(ielem,4)
224  i5=ikle(ielem,5)
225  i6=ikle(ielem,6)
226 ! DEPENDING ON NPOU0, II4 WILL BE I4 OR I1, ETC
227  ii4=i1+npou0
228  ii5=i2+npou0
229  ii6=i3+npou0
230 !
231  h1=z(i4)-z(i1)
232  h2=z(i5)-z(i2)
233  h3=z(i6)-z(i3)
234 !
235  sh1=h1+h2+h3
236  shh=h1*h1+h2*h2+h3*h3
237 !
238  d1=y(ielem,2)-y(ielem,3)
239 ! D2=Y(IELEM,3)-Y(IELEM,1)
240  d2=y(ielem,3)
241 ! D3=Y(IELEM,1)-Y(IELEM,2)
242  d3= -y(ielem,2)
243 !
244  d=xs2880/surfac(ielem)
245 !
246  ri=-(z(i1)*d1+z(i2)*d2+z(i3)*d3)
247  rs=-(z(i4)*d1+z(i5)*d2+z(i6)*d3)
248  r=ri+rs
249  rri=r+ri+ri
250  rrs=r+rs+rs
251  rrri=r*r+2*ri*ri
252  rrrs=r*r+2*rs*rs
253 !
254  d12=d1*d2
255  d13=d1*d3
256  d23=d2*d3
257 !
258  IF(max(z(i1),z(i2),z(i3)).GT.min(z(i4),z(i5),z(i6)).OR.
259  & h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia ) THEN
260  nf1=0.d0
261  nf2=0.d0
262  nf3=0.d0
263  nf4=0.d0
264  nf5=0.d0
265  nf6=0.d0
266  ng1=0.d0
267  ng2=0.d0
268  ng3=0.d0
269  ng4=0.d0
270  ng5=0.d0
271  ng6=0.d0
272  nh1=0.d0
273  nh2=0.d0
274  nh3=0.d0
275  ELSE
276  nf1=f(i1)/h1
277  nf2=f(i2)/h2
278  nf3=f(i3)/h3
279  nf4=f(i4)/h1
280  nf5=f(i5)/h2
281  nf6=f(i6)/h3
282  ng1=g(i1)/h1
283  ng2=g(i2)/h2
284  ng3=g(i3)/h3
285  ng4=g(i4)/h1
286  ng5=g(i5)/h2
287  ng6=g(i6)/h3
288 ! DEPENDING ON THE CASE (II4=I1 OR I4, ETC.)
289 ! ALTERNATIVE WITH VERTICAL LINEAR VISCOSITY (II4=I4)
290 ! ALTERNATIVE WITH P0 VERTICAL VISCOSITY ON THE VERTICAL (II4=I1)
291  nh1=(h(i1)+h(ii4))/h1
292  nh2=(h(i2)+h(ii5))/h2
293  nh3=(h(i3)+h(ii6))/h3
294  ENDIF
295 !
296  sni=nf1+nf2+nf3
297  sns=nf4+nf5+nf6
298  snhi=nf1*h1+nf2*h2+nf3*h3
299  snhs=nf4*h1+nf5*h2+nf6*h3
300  snhhi=(sni*sh1+snhi+snhi)*sh1+sni*shh
301  & +2*(nf1*h1*h1+nf2*h2*h2+nf3*h3*h3)
302  snhhs=(sns*sh1+snhs+snhs)*sh1+sns*shh
303  & +2*(nf4*h1*h1+nf5*h2*h2+nf6*h3*h3)
304  snhh=snhhi+snhhs
305  snhhi=snhh+snhhi+snhhi
306  snhhs=snhh+snhhs+snhhs
307 !
308  hhi12=snhhi*d12
309  hhi13=snhhi*d13
310  hhi23=snhhi*d23
311  hhs12=snhhs*d12
312  hhs13=snhhs*d13
313  hhs23=snhhs*d23
314  hh12=snhh*d12
315  hh13=snhh*d13
316  hh23=snhh*d23
317 !
318  snhi1=(sni+nf1)*(sh1+h1)+snhi+nf1*h1
319  snhs1=(sns+nf4)*(sh1+h1)+snhs+nf4*h1
320  snhi2=(sni+nf2)*(sh1+h2)+snhi+nf2*h2
321  snhs2=(sns+nf5)*(sh1+h2)+snhs+nf5*h2
322  snhi3=(sni+nf3)*(sh1+h3)+snhi+nf3*h3
323  snhs3=(sns+nf6)*(sh1+h3)+snhs+nf6*h3
324 !
325  hri1=rri*snhi1+r*snhs1
326  hrs1=rrs*snhs1+r*snhi1
327  hri2=rri*snhi2+r*snhs2
328  hrs2=rrs*snhs2+r*snhi2
329  hri3=rri*snhi3+r*snhs3
330  hrs3=rrs*snhs3+r*snhi3
331 !
332  rr11=2*(rrri*(sni+nf1+nf1)+rrrs*(sns+nf4+nf4))
333  rr22=2*(rrri*(sni+nf2+nf2)+rrrs*(sns+nf5+nf5))
334  rr33=2*(rrri*(sni+nf3+nf3)+rrrs*(sns+nf6+nf6))
335  rr12= rrri*(sni+nf1+nf2)+rrrs*(sns+nf4+nf5)
336  rr13= rrri*(sni+nf1+nf3)+rrrs*(sns+nf4+nf6)
337  rr23= rrri*(sni+nf2+nf3)+rrrs*(sns+nf5+nf6)
338 !
339 ! EXTRA-DIAGONAL TERMS
340 !
341  xm(ielem, 1)=d*( hhi12 -d1*hri2-d2*hri1 +rr12)
342  xm(ielem, 2)=d*( hhi13 -d1*hri3-d3*hri1 +rr13)
343  xm(ielem, 3)=d*(-hh12-hh13+d1*(hri1-hrs1)-rr11)
344  xm(ielem, 4)=d*( hh12 +d1*hri2-d2*hrs1 -rr12)
345  xm(ielem, 5)=d*( hh13 +d1*hri3-d3*hrs1 -rr13)
346  xm(ielem, 6)=d*( hhi23 -d2*hri3-d3*hri2 +rr23)
347  xm(ielem, 7)=d*( hh12 +d2*hri1-d1*hrs2 -rr12)
348  xm(ielem, 8)=d*(-hh12-hh23+d2*(hri2-hrs2)-rr22)
349  xm(ielem, 9)=d*( hh23 +d2*hri3-d3*hrs2 -rr23)
350  xm(ielem,10)=d*( hh13 +d3*hri1-d1*hrs3 -rr13)
351  xm(ielem,11)=d*( hh23 +d3*hri2-d2*hrs3 -rr23)
352  xm(ielem,12)=d*(-hh13-hh23+d3*(hri3-hrs3)-rr33)
353  xm(ielem,13)=d*( hhs12 +d1*hrs2+d2*hrs1 +rr12)
354  xm(ielem,14)=d*( hhs13 +d1*hrs3+d3*hrs1 +rr13)
355  xm(ielem,15)=d*( hhs23 +d2*hrs3+d3*hrs2 +rr23)
356 !
357 !-----------------------------------------------------------------------
358 !
359 ! DIFFUSION ALONG Y
360 !
361 !-----------------------------------------------------------------------
362 !
363  d1=x(ielem,3)-x(ielem,2)
364 ! D2=X(IELEM,1)-X(IELEM,3)
365  d2= -x(ielem,3)
366 ! D3=X(IELEM,2)-X(IELEM,1)
367  d3=x(ielem,2)
368 !
369  ri=-(z(i1)*d1+z(i2)*d2+z(i3)*d3)
370  rs=-(z(i4)*d1+z(i5)*d2+z(i6)*d3)
371  r=ri+rs
372  rri=r+ri+ri
373  rrs=r+rs+rs
374  rrri=r*r+2*ri*ri
375  rrrs=r*r+2*rs*rs
376 !
377  d12=d1*d2
378  d13=d1*d3
379  d23=d2*d3
380 !
381  sni=ng1+ng2+ng3
382  sns=ng4+ng5+ng6
383  snhi=ng1*h1+ng2*h2+ng3*h3
384  snhs=ng4*h1+ng5*h2+ng6*h3
385  snhhi=(sni*sh1+snhi+snhi)*sh1+sni*shh
386  & +2*(ng1*h1*h1+ng2*h2*h2+ng3*h3*h3)
387  snhhs=(sns*sh1+snhs+snhs)*sh1+sns*shh
388  & +2*(ng4*h1*h1+ng5*h2*h2+ng6*h3*h3)
389  snhh=snhhi+snhhs
390  snhhi=snhh+snhhi+snhhi
391  snhhs=snhh+snhhs+snhhs
392 !
393  hhi12=snhhi*d12
394  hhi13=snhhi*d13
395  hhi23=snhhi*d23
396  hhs12=snhhs*d12
397  hhs13=snhhs*d13
398  hhs23=snhhs*d23
399  hh12=snhh*d12
400  hh13=snhh*d13
401  hh23=snhh*d23
402 !
403  snhi1=(sni+ng1)*(sh1+h1)+snhi+ng1*h1
404  snhs1=(sns+ng4)*(sh1+h1)+snhs+ng4*h1
405  snhi2=(sni+ng2)*(sh1+h2)+snhi+ng2*h2
406  snhs2=(sns+ng5)*(sh1+h2)+snhs+ng5*h2
407  snhi3=(sni+ng3)*(sh1+h3)+snhi+ng3*h3
408  snhs3=(sns+ng6)*(sh1+h3)+snhs+ng6*h3
409 !
410  hri1=rri*snhi1+r*snhs1
411  hrs1=rrs*snhs1+r*snhi1
412  hri2=rri*snhi2+r*snhs2
413  hrs2=rrs*snhs2+r*snhi2
414  hri3=rri*snhi3+r*snhs3
415  hrs3=rrs*snhs3+r*snhi3
416 !
417  rr11=2*(rrri*(sni+ng1+ng1)+rrrs*(sns+ng4+ng4))
418  rr22=2*(rrri*(sni+ng2+ng2)+rrrs*(sns+ng5+ng5))
419  rr33=2*(rrri*(sni+ng3+ng3)+rrrs*(sns+ng6+ng6))
420  rr12= rrri*(sni+ng1+ng2)+rrrs*(sns+ng4+ng5)
421  rr13= rrri*(sni+ng1+ng3)+rrrs*(sns+ng4+ng6)
422  rr23= rrri*(sni+ng2+ng3)+rrrs*(sns+ng5+ng6)
423 !
424 ! EXTRA-DIAGONAL TERMS
425 !
426  xm(ielem, 1)=xm(ielem, 1)+d*( hhi12 -d1*hri2-d2*hri1 +rr12)
427  xm(ielem, 2)=xm(ielem, 2)+d*( hhi13 -d1*hri3-d3*hri1 +rr13)
428  xm(ielem, 3)=xm(ielem, 3)+d*(-hh12-hh13+d1*(hri1-hrs1)-rr11)
429  xm(ielem, 4)=xm(ielem, 4)+d*( hh12 +d1*hri2-d2*hrs1 -rr12)
430  xm(ielem, 5)=xm(ielem, 5)+d*( hh13 +d1*hri3-d3*hrs1 -rr13)
431  xm(ielem, 6)=xm(ielem, 6)+d*( hhi23 -d2*hri3-d3*hri2 +rr23)
432  xm(ielem, 7)=xm(ielem, 7)+d*( hh12 +d2*hri1-d1*hrs2 -rr12)
433  xm(ielem, 8)=xm(ielem, 8)+d*(-hh12-hh23+d2*(hri2-hrs2)-rr22)
434  xm(ielem, 9)=xm(ielem, 9)+d*( hh23 +d2*hri3-d3*hrs2 -rr23)
435  xm(ielem,10)=xm(ielem,10)+d*( hh13 +d3*hri1-d1*hrs3 -rr13)
436  xm(ielem,11)=xm(ielem,11)+d*( hh23 +d3*hri2-d2*hrs3 -rr23)
437  xm(ielem,12)=xm(ielem,12)+d*(-hh13-hh23+d3*(hri3-hrs3)-rr33)
438  xm(ielem,13)=xm(ielem,13)+d*( hhs12 +d1*hrs2+d2*hrs1 +rr12)
439  xm(ielem,14)=xm(ielem,14)+d*( hhs13 +d1*hrs3+d3*hrs1 +rr13)
440  xm(ielem,15)=xm(ielem,15)+d*( hhs23 +d2*hrs3+d3*hrs2 +rr23)
441 !
442 !-----------------------------------------------------------------------
443 !
444 ! DIFFUSION ALONG Z
445 !
446 !-----------------------------------------------------------------------
447 !
448 ! VERSION WITH SIMPLIFICATIONS TO ACHIEVE MONOTONY OF THE MATRIX
449 !
450  d=surfac(ielem)*xs06
451 !
452 ! EXTRA-DIAGONAL TERMS
453 !
454  xm(ielem, 3)=xm(ielem, 3)-d*nh1
455  xm(ielem, 8)=xm(ielem, 8)-d*nh2
456  xm(ielem,12)=xm(ielem,12)-d*nh3
457 !
458 !-----------------------------------------------------------------------
459 !
460 ! OLD VERSION OF DIFFUSION ALONG Z
461 !
462 ! R=NH1+NH2+NH3+NH4+NH5+NH6
463 ! D=((X(I2)-X(I1))*(Y(I3)-Y(I1))-(X(I3)-X(I1))*(Y(I2)-Y(I1)))
464 ! * *XMUL/240.D0
465 !
466 ! RR11=(R+NH1+NH1+NH4+NH4)*(D+D)
467 ! RR22=(R+NH2+NH2+NH5+NH5)*(D+D)
468 ! RR33=(R+NH3+NH3+NH6+NH6)*(D+D)
469 ! RR12=(R+NH1+NH2+NH4+NH5)*D
470 ! RR13=(R+NH1+NH3+NH4+NH6)*D
471 ! RR23=(R+NH2+NH3+NH5+NH6)*D
472 !
473 ! EXTRA-DIAGONAL TERMS
474 !
475 ! XM(IELEM, 1)=XM(IELEM, 1)+RR12
476 ! XM(IELEM, 2)=XM(IELEM, 2)+RR13
477 ! XM(IELEM, 3)=XM(IELEM, 3)-RR11
478 ! XM(IELEM, 4)=XM(IELEM, 4)-RR12
479 ! XM(IELEM, 5)=XM(IELEM, 5)-RR13
480 ! XM(IELEM, 6)=XM(IELEM, 6)+RR23
481 ! XM(IELEM, 7)=XM(IELEM, 7)-RR12
482 ! XM(IELEM, 8)=XM(IELEM, 8)-RR22
483 ! XM(IELEM, 9)=XM(IELEM, 9)-RR23
484 ! XM(IELEM,10)=XM(IELEM,10)-RR13
485 ! XM(IELEM,11)=XM(IELEM,11)-RR23
486 ! XM(IELEM,12)=XM(IELEM,12)-RR33
487 ! XM(IELEM,13)=XM(IELEM,13)+RR12
488 ! XM(IELEM,14)=XM(IELEM,14)+RR13
489 ! XM(IELEM,15)=XM(IELEM,15)+RR23
490 !
491 !
492 !-----------------------------------------------------------------------
493 !
494  ENDDO
495 !
496 !-----------------------------------------------------------------------
497 !
498  ELSE
499 !
500 ! VERSION WITHOUT TREATMENT OF HYDROSTATIC INCONSISTENCIES
501 !
502 !-----------------------------------------------------------------------
503 !
504 ! DIFFUSION ALONG X
505 !
506 !-----------------------------------------------------------------------
507 !
508 ! LOOP ON THE ELEMENTS
509 !
510  DO ielem=1,nelem
511 !
512  i1=ikle(ielem,1)
513  i2=ikle(ielem,2)
514  i3=ikle(ielem,3)
515  i4=ikle(ielem,4)
516  i5=ikle(ielem,5)
517  i6=ikle(ielem,6)
518 ! DEPENDING ON NPOU0, II4 WILL BE I4 OR I1, ETC
519  ii4=i1+npou0
520  ii5=i2+npou0
521  ii6=i3+npou0
522 !
523  h1=z(i4)-z(i1)
524  h2=z(i5)-z(i2)
525  h3=z(i6)-z(i3)
526 !
527  sh1=h1+h2+h3
528  shh=h1*h1+h2*h2+h3*h3
529 !
530  d1=y(ielem,2)-y(ielem,3)
531 ! D2=Y(IELEM,3)-Y(IELEM,1)
532  d2=y(ielem,3)
533 ! D3=Y(IELEM,1)-Y(IELEM,2)
534  d3= -y(ielem,2)
535 !
536  d=xs2880/surfac(ielem)
537 !
538  ri=-(z(i1)*d1+z(i2)*d2+z(i3)*d3)
539  rs=-(z(i4)*d1+z(i5)*d2+z(i6)*d3)
540  r=ri+rs
541  rri=r+ri+ri
542  rrs=r+rs+rs
543  rrri=r*r+2*ri*ri
544  rrrs=r*r+2*rs*rs
545 !
546  d12=d1*d2
547  d13=d1*d3
548  d23=d2*d3
549 !
550  IF(h1.LT.chouia.OR.h2.LT.chouia.OR.h3.LT.chouia) THEN
551  nf1=0.d0
552  nf2=0.d0
553  nf3=0.d0
554  nf4=0.d0
555  nf5=0.d0
556  nf6=0.d0
557  ng1=0.d0
558  ng2=0.d0
559  ng3=0.d0
560  ng4=0.d0
561  ng5=0.d0
562  ng6=0.d0
563  nh1=0.d0
564  nh2=0.d0
565  nh3=0.d0
566  ELSE
567  nf1=f(i1)/h1
568  nf2=f(i2)/h2
569  nf3=f(i3)/h3
570  nf4=f(i4)/h1
571  nf5=f(i5)/h2
572  nf6=f(i6)/h3
573  ng1=g(i1)/h1
574  ng2=g(i2)/h2
575  ng3=g(i3)/h3
576  ng4=g(i4)/h1
577  ng5=g(i5)/h2
578  ng6=g(i6)/h3
579 ! DEPENDING ON THE CASE (II4=I1 OR I4, ETC.)
580 ! ALTERNATIVE WITH VERTICAL LINEAR VISCOSITY (II4=I4)
581 ! ALTERNATIVE WITH P0 VERTICAL VISCOSITY ON THE VERTICAL (II4=I1)
582  nh1=(h(i1)+h(ii4))/h1
583  nh2=(h(i2)+h(ii5))/h2
584  nh3=(h(i3)+h(ii6))/h3
585  ENDIF
586 !
587  sni=nf1+nf2+nf3
588  sns=nf4+nf5+nf6
589  snhi=nf1*h1+nf2*h2+nf3*h3
590  snhs=nf4*h1+nf5*h2+nf6*h3
591  snhhi=(sni*sh1+snhi+snhi)*sh1+sni*shh
592  & +2*(nf1*h1*h1+nf2*h2*h2+nf3*h3*h3)
593  snhhs=(sns*sh1+snhs+snhs)*sh1+sns*shh
594  & +2*(nf4*h1*h1+nf5*h2*h2+nf6*h3*h3)
595  snhh=snhhi+snhhs
596  snhhi=snhh+snhhi+snhhi
597  snhhs=snhh+snhhs+snhhs
598 !
599  hhi12=snhhi*d12
600  hhi13=snhhi*d13
601  hhi23=snhhi*d23
602  hhs12=snhhs*d12
603  hhs13=snhhs*d13
604  hhs23=snhhs*d23
605  hh12=snhh*d12
606  hh13=snhh*d13
607  hh23=snhh*d23
608 !
609  snhi1=(sni+nf1)*(sh1+h1)+snhi+nf1*h1
610  snhs1=(sns+nf4)*(sh1+h1)+snhs+nf4*h1
611  snhi2=(sni+nf2)*(sh1+h2)+snhi+nf2*h2
612  snhs2=(sns+nf5)*(sh1+h2)+snhs+nf5*h2
613  snhi3=(sni+nf3)*(sh1+h3)+snhi+nf3*h3
614  snhs3=(sns+nf6)*(sh1+h3)+snhs+nf6*h3
615 !
616  hri1=rri*snhi1+r*snhs1
617  hrs1=rrs*snhs1+r*snhi1
618  hri2=rri*snhi2+r*snhs2
619  hrs2=rrs*snhs2+r*snhi2
620  hri3=rri*snhi3+r*snhs3
621  hrs3=rrs*snhs3+r*snhi3
622 !
623  rr11=2*(rrri*(sni+nf1+nf1)+rrrs*(sns+nf4+nf4))
624  rr22=2*(rrri*(sni+nf2+nf2)+rrrs*(sns+nf5+nf5))
625  rr33=2*(rrri*(sni+nf3+nf3)+rrrs*(sns+nf6+nf6))
626  rr12= rrri*(sni+nf1+nf2)+rrrs*(sns+nf4+nf5)
627  rr13= rrri*(sni+nf1+nf3)+rrrs*(sns+nf4+nf6)
628  rr23= rrri*(sni+nf2+nf3)+rrrs*(sns+nf5+nf6)
629 !
630 ! EXTRA-DIAGONAL TERMS
631 !
632  xm(ielem, 1)=d*( hhi12 -d1*hri2-d2*hri1 +rr12)
633  xm(ielem, 2)=d*( hhi13 -d1*hri3-d3*hri1 +rr13)
634  xm(ielem, 3)=d*(-hh12-hh13+d1*(hri1-hrs1)-rr11)
635  xm(ielem, 4)=d*( hh12 +d1*hri2-d2*hrs1 -rr12)
636  xm(ielem, 5)=d*( hh13 +d1*hri3-d3*hrs1 -rr13)
637  xm(ielem, 6)=d*( hhi23 -d2*hri3-d3*hri2 +rr23)
638  xm(ielem, 7)=d*( hh12 +d2*hri1-d1*hrs2 -rr12)
639  xm(ielem, 8)=d*(-hh12-hh23+d2*(hri2-hrs2)-rr22)
640  xm(ielem, 9)=d*( hh23 +d2*hri3-d3*hrs2 -rr23)
641  xm(ielem,10)=d*( hh13 +d3*hri1-d1*hrs3 -rr13)
642  xm(ielem,11)=d*( hh23 +d3*hri2-d2*hrs3 -rr23)
643  xm(ielem,12)=d*(-hh13-hh23+d3*(hri3-hrs3)-rr33)
644  xm(ielem,13)=d*( hhs12 +d1*hrs2+d2*hrs1 +rr12)
645  xm(ielem,14)=d*( hhs13 +d1*hrs3+d3*hrs1 +rr13)
646  xm(ielem,15)=d*( hhs23 +d2*hrs3+d3*hrs2 +rr23)
647 !
648 !-----------------------------------------------------------------------
649 !
650 ! DIFFUSION ALONG Y
651 !
652 !-----------------------------------------------------------------------
653 !
654  d1=x(ielem,3)-x(ielem,2)
655 ! D2=X(IELEM,1)-X(IELEM,3)
656  d2= -x(ielem,3)
657 ! D3=X(IELEM,2)-X(IELEM,1)
658  d3=x(ielem,2)
659 !
660  ri=-(z(i1)*d1+z(i2)*d2+z(i3)*d3)
661  rs=-(z(i4)*d1+z(i5)*d2+z(i6)*d3)
662  r=ri+rs
663  rri=r+ri+ri
664  rrs=r+rs+rs
665  rrri=r*r+2*ri*ri
666  rrrs=r*r+2*rs*rs
667 !
668  d12=d1*d2
669  d13=d1*d3
670  d23=d2*d3
671 !
672  sni=ng1+ng2+ng3
673  sns=ng4+ng5+ng6
674  snhi=ng1*h1+ng2*h2+ng3*h3
675  snhs=ng4*h1+ng5*h2+ng6*h3
676  snhhi=(sni*sh1+snhi+snhi)*sh1+sni*shh
677  & +2*(ng1*h1*h1+ng2*h2*h2+ng3*h3*h3)
678  snhhs=(sns*sh1+snhs+snhs)*sh1+sns*shh
679  & +2*(ng4*h1*h1+ng5*h2*h2+ng6*h3*h3)
680  snhh=snhhi+snhhs
681  snhhi=snhh+snhhi+snhhi
682  snhhs=snhh+snhhs+snhhs
683 !
684  hhi12=snhhi*d12
685  hhi13=snhhi*d13
686  hhi23=snhhi*d23
687  hhs12=snhhs*d12
688  hhs13=snhhs*d13
689  hhs23=snhhs*d23
690  hh12=snhh*d12
691  hh13=snhh*d13
692  hh23=snhh*d23
693 !
694  snhi1=(sni+ng1)*(sh1+h1)+snhi+ng1*h1
695  snhs1=(sns+ng4)*(sh1+h1)+snhs+ng4*h1
696  snhi2=(sni+ng2)*(sh1+h2)+snhi+ng2*h2
697  snhs2=(sns+ng5)*(sh1+h2)+snhs+ng5*h2
698  snhi3=(sni+ng3)*(sh1+h3)+snhi+ng3*h3
699  snhs3=(sns+ng6)*(sh1+h3)+snhs+ng6*h3
700 !
701  hri1=rri*snhi1+r*snhs1
702  hrs1=rrs*snhs1+r*snhi1
703  hri2=rri*snhi2+r*snhs2
704  hrs2=rrs*snhs2+r*snhi2
705  hri3=rri*snhi3+r*snhs3
706  hrs3=rrs*snhs3+r*snhi3
707 !
708  rr11=2*(rrri*(sni+ng1+ng1)+rrrs*(sns+ng4+ng4))
709  rr22=2*(rrri*(sni+ng2+ng2)+rrrs*(sns+ng5+ng5))
710  rr33=2*(rrri*(sni+ng3+ng3)+rrrs*(sns+ng6+ng6))
711  rr12= rrri*(sni+ng1+ng2)+rrrs*(sns+ng4+ng5)
712  rr13= rrri*(sni+ng1+ng3)+rrrs*(sns+ng4+ng6)
713  rr23= rrri*(sni+ng2+ng3)+rrrs*(sns+ng5+ng6)
714 !
715 ! EXTRA-DIAGONAL TERMS
716 !
717  xm(ielem, 1)=xm(ielem, 1)+d*( hhi12 -d1*hri2-d2*hri1 +rr12)
718  xm(ielem, 2)=xm(ielem, 2)+d*( hhi13 -d1*hri3-d3*hri1 +rr13)
719  xm(ielem, 3)=xm(ielem, 3)+d*(-hh12-hh13+d1*(hri1-hrs1)-rr11)
720  xm(ielem, 4)=xm(ielem, 4)+d*( hh12 +d1*hri2-d2*hrs1 -rr12)
721  xm(ielem, 5)=xm(ielem, 5)+d*( hh13 +d1*hri3-d3*hrs1 -rr13)
722  xm(ielem, 6)=xm(ielem, 6)+d*( hhi23 -d2*hri3-d3*hri2 +rr23)
723  xm(ielem, 7)=xm(ielem, 7)+d*( hh12 +d2*hri1-d1*hrs2 -rr12)
724  xm(ielem, 8)=xm(ielem, 8)+d*(-hh12-hh23+d2*(hri2-hrs2)-rr22)
725  xm(ielem, 9)=xm(ielem, 9)+d*( hh23 +d2*hri3-d3*hrs2 -rr23)
726  xm(ielem,10)=xm(ielem,10)+d*( hh13 +d3*hri1-d1*hrs3 -rr13)
727  xm(ielem,11)=xm(ielem,11)+d*( hh23 +d3*hri2-d2*hrs3 -rr23)
728  xm(ielem,12)=xm(ielem,12)+d*(-hh13-hh23+d3*(hri3-hrs3)-rr33)
729  xm(ielem,13)=xm(ielem,13)+d*( hhs12 +d1*hrs2+d2*hrs1 +rr12)
730  xm(ielem,14)=xm(ielem,14)+d*( hhs13 +d1*hrs3+d3*hrs1 +rr13)
731  xm(ielem,15)=xm(ielem,15)+d*( hhs23 +d2*hrs3+d3*hrs2 +rr23)
732 !
733 !-----------------------------------------------------------------------
734 !
735 ! DIFFUSION ALONG Z
736 !
737 !-----------------------------------------------------------------------
738 !
739 ! VERSION WITH SIMPLIFICATIONS TO ACHIEVE MONOTONY OF THE MATRIX
740 !
741  d=surfac(ielem)*xs06
742 !
743 ! EXTRA-DIAGONAL TERMS
744 !
745  xm(ielem, 3)=xm(ielem, 3)-d*nh1
746  xm(ielem, 8)=xm(ielem, 8)-d*nh2
747  xm(ielem,12)=xm(ielem,12)-d*nh3
748 !
749 !-----------------------------------------------------------------------
750 !
751  ENDDO
752 !
753 ! IF(INCHYD) THEN
754  ENDIF
755 !
756 !-----------------------------------------------------------------------
757 !
758 ! POSITIVE EXTRA-DIAGONAL TERMS REMOVED
759 ! TO ENSURE MONOTONY
760 !
761  IF(formul(14:16).EQ.'MON') THEN
762 !
763  IF(xmul.GT.0.d0) THEN
764  DO ielem=1,nelem
765  xm(ielem, 1)=min(xm(ielem, 1),0.d0)
766  xm(ielem, 2)=min(xm(ielem, 2),0.d0)
767  xm(ielem, 3)=min(xm(ielem, 3),0.d0)
768  xm(ielem, 4)=min(xm(ielem, 4),0.d0)
769  xm(ielem, 5)=min(xm(ielem, 5),0.d0)
770  xm(ielem, 6)=min(xm(ielem, 6),0.d0)
771  xm(ielem, 7)=min(xm(ielem, 7),0.d0)
772  xm(ielem, 8)=min(xm(ielem, 8),0.d0)
773  xm(ielem, 9)=min(xm(ielem, 9),0.d0)
774  xm(ielem,10)=min(xm(ielem,10),0.d0)
775  xm(ielem,11)=min(xm(ielem,11),0.d0)
776  xm(ielem,12)=min(xm(ielem,12),0.d0)
777  xm(ielem,13)=min(xm(ielem,13),0.d0)
778  xm(ielem,14)=min(xm(ielem,14),0.d0)
779  xm(ielem,15)=min(xm(ielem,15),0.d0)
780  ENDDO
781  ELSE
782  DO ielem=1,nelem
783  xm(ielem, 1)=max(xm(ielem, 1),0.d0)
784  xm(ielem, 2)=max(xm(ielem, 2),0.d0)
785  xm(ielem, 3)=max(xm(ielem, 3),0.d0)
786  xm(ielem, 4)=max(xm(ielem, 4),0.d0)
787  xm(ielem, 5)=max(xm(ielem, 5),0.d0)
788  xm(ielem, 6)=max(xm(ielem, 6),0.d0)
789  xm(ielem, 7)=max(xm(ielem, 7),0.d0)
790  xm(ielem, 8)=max(xm(ielem, 8),0.d0)
791  xm(ielem, 9)=max(xm(ielem, 9),0.d0)
792  xm(ielem,10)=max(xm(ielem,10),0.d0)
793  xm(ielem,11)=max(xm(ielem,11),0.d0)
794  xm(ielem,12)=max(xm(ielem,12),0.d0)
795  xm(ielem,13)=max(xm(ielem,13),0.d0)
796  xm(ielem,14)=max(xm(ielem,14),0.d0)
797  xm(ielem,15)=max(xm(ielem,15),0.d0)
798  ENDDO
799  ENDIF
800  ENDIF
801 !
802 !-----------------------------------------------------------------------
803 !
804 ! DIAGONAL TERMS OBTAINED GIVEN THAT :
805 ! SUM OF THE TERMS IN A ROW =0
806 !
807 !-----------------------------------------------------------------------
808 !
809  DO ielem=1,nelem
810 !
811  t(ielem,1)= -xm(ielem,01)
812  & -xm(ielem,02)
813  & -xm(ielem,03)
814  & -xm(ielem,04)
815  & -xm(ielem,05)
816  t(ielem,2)= -xm(ielem,01)
817  & -xm(ielem,06)
818  & -xm(ielem,07)
819  & -xm(ielem,08)
820  & -xm(ielem,09)
821  t(ielem,3)= -xm(ielem,02)
822  & -xm(ielem,06)
823  & -xm(ielem,10)
824  & -xm(ielem,11)
825  & -xm(ielem,12)
826  t(ielem,4)= -xm(ielem,03)
827  & -xm(ielem,07)
828  & -xm(ielem,10)
829  & -xm(ielem,13)
830  & -xm(ielem,14)
831  t(ielem,5)= -xm(ielem,04)
832  & -xm(ielem,08)
833  & -xm(ielem,11)
834  & -xm(ielem,13)
835  & -xm(ielem,15)
836  t(ielem,6)= -xm(ielem,05)
837  & -xm(ielem,09)
838  & -xm(ielem,12)
839  & -xm(ielem,14)
840  & -xm(ielem,15)
841 !
842  ENDDO
843 !
844 !-----------------------------------------------------------------------
845 !
846  RETURN
847  END
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
mt02pp
subroutine mt02pp(T, XM, XMUL, SF, SG, SH, F, G, H, X, Y, Z, SURFAC, IKLE, NELEM, NELMAX, INCHYD, FORMUL, NPLAN)
Definition: mt02pp.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