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