1 : C ZCVI SOURCE CHAT 05/01/13 04:21:47 5004
2 : SUBROUTINE ZCVI
3 : C
4 : C
5 : C VERSION VECTORISEE
6 : C
7 : C Les {l{ments sont group{s en paquets de LRV {l{ments, LRV {tant
8 : C la longueur des registres vectoriels de la machine cible, i.e
9 : C 64 sur Cray, 128 ou 256 sur IBM 3090VF. On prom}ne une fenetre
10 : C de longueur LRV sur la boucle g{n{rale de longueur NEL.
11 : C
12 : C
13 : & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI,IKOMP,IKAS,
14 : & COEFF,IK1,RGE,IKG,NELG,TN,IKT,TREF,IKREF,IPADS,
15 : & UN,IPADU,NPTU,GN,F,IPADI,VF,IPADF,NPTF,
16 : & VOLU,COTE,NELZ,IDCEN,IPG,
17 : & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,FN)
18 :
19 : IMPLICIT INTEGER(I-N)
20 : IMPLICIT REAL*8 (A-H,O-Z)
21 :
22 : C***********************************************************************
23 : C
24 : C CE SP DISCRETISE LES EQUATIONS DE NAVIER STOKES
25 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI
26 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6
27 : C LES OPERATEURS SONT "SOUS-INTEGRES"
28 : C
29 : C SYNTAXE :
30 : C
31 : C NS(NU,UN,RGE,DE) INCO GN :
32 : C
33 : C COEFFICIENTS :
34 : C --------------
35 : C
36 : C UN(NPTD,IES) CHAMPS DE VITESSE TRANSPORTANT
37 : C COEFF(SCAL DOMA) VISCOSITE CINEMATIQE MOLECULAIRE( NU )
38 : C (SCAL ELEM)
39 : C RGE(NELG,IES) TERME SOURCE
40 : C
41 : C
42 : C INCONNUES :
43 : C -----------
44 : C
45 : C UN(NPTD,IES) CHAMPS DE VITESSE TRANSPORTANT
46 : C GN(NPTD,IES) CHAMPS DE VITESSE TRANSPORTE
47 : C VN(NPTD,IES) CHAMPS DE VITESSE DU FLUIDE
48 : C
49 : C
50 : C
51 : C
52 : C TABLEAUX DE TRAVAIL :
53 : C --------------------- N N-1
54 : C (D + D1) U - D U N-1 T N
55 : C ------------------- = F - A U - C P
56 : C DT
57 : C
58 : C N-1
59 : C F(NPTD,IES) CONTIENT A U - F (VITESSE)
60 : C
61 : C
62 : C
63 : C***********************************************************************
64 :
65 : -INC CCVQUA4
66 : -INC CCREEL
67 : *-
68 : -INC CCOPTIO
69 : -INC SMCOORD
70 : C
71 : C Longueur des registres vectoriels de la machine cible
72 : C On prend 64 pour ne pas augmenter la taille des tableaux
73 : C n{cessaires @ la vectorisation.
74 : C
75 : PARAMETER(LRV=64)
76 :
77 : DIMENSION UN(NPTU,IES),GN(NPTD,IES),VF(NPTF,IES)
78 : DIMENSION TN(*),TREF(*)
79 : DIMENSION COEFF(*),RGE(NELG,IES)
80 : DIMENSION COTE(NELZ,IES),VOLU(NELZ),KLIP(100)
81 :
82 : DIMENSION IPADI(*),LE(NP,1),IPADU(*),IPADF(*),IPADS(*)
83 : DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL)
84 :
85 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
86 :
87 : DIMENSION COEF(LRV),AIRE(LRV)
88 : DIMENSION WX(LRV,9),WY(LRV,9),WZ(LRV,9)
89 : DIMENSION AL(LRV),AH(LRV),AP(LRV)
90 : C UIX,... vitesse transportante
91 : DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9)
92 : C GIX,... vitesse massique ou transportée ou inconnue du/dt
93 : DIMENSION GIX(LRV,9),GIY(LRV,9),GIZ(LRV,9)
94 : C VIX,... vitesse du fluide
95 : DIMENSION VIX(LRV,9),VIY(LRV,9),VIZ(LRV,9)
96 :
97 : DIMENSION UMI(LRV,3),VMI(LRV,3)
98 : DIMENSION COEFT(LRV),RGX(LRV),RGY(LRV),RGZ(LRV)
99 : DIMENSION TO1(LRV),TO2X(LRV),TO2Y(LRV)
100 : DIMENSION SAF1(LRV,9),SAF2(LRV,9),SAF3(LRV,9)
101 : DIMENSION CHGLD(LRV),CHGLPX(LRV),CHGLPY(LRV),CHGLPZ(LRV)
102 :
103 : DIMENSION F(NPTD,*),FN(NP,*)
104 :
105 : SAVE IPAS,QGGT,Q1,Q2,Q3
106 : DATA IPAS/0/
107 :
108 : C
109 : C INITIALISATIONS DIVERSES
110 : C
111 : C write(6,*)' DEBUT YCVI ',' IKAS=',ikas,' IKOMP=',ikomp,idcen
112 : C write(6,*)' NPTD=',nptd,IPAS
112 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
113 :
114 : C ********
115 : C * 2D *
116 : C ********
117 :
118 : IF(IES.EQ.3)GO TO 10
119 :
120 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
121 : C DIFFERENCES TRIANGLE / QUADRANGLE
122 : IF(NP.EQ.4)THEN
123 : QUA4=1.D0
124 : ELSE
125 : QUA4=0.D0
126 : ENDIF
127 : C
128 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
129 :
130 : C
131 : C Calcul du nombre de paquets de LRV {l{ments
132 : C
133 : NNN=MOD(NEL,LRV)
134 : IF(NNN.EQ.0) NPACK=NEL/LRV
135 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
136 : KPACKD=1
137 : KPACKF=NPACK
138 : C
139 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
140 : C
141 : C write(6,*)' DEBUT YCVI 7001'
142 : DO 7001 KPACK=KPACKD,KPACKF
143 : C
144 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
145 : C
146 : C 1. Calcul des limites du paquet courant.
147 : KDEB=1+(KPACK-1)*LRV
148 : KFIN=MIN(NEL,KDEB+LRV-1)
149 : C
150 : C Ben voil@, on peut y aller ... i.e. traiter le paquet courant.
151 : C
152 : DO 7002 K=KDEB,KFIN
153 : KP=K-KDEB+1
154 : NK=K+K0
155 : K1=1+(1-IK1)*(NK-1)
156 : COEF(KP)=COEFF(K1)
157 : AIRE(KP)=VOLU(NK)
158 : AL(KP)=COTE(NK,1)+XPETIT
159 : AH(KP)=COTE(NK,2)+XPETIT
160 : 7002 CONTINUE
161 :
162 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
163 : DO 7006 I=1,NP
164 : DO 7006 K=KDEB,KFIN
165 : KP=K-KDEB+1
166 : NU=IPADU(LE(I,K))
167 : NG=IPADI(LE(I,K))
168 : NF=IPADF(LE(I,K))
169 : UIX(KP,I)=UN(NU,1)
170 : UIY(KP,I)=UN(NU,2)
171 : GIX(KP,I)=GN(NG,1)
172 : GIY(KP,I)=GN(NG,2)
173 : VIX(KP,I)=VF(NF,1)
174 : VIY(KP,I)=VF(NF,2)
175 : 7006 CONTINUE
176 :
176 : CALL KSUPG1(WX,UMI,CHGLD,CHGLPX,KDEB,KFIN,LRV,
177 : &GN(1,1),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
178 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,QUA4,IKOMP,
179 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
180 :
180 : CALL KSUPG1(WY,UMI,CHGLD,CHGLPY,KDEB,KFIN,LRV,
181 : &GN(1,2),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
182 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,QUA4,IKOMP,
183 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
184 :
185 : C
186 : C Initialisation des variables d'accumulation SAF1,SAF2,SBF
187 : C
188 :
189 : IF(IAXI.NE.0)THEN
190 : DO 70050 N=1,IDIM
191 : DO 70050 K=KDEB,KFIN
192 : KP=K-KDEB+1
193 : VMI(KP,N)=XPETIT
194 : 70050 CONTINUE
195 : DO 70051 N=1,IDIM
196 : DO 70051 I=1,NP
197 : DO 70051 K=KDEB,KFIN
198 : KP=K-KDEB+1
199 : NF=IPADF(LE(I,K))
200 : VMI(KP,N)=VMI(KP,N)+VF(NF,N)*DRR(I,K)
201 : 70051 CONTINUE
202 : DO 70052 K=KDEB,KFIN
203 : KP=K-KDEB+1
204 : VMI(KP,1)=VMI(KP,1)/AIRE(KP)
205 : VMI(KP,2)=VMI(KP,2)/AIRE(KP)
206 : 70052 CONTINUE
207 : ENDIF
208 :
209 : IF(IKOMP.EQ.0)THEN
210 :
211 : IF(IKAS.EQ.1)THEN
212 : DO 70061 I=1,NP
213 : DO 70061 K=KDEB,KFIN
214 : KP=K-KDEB+1
215 : SAF1(KP,I)=0.D0
216 : SAF2(KP,I)=0.D0
217 : 70061 CONTINUE
218 : ELSEIF(IKAS.EQ.2)THEN
219 : DO 70021 K=KDEB,KFIN
220 : KP=K-KDEB+1
221 : NK=K+K0
222 : NKG=1+(1-IKG)*(NK-1)
223 : RGX(KP)=RGE(NKG,1)
224 : RGY(KP)=RGE(NKG,2)
225 : 70021 CONTINUE
226 :
227 : IF(IPG.EQ.0)THEN
228 : DO 70062 I=1,NP
229 : DO 70062 K=KDEB,KFIN
230 : KP=K-KDEB+1
231 : SAF1(KP,I)=(-RGX(KP))*DRR(I,K)
232 : SAF2(KP,I)=(-RGY(KP))*DRR(I,K)
233 : 70062 CONTINUE
234 : ELSE
235 : DO 71062 I=1,NP
236 : DO 71062 K=KDEB,KFIN
237 : KP=K-KDEB+1
238 : SAF1(KP,I)=(-RGX(KP))*WX(KP,I)
239 : SAF2(KP,I)=(-RGY(KP))*WY(KP,I)
240 : 71062 CONTINUE
241 : ENDIF
242 :
243 : ELSEIF(IKAS.EQ.4)THEN
244 : DO 70022 K=KDEB,KFIN
245 : KP=K-KDEB+1
246 : NK=K+K0
247 : NKG=1+(1-IKG)*(NK-1)
248 : RGX(KP)=RGE(NKG,1)
249 : RGY(KP)=RGE(NKG,2)
250 : 70022 CONTINUE
251 :
252 : IF(IPG.EQ.0)THEN
253 : DO 70063 I=1,NP
254 : DO 70063 K=KDEB,KFIN
255 : KP=K-KDEB+1
256 : NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
257 : NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
258 : C? write(6,*)' NF=',NF,' NFR=',nfr
259 : SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
260 : SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
261 : 70063 CONTINUE
262 : ELSE
263 : DO 71063 I=1,NP
264 : DO 71063 K=KDEB,KFIN
265 : KP=K-KDEB+1
266 : NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
267 : NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
268 : SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*WX(KP,I)
269 : SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*WY(KP,I)
270 : 71063 CONTINUE
271 : ENDIF
272 :
273 : ENDIF
274 :
275 : ELSEIF(IKOMP.EQ.1)THEN
276 :
277 : IF(IKAS.EQ.2)THEN
278 : DO 70064 I=1,NP
279 : DO 70064 K=KDEB,KFIN
280 : KP=K-KDEB+1
281 : SAF1(KP,I)=0.D0
282 : SAF2(KP,I)=0.D0
283 : 70064 CONTINUE
284 :
285 : ELSEIF(IKAS.EQ.3)THEN
286 : DO 70024 K=KDEB,KFIN
287 : KP=K-KDEB+1
288 : NK=K+K0
289 : NKG=1+(1-IKG)*(NK-1)
290 : RGX(KP)=RGE(NKG,1)
291 : RGY(KP)=RGE(NKG,2)
292 : 70024 CONTINUE
293 :
294 : IF(IPG.EQ.0)THEN
295 : DO 70065 I=1,NP
296 : DO 70065 K=KDEB,KFIN
297 : KP=K-KDEB+1
298 : SAF1(KP,I)=-RGX(KP)*DRR(I,K)
299 : SAF2(KP,I)=-RGY(KP)*DRR(I,K)
300 : 70065 CONTINUE
301 : ELSE
302 : DO 71065 I=1,NP
303 : DO 71065 K=KDEB,KFIN
304 : KP=K-KDEB+1
305 : SAF1(KP,I)=-RGX(KP)*WX(KP,I)
306 : SAF2(KP,I)=-RGY(KP)*WY(KP,I)
307 : 71065 CONTINUE
308 : ENDIF
309 :
310 : ENDIF
311 :
312 : ENDIF
313 :
314 : C Le coeur du calcul ...
315 :
316 : IF(IKOMP.EQ.0)THEN
317 :
318 : DO 7014 I=1,NP
319 : DO 7014 J= 1,NP
320 : DO 7014 K=KDEB,KFIN
321 : KP=K-KDEB+1
322 :
323 : ZVGGX=AIRE(KP)*CHGLPX(KP)*VGGT(J,I)
324 :
325 : ZVGGY=AIRE(KP)*CHGLPY(KP)*VGGT(J,I)
326 :
327 : ZVGT=AIRE(KP)*(
328 : & HR(K,I,1)*HR(K,J,1)*COEF(KP)
329 : &+ HR(K,I,2)*HR(K,J,2)*COEF(KP)
330 : &+ CHGLD(KP)*VGGT(J,I) )
331 :
332 : V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)
333 :
334 : SAF1(KP,I)=SAF1(KP,I)+(V2*WX(KP,I)+ZVGGX)*GIX(KP,J)+ZVGT*VIX(KP,J)
335 : SAF2(KP,I)=SAF2(KP,I)+(V2*WY(KP,I)+ZVGGY)*GIY(KP,J)+ZVGT*VIY(KP,J)
336 :
337 : 7014 CONTINUE
338 :
339 : ELSEIF(IKOMP.EQ.1)THEN
340 :
341 : DO 7015 I=1,NP
342 : DO 7015 J= 1,NP
343 : DO 7015 K=KDEB,KFIN
344 : KP=K-KDEB+1
345 : ZVGGX=AIRE(KP)*CHGLPX(KP)*VGGT(J,I)
346 :
347 : ZVGGY=AIRE(KP)*CHGLPY(KP)*VGGT(J,I)
348 :
349 : ZVGT=AIRE(KP)*(
350 : & HR(K,I,1)*HR(K,J,1)*COEF(KP)
351 : &+ HR(K,I,2)*HR(K,J,2)*COEF(KP)
352 : &+ CHGLD(KP)*VGGT(J,I) )
353 :
354 : COEFT3=COEF(KP)/3.D0
355 : ZVGUU=AIRE(KP)*(1.D0/AL(KP)/AL(KP)/12.D0
356 : & *VGGT(J,I)*QUA4 )*COEFT3*UIX(KP,J)
357 : ZVGUV=AIRE(KP)*( HR(K,I,1)*HR(K,J,2))*COEFT3*UIY(KP,J)
358 : ZVGVU=AIRE(KP)*( HR(K,I,2)*HR(K,J,1))*COEFT3*UIX(KP,J)
359 : ZVGVV=AIRE(KP)*(1.D0/AH(KP)/AH(KP)/12.D0
360 : & *VGGT(J,I)*QUA4 )*COEFT3*UIY(KP,J)
361 :
362 : V2=UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)
363 :
364 : SAF1(KP,I)=SAF1(KP,I)+(V2*WX(KP,I)+ZVGGX)*GIX(KP,J)+ZVGT*UIX(KP,J)
365 : &+ ZVGUU + ZVGUV
366 : SAF2(KP,I)=SAF2(KP,I)+(V2*WY(KP,I)+ZVGGY)*GIY(KP,J)+ZVGT*UIY(KP,J)
367 : &+ ZVGVU + ZVGVV
368 :
369 : 7015 CONTINUE
370 :
371 : ENDIF
372 :
373 : IF(IAXI.NE.0) THEN
374 :
375 : DO 7016 I=1,NP
376 : DO 7016 K=KDEB,KFIN
377 : KP=K-KDEB+1
378 :
379 : R2=1.D0/RPG(K)/RPG(K)*WX(KP,I)
380 : SAF1(KP,I)=SAF1(KP,I)+R2*COEF(KP)*VMI(KP,1)
381 :
382 : 7016 CONTINUE
383 :
384 : IF(IKOMP.EQ.1)THEN
385 : DO 7118 I=1,NP
386 : DO 7118 K=KDEB,KFIN
387 : KP=K-KDEB+1
388 : R1=1.D0/RPG(K)*WY(KP,I)
389 : SAF1(KP,I)=SAF1(KP,I)+R1*UMI(KP,1)*GIX(KP,I)
390 : SAF2(KP,I)=SAF2(KP,I)+R1*UMI(KP,1)*GIY(KP,I)
391 : 7118 CONTINUE
392 : ENDIF
393 :
394 : ENDIF
395 : C
396 : C Fin de l'accumulation dans SAF1,SAF2.
397 : C On ajoute ces incr{ments @ F.
398 : C
399 : DO 7017 I=1,NP
400 : DO 7017 K=KDEB,KFIN
401 : KP=K-KDEB+1
402 : NF=IPADI(LE(I,K))
403 : F(NF,1)=F(NF,1)+SAF1(KP,I)
404 : F(NF,2)=F(NF,2)+SAF2(KP,I)
405 : 7017 CONTINUE
406 :
407 : 1960 FORMAT(/,' ***** SUB XCVTIT : IPAT=',I5,' K=',I5,' *****')
408 : 1961 FORMAT(2X,I5,' * ',4(1X,I5))
409 : 1962 FORMAT(2X,8(1X,1PE11.4))
410 : 1964 FORMAT(4(1X,1PE11.4))
411 : 7001 CONTINUE
412 :
413 : C WRITE(6,*)' ********** FIN YCVI 2D *****************'
413 : C CALL ARRET(0)
414 : IPAS=1
415 : RETURN
416 :
417 : C ********
418 : C * 3D *
419 : C ********
420 :
421 : 10 CONTINUE
422 :
423 : C::::::BENET:::SUPPRESION CORRECTION HOURGLASS POUR LES PRISME::29:01:91
424 :
425 : CUB8=0.D0
426 : IF(NP.EQ.8)CUB8=1.D0
427 :
428 : C
429 : C Calcul du nombre de paquets de LRV {l{ments
430 : C
431 : NNN=MOD(NEL,LRV)
432 : IF(NNN.EQ.0) NPACK=NEL/LRV
433 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
434 : KPACKD=1
435 : KPACKF=NPACK
436 : C
437 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
438 : C
439 : DO 8001 KPACK=KPACKD,KPACKF
440 : C
441 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
442 : C
443 : C 1. Calcul des limites du paquet courant.
444 : KDEB=1+(KPACK-1)*LRV
445 : KFIN=MIN(NEL,KDEB+LRV-1)
446 : C
447 : C Ben voil@, on peut y aller ... i.e. traiter le paquet courant.
448 : C
449 : DO 8002 K=KDEB,KFIN
450 : KP=K-KDEB+1
451 : NK=K+K0
452 : K1=1+(1-IK1)*(NK-1)
453 : COEF(KP)=COEFF(K1)
454 : AIRE(KP)=VOLU(NK)
455 : AL(KP)=COTE(NK,1)+XPETIT
456 : AH(KP)=COTE(NK,2)+XPETIT
457 : AP(KP)=COTE(NK,3)+XPETIT
458 : 8002 CONTINUE
459 :
460 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
461 : C Initialisation des UMI avant accumulation
462 :
463 : DO 8006 I=1,NP
464 : DO 8006 K=KDEB,KFIN
465 : KP=K-KDEB+1
466 : NU=IPADU(LE(I,K))
467 : NF=IPADF(LE(I,K))
468 : NG=IPADI(LE(I,K))
469 : UIX(KP,I)=UN(NU,1)
470 : UIY(KP,I)=UN(NU,2)
471 : UIZ(KP,I)=UN(NU,3)
472 : VIX(KP,I)=VF(NF,1)
473 : VIY(KP,I)=VF(NF,2)
474 : VIZ(KP,I)=VF(NF,3)
475 : GIX(KP,I)=GN(NG,1)
476 : GIY(KP,I)=GN(NG,2)
477 : GIZ(KP,I)=GN(NG,3)
478 : 8006 CONTINUE
479 :
480 : C write(6,*)'****************************'
481 : C write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
481 : CALL KSUPG1(WX,UMI,CHGLD,CHGLPX,KDEB,KFIN,LRV,
482 : &GN(1,1),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
483 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP,
484 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
485 : C write(6,1002)umi
486 :
487 : C write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
487 : CALL KSUPG1(WY,UMI,CHGLD,CHGLPY,KDEB,KFIN,LRV,
488 : &GN(1,2),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
489 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP,
490 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
491 :
492 : C write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
492 : CALL KSUPG1(WZ,UMI,CHGLD,CHGLPZ,KDEB,KFIN,LRV,
493 : &GN(1,3),IPADI,UN,COEF,NPTD,NEL,NP,DRR,HR,FN,
494 : &AIRE,AL,AH,AP,IDCEN,IPADU,LE,CUB8,IKOMP,
495 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
496 : C write(6,*)' DT,DTT1,DTT2=',DT,DTT1,DTT2
497 :
498 : C
499 : C Initialisation des variables d'accumulation SAF1,SAF2,SBF
500 : C
501 : IF(IKOMP.EQ.0)THEN
502 :
503 : IF(IKAS.EQ.1)THEN
504 : DO 80061 I=1,NP
505 : DO 80061 K=KDEB,KFIN
506 : KP=K-KDEB+1
507 : SAF1(KP,I)=0.D0
508 : SAF2(KP,I)=0.D0
509 : SAF3(KP,I)=0.D0
510 : 80061 CONTINUE
511 : ELSEIF(IKAS.EQ.2)THEN
512 : DO 80021 K=KDEB,KFIN
513 : KP=K-KDEB+1
514 : NK=K+K0
515 : NKG=1+(1-IKG)*(NK-1)
516 : RGX(KP)=RGE(NKG,1)
517 : RGY(KP)=RGE(NKG,2)
518 : RGZ(KP)=RGE(NKG,3)
519 : 80021 CONTINUE
520 :
521 : IF(IPG.EQ.0)THEN
522 : DO 80062 I=1,NP
523 : DO 80062 K=KDEB,KFIN
524 : KP=K-KDEB+1
525 : SAF1(KP,I)=(-RGX(KP))*DRR(I,K)
526 : SAF2(KP,I)=(-RGY(KP))*DRR(I,K)
527 : SAF3(KP,I)=(-RGZ(KP))*DRR(I,K)
528 : 80062 CONTINUE
529 : ELSE
530 : DO 81062 I=1,NP
531 : DO 81062 K=KDEB,KFIN
532 : KP=K-KDEB+1
533 : SAF1(KP,I)=(-RGX(KP))*WX(KP,I)
534 : SAF2(KP,I)=(-RGY(KP))*WY(KP,I)
535 : SAF3(KP,I)=(-RGZ(KP))*WZ(KP,I)
536 : 81062 CONTINUE
537 : ENDIF
538 :
539 : ELSEIF(IKAS.EQ.4)THEN
540 : DO 80022 K=KDEB,KFIN
541 : KP=K-KDEB+1
542 : NK=K+K0
543 : NKG=1+(1-IKG)*(NK-1)
544 : RGX(KP)=RGE(NKG,1)
545 : RGY(KP)=RGE(NKG,2)
546 : RGZ(KP)=RGE(NKG,3)
547 : 80022 CONTINUE
548 :
549 : IF(IPG.EQ.0)THEN
550 : DO 80063 I=1,NP
551 : DO 80063 K=KDEB,KFIN
552 : KP=K-KDEB+1
553 : NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
554 : NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
555 : SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
556 : SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
557 : SAF3(KP,I)=(-RGZ(KP)*(TREF(NFR)-TN(NF)))*DRR(I,K)
558 : 80063 CONTINUE
559 : ELSE
560 : DO 81063 I=1,NP
561 : DO 81063 K=KDEB,KFIN
562 : KP=K-KDEB+1
563 : NF=1+(1-IKT)*(IPADS(LE(I,K))-1)
564 : NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1)
565 : SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*WX(KP,I)
566 : SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*WY(KP,I)
567 : SAF3(KP,I)=(-RGZ(KP)*(TREF(NFR)-TN(NF)))*WZ(KP,I)
568 : 81063 CONTINUE
569 : ENDIF
570 :
571 : ENDIF
572 :
573 : ELSEIF(IKOMP.EQ.1)THEN
574 :
575 : IF(IKAS.EQ.2)THEN
576 : DO 80064 I=1,NP
577 : DO 80064 K=KDEB,KFIN
578 : KP=K-KDEB+1
579 : SAF1(KP,I)=0.D0
580 : SAF2(KP,I)=0.D0
581 : SAF3(KP,I)=0.D0
582 : 80064 CONTINUE
583 :
584 : ELSEIF(IKAS.EQ.3)THEN
585 : DO 80024 K=KDEB,KFIN
586 : KP=K-KDEB+1
587 : NK=K+K0
588 : NKG=1+(1-IKG)*(NK-1)
589 : RGX(KP)=RGE(NKG,1)
590 : RGY(KP)=RGE(NKG,2)
591 : RGZ(KP)=RGE(NKG,3)
592 : 80024 CONTINUE
593 :
594 : IF(IPG.EQ.0)THEN
595 : DO 80065 I=1,NP
596 : DO 80065 K=KDEB,KFIN
597 : KP=K-KDEB+1
598 : SAF1(KP,I)=-RGX(KP)*DRR(I,K)
599 : SAF2(KP,I)=-RGY(KP)*DRR(I,K)
600 : SAF3(KP,I)=-RGZ(KP)*DRR(I,K)
601 : 80065 CONTINUE
602 : ELSE
603 : DO 81065 I=1,NP
604 : DO 81065 K=KDEB,KFIN
605 : KP=K-KDEB+1
606 : SAF1(KP,I)=-RGX(KP)*WX(KP,I)
607 : SAF2(KP,I)=-RGY(KP)*WY(KP,I)
608 : SAF3(KP,I)=-RGZ(KP)*WZ(KP,I)
609 : 81065 CONTINUE
610 : ENDIF
611 :
612 : ENDIF
613 :
614 : ENDIF
615 :
616 : C Le coeur du calcul ...
617 :
618 : IF(IKOMP.EQ.0)THEN
619 :
620 : DO 8014 I=1,NP
621 : DO 8014 J= 1,NP
622 : DO 8014 K=KDEB,KFIN
623 : KP=K-KDEB+1
624 :
625 : ZVGGX=AIRE(KP)*CHGLPX(KP)*QGGT(J,I)
626 :
627 : ZVGGY=AIRE(KP)*CHGLPY(KP)*QGGT(J,I)
628 :
629 : ZVGGZ=AIRE(KP)*CHGLPZ(KP)*QGGT(J,I)
630 :
631 : ZVGT=AIRE(KP)*(
632 : & HR(K,I,1)*HR(K,J,1)*COEF(KP)
633 : &+ HR(K,I,2)*HR(K,J,2)*COEF(KP)
634 : &+ HR(K,I,3)*HR(K,J,3)*COEF(KP) )
635 : &+ CHGLD(KP)*QGGT(J,I)
636 :
637 : V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)
638 : & +UMI(KP,3)*HR(K,J,3)
639 :
640 : SAF1(KP,I)=SAF1(KP,I)+(V2*WX(KP,I)+ZVGGX)*GIX(KP,J)+ZVGT*VIX(KP,J)
641 : SAF2(KP,I)=SAF2(KP,I)+(V2*WY(KP,I)+ZVGGY)*GIY(KP,J)+ZVGT*VIY(KP,J)
642 : SAF3(KP,I)=SAF3(KP,I)+(V2*WZ(KP,I)+ZVGGZ)*GIZ(KP,J)+ZVGT*VIZ(KP,J)
643 :
644 : 8014 CONTINUE
645 :
646 : ELSEIF(IKOMP.EQ.1)THEN
647 :
648 : DO 8015 I=1,NP
649 : DO 8015 J= 1,NP
650 : DO 8015 K=KDEB,KFIN
651 : KP=K-KDEB+1
652 :
653 : ZVGGX=AIRE(KP)*CHGLPX(KP)*QGGT(J,I)
654 :
655 : ZVGGY=AIRE(KP)*CHGLPY(KP)*QGGT(J,I)
656 :
657 : ZVGGZ=AIRE(KP)*CHGLPZ(KP)*QGGT(J,I)
658 :
659 : ZVGT=AIRE(KP)*(
660 : & HR(K,I,1)*HR(K,J,1)*COEF(KP)
661 : &+ HR(K,I,2)*HR(K,J,2)*COEF(KP)
662 : &+ HR(K,I,3)*HR(K,J,3)*COEF(KP) )
663 : &+ CHGLD(KP)*QGGT(J,I)
664 :
665 : COEFT3=COEF(KP)/3.D0
666 : GEO1=0.D0
667 : ZVGUU=(AIRE(KP)*(HR(K,I,1)*HR(K,J,1))+GEO1)*COEFT3*UIX(KP,J)
668 : ZVGUV=AIRE(KP)*( HR(K,I,1)*HR(K,J,2))*COEFT3*UIY(KP,J)
669 : ZVGUW=AIRE(KP)*( HR(K,I,1)*HR(K,J,3))*COEFT3*UIZ(KP,J)
670 : ZVGVU=AIRE(KP)*( HR(K,I,2)*HR(K,J,1))*COEFT3*UIX(KP,J)
671 : ZVGVV=(AIRE(KP)*(HR(K,I,2)*HR(K,J,2))+GEO1)*COEFT3*UIY(KP,J)
672 : ZVGVW=AIRE(KP)*( HR(K,I,2)*HR(K,J,3))*COEFT3*UIZ(KP,J)
673 : ZVGWU=AIRE(KP)*( HR(K,I,3)*HR(K,J,1))*COEFT3*UIX(KP,J)
674 : ZVGWV=AIRE(KP)*( HR(K,I,3)*HR(K,J,2))*COEFT3*UIY(KP,J)
675 : ZVGWW=(AIRE(KP)*(HR(K,I,3)*HR(K,J,3))+GEO1)*COEFT3*UIZ(KP,J)
676 :
677 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+UIZ(KP,J)*HR(K,J,3))
678 : & *WZ(KP,I)
679 :
680 : SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGT*UIX(KP,J)
681 : &+ ZVGUU + ZVGUV + ZVGUW
682 : SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGT*UIY(KP,J)
683 : &+ ZVGVU + ZVGUV + ZVGVW
684 : SAF3(KP,I)=SAF3(KP,I)+(V2+ZVGGZ)*GIZ(KP,J)+ZVGT*UIZ(KP,J)
685 : &+ ZVGWU + ZVGWV + ZVGWW
686 :
687 : 8015 CONTINUE
688 :
689 : ENDIF
690 :
691 : C
692 : C Fin de l'accumulation dans SAF1,SAF2.
693 : C On ajoute ces incr{ments @ F.
694 : C
695 : DO 8017 I=1,NP
696 : DO 8017 K=KDEB,KFIN
697 : KP=K-KDEB+1
698 : NF=IPADI(LE(I,K))
699 : F(NF,1)=F(NF,1)+SAF1(KP,I)
700 : F(NF,2)=F(NF,2)+SAF2(KP,I)
701 : F(NF,3)=F(NF,3)+SAF3(KP,I)
702 : 8017 CONTINUE
703 :
704 : 8001 CONTINUE
705 :
706 : C WRITE(6,*)' ********** FIN YCVI 3D *****************'
707 :
708 : IPAS=1
709 : RETURN
710 : 1002 FORMAT(10(1X,1PE11.4))
711 : 1001 FORMAT(20(1X,I5))
712 : END
713 :
714 :
715 :
716 :
717 :
718 :
719 :
© Cast3M 2003 - All rights reserved.
Disclaimer