1 : C ZJOHNS SOURCE CHAT 05/01/13 04:23:03 5004 2 : SUBROUTINE ZJOHNS 3 : C 4 : C VERSION VECTORISEE 5 : C 6 : C Les éléments sont groupés en paquets de LRV éléments, LRV étant 7 : C la longueur des registres vectoriels de la machine cible, i.e 8 : C 64 sur Cray, 128 ou 256 sur IBM 3090VF. On promène une fenêtre 9 : C de longueur LRV sur la boucle générale de longueur NEL. 10 : C 11 : 12 : C KDESIGN n'est pas défini !!!!!! 13 : C 14 : C 15 : & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI, 16 : & IPADL,IKOMP,IKAS, 17 : & ALFE,IND1,UN,INDU,NPTU,IPADU, 18 : & TN,QE,IKS, 19 : & HRN,G,ALT,SGT, 20 : & VOLU,COTE,NELZ, 21 : & DTM1,DT,DTT1,DTT2,NUEL,DIAEL) 22 : 23 : IMPLICIT INTEGER(I-N) 24 : IMPLICIT REAL*8 (A-H,O-Z) 25 : 26 : C*********************************************************************** 27 : C 28 : C CE SP DISCRETISE UNE EQUATION GENERALE DE TRANSPORT-DIFFUSION AVEC 29 : C SOURCE. 30 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI 31 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 32 : C LES OPERATEURS SONT "SOUS-INTEGRES" 33 : C 34 : C 35 : C APPELE PAR YTSCAL 36 : C 37 : C 38 : C*********************************************************************** 39 : 40 : -INC CCVQUA4 41 : -INC CCREEL 42 : -INC SMCOORD 43 : -INC CCOPTIO 44 : 45 : C 46 : C Longueur des registres vectoriels de la machine cible 47 : C On prend 64 pour ne pas augmenter la taille des tableaux 48 : C nécessaires à la vectorisation. 49 : C 50 : PARAMETER(LRV=64) 51 : 52 : DIMENSION UN(NPTU,IES),HRN(NPTD),TN(NPTD) 53 : DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*) 54 : DIMENSION ALFE(*),ALT(*) 55 : 56 : DIMENSION IPADL(1),LE(NP,1),IPADU(*) 57 : DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL) 58 : 59 : DIMENSION BF(9,9) 60 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 61 : 62 : DIMENSION COEF(LRV),AIRE(LRV) 63 : DIMENSION AL2(LRV),AH2(LRV),AP2(LRV) 64 : DIMENSION AL(LRV),AH(LRV),AP(LRV) 65 : DIMENSION ALFT(LRV),QT(LRV) 66 : DIMENSION XMB(LRV),XMH(LRV) 67 : DIMENSION CFM(LRV) 68 : C DIMENSION CF1(LRV),CF2(LRV),CF3(LRV) 69 : DIMENSION DR(LRV,9) 70 : DIMENSION UM(LRV),UP(LRV) 71 : DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9) 72 : DIMENSION TETAC(LRV,9),TETAD(LRV,9),TETA(LRV,9) 73 : DIMENSION UMI(LRV,3),UPI(LRV,3) 74 : DIMENSION CXT(LRV),CYT(LRV),CXY(LRV) 75 : DIMENSION DXT(LRV),DYT(LRV),DXY(LRV) 76 : DIMENSION CZT(LRV),CXZ(LRV),CYZ(LRV) 77 : DIMENSION DZT(LRV),DXZ(LRV),DYZ(LRV) 78 : DIMENSION SBF(LRV,9) 79 : DIMENSION GRADT(LRV,3) 80 : DIMENSION BM(LRV),BP(LRV) 81 : 82 : DIMENSION BMX(LRV),BMY(LRV),BPX(LRV),BPY(LRV) 83 : 84 : REAL*8 G(1),n(2,4),MMAX(4),lll,Fi 85 : REAL*8 b(2),Nm,x(4),y(4),kkk(4),Kchap 86 : INTEGER zz,p,kdesign 87 : PARAMETER (zz=1) 88 : 89 : SAVE IPAS,QGGT,Q1,Q2,Q3 90 : 91 : DATA CD/1.D0/ 92 : 93 : DATA IPAS/0/ 94 : 95 : DATA IDCENN/2/ 96 : C************************************************************************ 97 : C 98 : C INITIALISATIONS DIVERSES 99 : C 100 : KDESIGN=0 101 : NK=K0 102 : C ******** 103 : C * 2D * 104 : C ******** 105 : 106 : 107 : IF (NP.EQ.3) THEN 108 : x(4)=0.d0 109 : y(4)=0.d0 110 : ENDIF 111 : 112 : IF(IES.EQ.3)GO TO 10 113 : 114 : IAX1=0 115 : IAX2=0 116 : IF(IAXI.EQ.1)IAX2=1 117 : IF(IAXI.EQ.2)IAX1=1 118 : 119 : HMIN = 1.D20 120 : HMAX = 0.D0 121 : 122 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 123 : 124 : C DIFFERENCES TRIANGLE / QUADRANGLE 125 : QUA4=0.D0 126 : IF(NP.EQ.4)QUA4=1.D0 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 : DO 7001 KPACK=KPACKD,KPACKF 142 : C 143 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 144 : C 145 : C 1. Calcul des limites du paquet courant. 146 : KDEB=1+(KPACK-1)*LRV 147 : KFIN=MIN(NEL,KDEB+LRV-1) 148 : C 149 : DO 7002 K=KDEB,KFIN 150 : KP=K-KDEB+1 151 : NK=K+K0 152 : NK1=(1-IND1)*(NK-1)+1 153 : ALFT(KP)=ALFE(NK1)+XPETIT 154 : AIRE(KP)=VOLU(NK) 155 : AL(KP)=COTE(NK,1) 156 : AH(KP)=COTE(NK,2) 157 : AL2(KP)=1.D0/AL(KP)/AL(KP) 158 : AH2(KP)=1.D0/AH(KP)/AH(KP) 159 : XMH(KP)=(AL(KP)+AH(KP))/2.D0 160 : 7002 CONTINUE 161 : 162 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 163 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 164 : DO 7003 K=KDEB,KFIN 165 : KP=K-KDEB+1 166 : NK=K+K0 167 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 168 : 7003 CONTINUE 169 : ENDIF 170 : 171 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 172 : 173 : C Initialisation des UMI avant accumulation 174 : DO 7005 K=KDEB,KFIN 175 : KP=K-KDEB+1 176 : UMI(KP,1)=XPETIT 177 : UMI(KP,2)=XPETIT 178 : UPI(KP,1)=XPETIT 179 : UPI(KP,2)=XPETIT 180 : GRADT(KP,1)=XPETIT 181 : GRADT(KP,2)=XPETIT 182 : 7005 CONTINUE 183 : 184 : DO 7006 I=1,NP 185 : C*IBMDIR* PREFER VECTOR 186 : DO 7006 K=KDEB,KFIN 187 : KP=K-KDEB+1 188 : NF=IPADL(LE(I,K)) 189 : NFU=IPADU(LE(I,K)) 190 : NFU=(1-INDU)*(NFU-1)+1 191 : DR(KP,I)=DRR(I,K) 192 : UIX(KP,I)=UN(NFU,1) 193 : UIY(KP,I)=UN(NFU,2) 194 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 195 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 196 : 197 : TETAC(KP,I)=HRN(NF) 198 : TETAD(KP,I)=TN(NF) 199 : GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I) 200 : GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I) 201 : 7006 CONTINUE 202 : 203 : C 204 : C Initialisation de la variable d'accumulation SBF au terme source 205 : C 206 : C write(6,*)' IKomp,ias=',IKomp,ikas 207 : C write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU 208 : 209 : IF(IKOMP.EQ.0)THEN 210 : 211 : DO 70021 K=KDEB,KFIN 212 : KP=K-KDEB+1 213 : NK=K+K0 214 : NKS=(1-IKS)*(NK-1)+1 215 : QT(KP)=QE(NKS) 216 : 70021 CONTINUE 217 : 218 : DO 70062 I=1,NP 219 : C*IBMDIR* PREFER VECTOR 220 : DO 70062 K=KDEB,KFIN 221 : KP=K-KDEB+1 222 : SBF(KP,I)=-QT(KP)*DR(KP,I) 223 : 70062 CONTINUE 224 : 225 : ELSEIF(IKOMP.EQ.1)THEN 226 : 227 : DO 70023 K=KDEB,KFIN 228 : KP=K-KDEB+1 229 : NK=K+K0 230 : NKS=(1-IKS)*(NK-1)+1 231 : QT(KP)=QE(NKS) 232 : 70023 CONTINUE 233 : 234 : DO 70064 I=1,NP 235 : C*IBMDIR* PREFER VECTOR 236 : DO 70064 K=KDEB,KFIN 237 : KP=K-KDEB+1 238 : SBF(KP,I)=-QT(KP)*DR(KP,I) 239 : 70064 CONTINUE 240 : 241 : ENDIF 242 : 243 : DO 7007 K=KDEB,KFIN 244 : KP=K-KDEB+1 245 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 246 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 247 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 248 : UM(KP)=SQRT(UM(KP)) + XPETIT 249 : 250 : A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2) 251 : A=SQRT(A)+XPETIT 252 : GRADT(KP,1)=GRADT(KP,1)/A 253 : GRADT(KP,2)=GRADT(KP,2)/A 254 : UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2) 255 : UPI(KP,1)=UP(KP)*GRADT(KP,1) 256 : UPI(KP,2)=UP(KP)*GRADT(KP,2) 257 : UP(KP) = ABS(UP(KP)) + XPETIT 258 : 7007 CONTINUE 259 : 260 : * DO 70074 K=KDEB,KFIN 261 : * KP=K-KDEB+1 262 : * BMX=UMI(KP,1)/AL(KP) 263 : * BMY=UMI(KP,2)/AH(KP) 264 : * BM(KP)=BMX*BMX+BMY*BMY 265 : * BM(KP)=SQRT(BM(KP))+XPETIT 266 : * BPX=UPI(KP,1)/AL(KP) 267 : * BPY=UPI(KP,2)/AH(KP) 268 : * BP(KP)=BPX*BPX+BPY*BPY 269 : * BP(KP)=SQRT(BP(KP))+XPETIT 270 : *70074 CONTINUE 271 : 272 : 273 : DO 70074 K=KDEB,KFIN 274 : KP=K-KDEB+1 275 : BMX(KP)=UMI(KP,1)/AL(KP) 276 : BMY(KP)=UMI(KP,2)/AH(KP) 277 : 278 : BM(KP)=BMX(KP)*BMX(KP)+BMY(KP)*BMY(KP) 279 : BM(KP)=SQRT(BM(KP))+XPETIT 280 : BPX(KP)=UPI(KP,1)/AL(KP) 281 : BPY(KP)=UPI(KP,2)/AH(KP) 282 : BP(KP)=BPX(KP)*BPX(KP)+BPY(KP)*BPY(KP) 283 : BP(KP)=SQRT(BP(KP))+XPETIT 284 : 70074 CONTINUE 285 : 286 : 287 : C 288 : C AVANT DECENTREMENT 289 : C 290 : 291 : 292 : IF(IKOMP.EQ.0)THEN 293 : 294 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008 295 : IF(IDCENN.EQ.1)THEN 296 : DO 70081 K=KDEB,KFIN 297 : KP=K-KDEB+1 298 : XMB(KP)=XMH(KP) 299 : CXT(KP)=0.D0 300 : CYT(KP)=0.D0 301 : CXY(KP)=0.D0 302 : 70081 CONTINUE 303 : 304 : ELSEIF(IDCENN.EQ.2)THEN 305 : 306 : DO 7008 K=KDEB,KFIN 307 : KP=K-KDEB+1 308 : 309 : HMK=2.D0*UM(KP)/BM(KP) 310 : HPK=2.D0*UP(KP)/BP(KP) 311 : 312 : IF (HMK.LT.HMIN) HMIN=HMK 313 : IF (HMK.GT.HMAX) HMAX=HMK 314 : 315 : XMB(KP)=HMK 316 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 317 : AKSI=ALFA/3.D0 318 : IF(ALFA.GT.3.D0)AKSI=1.D0 319 : CCT=AKSI*HMK/2.d0/UM(KP) 320 : 321 : Fi= (GRADT(KP,1)*UMI(KP,1) 322 : &+ GRADT(KP,2)*UMI(KP,2)) 323 : 324 : Kchap=0.5*(XMB(KP)*abs(Fi))/(sqrt((GRADT(KP,1)**2+ 325 : & GRADT(KP,2)**2))+XMB(KP)) 326 : 327 : CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+Kchap) 328 : CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+Kchap) 329 : CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT) 330 : 331 : 7008 CONTINUE 332 : 333 : ELSEIF(IDCENN.EQ.3)THEN 334 : DO 7018 K=KDEB,KFIN 335 : KP=K-KDEB+1 336 : C DESIGN TEZDUYAR 337 : 338 : IF (KDESIGN.EQ.1) THEN 339 : * IF (NP.NE.3) RETURN 340 : S=0.d0 341 : do 119 I=1,NP 342 : 119 S=S+abs(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))/UM(KP) 343 : HMK=2.d0/S 344 : C DESIGN ORIGINAL 345 : ELSEIF (KDESIGN.EQ.0) THEN 346 : HMK=2.D0*UM(KP)/BM(KP) 347 : 348 : 349 : ELSEIF (KDESIGN.EQ.2) THEN 350 : do 774 i=1,NP 351 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 352 : 774 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 353 : 354 : DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3)) 355 : * -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3)) 356 : 357 : IF (NP.EQ.3) lll=1.d0 358 : IF (NP.EQ.4) lll=4.d0 359 : DJ=DJ/(lll**2) 360 : 361 : b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)- 362 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 363 : b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)- 364 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 365 : 366 : p=1 367 : s=0.d0 368 : do 304 kk=1,2 369 : 304 s=s+(abs(b(kk)))**p 370 : Nm=s**(1.d0/p) 371 : HMK=2.d0*UM(KP)/Nm 372 : 373 : ELSEIF (KDESIGN.EQ.3) THEN 374 : do 775 i=1,NP 375 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 376 : 775 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 377 : 378 : DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3)) 379 : * -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3)) 380 : 381 : IF (NP.EQ.3) lll=1.d0 382 : IF (NP.EQ.4) lll=4.d0 383 : DJ=DJ/(lll**2) 384 : 385 : b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)- 386 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 387 : b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)- 388 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 389 : 390 : p=2 391 : s=0.d0 392 : do 305 kk=1,2 393 : 305 s=s+(abs(b(kk)))**p 394 : Nm=s**(1.d0/p) 395 : HMK=2.d0*UM(KP)/Nm 396 : 397 : ELSEIF (KDESIGN.EQ.4) THEN 398 : do 776 i=1,NP 399 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 400 : 776 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 401 : 402 : do 318 i=1,NP 403 : 318 kkk(i)=Aire(KP)*(UMI(KP,1)*HR(K,i,1)+UMI(KP,2)*HR(K,i,2)) 404 : 405 : s0=0.d0 406 : s1=0.d0 407 : do 319 i=1,NP 408 : s0=s0+min(0.d0,kkk(i)) 409 : 319 s1=s1+max(0.d0,kkk(i)) 410 : 411 : s2=0.d0 412 : s3=0.d0 413 : do 320 i=1,NP 414 : s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1 415 : 320 s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1 416 : HMK=sqrt(s2**2+s3**2) 417 : 418 : ELSEIF (KDESIGN.EQ.5) THEN 419 : do 800 i=1,NP 420 : n(1,i)=abs(2.d0*Aire(KP)*HR(K,i,1)) 421 : n(2,i)=abs(2.d0*Aire(KP)*HR(K,i,2)) 422 : MMAX(i)=n(1,i) 423 : 800 IF (n(1,i).lt.n(2,i)) MMAX(i)=n(2,i) 424 : 425 : HMK=0.d0 426 : do 801 i=1,NP 427 : 801 IF (MMAX(i).gt.HMK) HMK=MMAX(i) 428 : 429 : ENDIF 430 : 431 : 432 : IF (HMK.LT.HMIN) HMIN=HMK 433 : IF (HMK.GT.HMAX) HMAX=HMK 434 : 435 : XMB(KP)=HMK 436 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 437 : AKSI=ALFA/3.D0 438 : IF(ALFA.GT.3.D0)AKSI=1.D0 439 : CCT=AKSI*HMK/2.d0/UM(KP) 440 : 441 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 442 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 443 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 444 : 445 : 7018 CONTINUE 446 : 447 : ELSEIF(IDCENN.EQ.4)THEN 448 : DT19=DTM1*0.5D0 449 : DO 7009 K=KDEB,KFIN 450 : KP=K-KDEB+1 451 : XMB(KP)=XMH(KP) 452 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 453 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 454 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 455 : 7009 CONTINUE 456 : 457 : 458 : ENDIF 459 : ELSEIF(IKOMP.EQ.1)THEN 460 : 461 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008 462 : IF(IDCENN.EQ.1)THEN 463 : DO 71081 K=KDEB,KFIN 464 : KP=K-KDEB+1 465 : XMB(KP)=XMH(KP) 466 : 467 : CXT(KP)=0.D0 468 : CYT(KP)=0.D0 469 : CXY(KP)=0.D0 470 : 471 : DXT(KP)=0.D0 472 : DYT(KP)=0.D0 473 : DXY(KP)=0.D0 474 : 475 : 71081 CONTINUE 476 : 477 : ELSEIF(IDCENN.EQ.2)THEN 478 : DO 7108 K=KDEB,KFIN 479 : KP=K-KDEB+1 480 : C DESIGN TEZDUYAR 481 : IF (KDESIGN.EQ.1) THEN 482 : * IF (NP.NE.3) RETURN 483 : S=0.d0 484 : do 120 I=1,NP 485 : 120 S=S+abs(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))/UM(KP) 486 : HMK=2.d0/S 487 : S=0.d0 488 : do 121 I=1,NP 489 : 121 S=S+abs(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2))/UP(KP) 490 : HPK=2.d0/S 491 : C DESIGN ORIGINAL 492 : ELSEIF (KDESIGN.EQ.0) THEN 493 : HMK=2.D0*UM(KP)/BM(KP) 494 : HPK=2.D0*UP(KP)/BP(KP) 495 : 496 : ELSEIF (KDESIGN.EQ.2) THEN 497 : do 777 i=1,NP 498 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 499 : 777 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 500 : 501 : DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3)) 502 : * -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3)) 503 : 504 : IF (NP.EQ.3) lll=1.d0 505 : IF (NP.EQ.4) lll=4.d0 506 : DJ=DJ/(lll**2) 507 : 508 : b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)- 509 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 510 : b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)- 511 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 512 : 513 : p=1 514 : s=0.d0 515 : do 306 kk=1,2 516 : 306 s=s+(abs(b(kk)))**p 517 : Nm=s**(1.d0/p) 518 : HMK=2.d0*UM(KP)/Nm 519 : C Calcul de HPK: 520 : 521 : b(1)=(UPI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UPI(KP,2)*(x(3)- 522 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 523 : b(2)=(-UPI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UPI(KP,2)*(x(2)- 524 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 525 : 526 : s=0.d0 527 : do 308 kk=1,2 528 : 308 s=s+(abs(b(kk)))**p 529 : Nm=s**(1.d0/p) 530 : HPK=2.d0*UP(KP)/Nm 531 : 532 : ELSEIF (KDESIGN.EQ.3) THEN 533 : do 778 i=1,NP 534 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 535 : 778 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 536 : 537 : DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3)) 538 : * -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3)) 539 : 540 : IF (NP.EQ.3) lll=1.d0 541 : IF (NP.EQ.4) lll=4.d0 542 : DJ=DJ/(lll**2) 543 : 544 : b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)- 545 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 546 : b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)- 547 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 548 : 549 : p=2 550 : s=0.d0 551 : do 309 kk=1,2 552 : 309 s=s+(abs(b(kk)))**p 553 : Nm=s**(1.d0/p) 554 : HMK=2.d0*UM(KP)/Nm 555 : C Calcul de HPK: 556 : 557 : b(1)=(UPI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UPI(KP,2)*(x(3)- 558 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 559 : b(2)=(-UPI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UPI(KP,2)*(x(2)- 560 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 561 : 562 : s=0.d0 563 : do 310 kk=1,2 564 : 310 s=s+(abs(b(kk)))**p 565 : Nm=s**(1.d0/p) 566 : HPK=2.d0*UP(KP)/Nm 567 : 568 : ELSEIF (KDESIGN.EQ.4) THEN 569 : do 779 i=1,NP 570 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 571 : 779 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 572 : 573 : * ZZ=1 :HMK=HPK 574 : * ZZ=2 :HMK<>HPK 575 : 576 : do 321 i=1,NP 577 : 321 kkk(i)=Aire(KP)*(UMI(KP,1)*HR(K,i,1)+UMI(KP,2)*HR(K,i,2)) 578 : 579 : s0=0.d0 580 : s1=0.d0 581 : do 322 i=1,NP 582 : s0=s0+min(0.d0,kkk(i)) 583 : 322 s1=s1+max(0.d0,kkk(i)) 584 : 585 : s2=0.d0 586 : s3=0.d0 587 : do 323 i=1,NP 588 : s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1 589 : 323 s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1 590 : HMK=sqrt(s2**2+s3**2) 591 : HPK=HMK 592 : 593 : IF (zz.eq.2) THEN 594 : do 324 i=1,NP 595 : 324 kkk(i)=Aire(KP)*(UPI(KP,1)*HR(K,i,1)+UPI(KP,2)*HR(K,i,2)) 596 : 597 : s0=0.d0 598 : s1=0.d0 599 : do 325 i=1,NP 600 : s0=s0+min(0.d0,kkk(i)) 601 : 325 s1=s1+max(0.d0,kkk(i)) 602 : 603 : s2=0.d0 604 : s3=0.d0 605 : do 3266 i=1,NP 606 : s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1 607 : 3266 s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1 608 : HPK=sqrt(s2**2+s3**2) 609 : ENDIF 610 : 611 : ELSEIF (KDESIGN.EQ.5) THEN 612 : do 802 i=1,NP 613 : n(1,i)=abs(2.d0*Aire(KP)*HR(K,i,1)) 614 : n(2,i)=abs(2.d0*Aire(KP)*HR(K,i,2)) 615 : MMAX(i)=n(1,i) 616 : 802 IF (n(1,i).lt.n(2,i)) MMAX(i)=n(2,i) 617 : 618 : HMK=0.d0 619 : do 803 i=1,NP 620 : 803 IF (MMAX(i).gt.HMK) HMK=MMAX(i) 621 : 622 : ENDIF 623 : 624 : IF (HMK.LT.HMIN) HMIN=HMK 625 : IF (HMK.GT.HMAX) HMAX=HMK 626 : 627 : XMB(KP)=HMK 628 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 629 : AKSI=ALFA/3.D0 630 : IF(ALFA.GT.3.D0)AKSI=1.D0 631 : CCT=AKSI*HMK/2.d0/UM(KP) 632 : 633 : * si on a le meme h 634 : HPK=HMK 635 : ALFA=UM(KP)*HPK/ALFT(KP)/2.D0 636 : AKSI=ALFA/3.D0 637 : IF(ALFA.GT.3.D0)AKSI=1.D0 638 : CCP=AKSI*HPK/2.d0/UP(KP) 639 : CPT=CCP-CCT 640 : CC2=0.D0 641 : IF(CPT.GE.0.D0)CC2=CPT 642 : 643 : CXT(KP)=UMI(KP,1)*CCT 644 : CYT(KP)=UMI(KP,2)*CCT 645 : CXY(KP)=CCT 646 : 647 : DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2 648 : DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2 649 : DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2 650 : 651 : 7108 CONTINUE 652 : 653 : ELSEIF(IDCENN.EQ.3)THEN 654 : DO 71018 K=KDEB,KFIN 655 : KP=K-KDEB+1 656 : C DESIGN TEZDUYAR 657 : IF (KDESIGN.EQ.1) THEN 658 : * IF (NP.NE.3) RETURN 659 : S=0.d0 660 : do 122 I=1,NP 661 : 122 S=S+abs(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))/UM(KP) 662 : HMK=2.d0/S 663 : C DESIGN ORIGINAL 664 : ELSEIF (KDESIGN.EQ.0) THEN 665 : HMK=2.D0*UM(KP)/BM(KP) 666 : 667 : ELSEIF (KDESIGN.EQ.2) THEN 668 : do 780 i=1,NP 669 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 670 : 780 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 671 : 672 : DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3)) 673 : * -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3)) 674 : 675 : IF (NP.EQ.3) lll=1.d0 676 : IF (NP.EQ.4) lll=4.d0 677 : DJ=DJ/(lll**2) 678 : 679 : b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)- 680 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 681 : b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)- 682 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 683 : 684 : p=1 685 : s=0.d0 686 : do 311 kk=1,2 687 : 311 s=s+(abs(b(kk)))**p 688 : Nm=s**(1.d0/p) 689 : HMK=2.d0*UM(KP)/Nm 690 : 691 : 692 : ELSEIF (KDESIGN.EQ.3) THEN 693 : do 781 i=1,NP 694 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 695 : 781 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 696 : 697 : DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3)) 698 : * -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3)) 699 : 700 : IF (NP.EQ.3) lll=1.d0 701 : IF (NP.EQ.4) lll=4.d0 702 : DJ=DJ/(lll**2) 703 : 704 : b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)- 705 : * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll 706 : b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)- 707 : * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll 708 : 709 : p=2 710 : s=0.d0 711 : do 312 kk=1,2 712 : 312 s=s+(abs(b(kk)))**p 713 : Nm=s**(1.d0/p) 714 : HMK=2.d0*UM(KP)/Nm 715 : 716 : ELSEIF (KDESIGN.EQ.4) THEN 717 : do 782 i=1,NP 718 : x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1) 719 : 782 y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2) 720 : 721 : do 326 i=1,NP 722 : 326 kkk(i)=Aire(KP)*(UMI(KP,1)*HR(K,i,1)+UMI(KP,2)*HR(K,i,2)) 723 : 724 : s0=0.d0 725 : s1=0.d0 726 : do 3277 i=1,KP 727 : s0=s0+min(0.d0,kkk(i)) 728 : 3277 s1=s1+max(0.d0,kkk(i)) 729 : 730 : s2=0.d0 731 : s3=0.d0 732 : do 3288 i=1,KP 733 : s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1 734 : 3288 s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1 735 : HMK=sqrt(s2**2+s3**2) 736 : 737 : ELSEIF (KDESIGN.EQ.5) THEN 738 : do 804 i=1,NP 739 : n(1,i)=abs(2.d0*Aire(KP)*HR(K,i,1)) 740 : n(2,i)=abs(2.d0*Aire(KP)*HR(K,i,2)) 741 : MMAX(i)=n(1,i) 742 : 804 IF (n(1,i).lt.n(2,i)) MMAX(i)=n(2,i) 743 : 744 : HMK=0.d0 745 : do 805 i=1,NP 746 : 805 IF (MMAX(i).gt.HMK) HMK=MMAX(i) 747 : 748 : ENDIF 749 : 750 : IF (HMK.LT.HMIN) HMIN=HMK 751 : IF (HMK.GT.HMAX) HMAX=HMK 752 : 753 : XMB(KP)=HMK 754 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 755 : AKSI=ALFA/3.D0 756 : IF(ALFA.GT.3.D0)AKSI=1.D0 757 : CCT=AKSI*HMK/2.d0/UM(KP) 758 : 759 : CXT(KP)=UMI(KP,1)*CCT 760 : CYT(KP)=UMI(KP,2)*CCT 761 : CXY(KP)=CCT 762 : 763 : DXT(KP)=0.D0 764 : DYT(KP)=0.D0 765 : DXY(KP)=0.D0 766 : 767 : 71018 CONTINUE 768 : 769 : ELSEIF(IDCENN.EQ.4)THEN 770 : DT19=DTM1*0.5D0 771 : DO 71009 K=KDEB,KFIN 772 : KP=K-KDEB+1 773 : XMB(KP)=XMH(KP) 774 : CXT(KP)=UMI(KP,1)*DT19 775 : CYT(KP)=UMI(KP,2)*DT19 776 : CXY(KP)=DT19 777 : 778 : DXT(KP)=0.D0 779 : DYT(KP)=0.D0 780 : DXY(KP)=0.D0 781 : 782 : 71009 CONTINUE 783 : 784 : ENDIF 785 : 786 : ENDIF 787 : C *********************** 788 : C 789 : C AVANT CALCUL DT 790 : C 791 : IF(IKOMP.EQ.0)THEN 792 : 793 : C*IBMDIR* PREFER SCALAR 794 : DO 7010 K=KDEB,KFIN 795 : KP=K-KDEB+1 796 : DT0=DT 797 : 798 : DT1=2.D0/ 799 : & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)) 800 : & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)) ) 801 : 802 : DT2=0.5D0/ 803 : & ( (ALFT(KP)+CXT(KP))*AL2(KP) 804 : & +(ALFT(KP)+CYT(KP))*AH2(KP) ) 805 : 806 : 807 : IF(DT1.LT.DT)DT=DT1 808 : IF(DT2.LT.DT)DT=DT2 809 : C IF(DT3.LT.DT)DT=DT3 810 : IF(DT.NE.DT0) THEN 811 : DTT1=DT1 812 : DTT2=DT2 813 : DTT3=0.D0 814 : DIAEL=XMB(KP) 815 : NUEL=K 816 : END IF 817 : CXY(KP)=CXY(KP)*CD 818 : 7010 CONTINUE 819 : 820 : ELSEIF(IKOMP.EQ.1)THEN 821 : 822 : C*IBMDIR* PREFER SCALAR 823 : DO 7110 K=KDEB,KFIN 824 : KP=K-KDEB+1 825 : DT0=DT 826 : 827 : DT1=2.D0/ 828 : & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP)) 829 : & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP)) ) 830 : 831 : DT2=0.5D0/ 832 : & ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP) 833 : & +(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP) ) 834 : 835 : IF(DT1.LT.DT)DT=DT1 836 : IF(DT2.LT.DT)DT=DT2 837 : C IF(DT3.LT.DT)DT=DT3 838 : IF(DT.NE.DT0) THEN 839 : DTT1=DT1 840 : DTT2=DT2 841 : DTT3=0.D0 842 : DIAEL=XMB(KP) 843 : NUEL=K 844 : END IF 845 : CXY(KP)=CXY(KP)*CD 846 : 7110 CONTINUE 847 : 848 : ENDIF 849 : C Le coeur du calcul ... 850 : 851 : IF(IKOMP.EQ.0)THEN 852 : 853 : DO 7014 I=1,NP 854 : DO 7014 J= 1,NP 855 : C*IBMDIR* PREFER VECTOR 856 : 857 : DO 7014 K=KDEB,KFIN 858 : KP=K-KDEB+1 859 : 860 : ZVGG=AIRE(KP)*( 861 : & HR(K,I,1)*HR(K,J,1)*CXT(KP) 862 : &+ HR(K,I,1)*HR(K,J,2)*CXY(KP) 863 : &+ HR(K,I,2)*HR(K,J,1)*CXY(KP) 864 : &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) 865 : &+ (CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0 866 : & *VGGT(J,I)*QUA4 ) 867 : 868 : C &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) ) 869 : C &+ (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4 870 : 871 : ZVGT=AIRE(KP)*( 872 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 873 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 874 : &+ (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0 875 : & *VGGT(J,I)*QUA4 ) 876 : 877 : C &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) ) 878 : C &+ ALFT(KP)/12.D0*VGGT(J,I)*QUA4 879 : 880 : C? V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I) 881 : V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)*HR(K,J,2) 882 : 883 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 884 : 885 : 7014 CONTINUE 886 : ELSEIF(IKOMP.EQ.1)THEN 887 : 888 : DO 7015 I=1,NP 889 : DO 7015 J= 1,NP 890 : C*IBMDIR* PREFER VECTOR 891 : DO 7015 K=KDEB,KFIN 892 : KP=K-KDEB+1 893 : ZVGG=AIRE(KP)*( 894 : & HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP)) 895 : &+ HR(K,I,1)*HR(K,J,2) 896 : & *(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP)) 897 : &+ HR(K,I,2)*HR(K,J,1) 898 : & *(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP)) 899 : &+ HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP)) 900 : &+ ((CXT(KP)*UIX(KP,J)+DXT(KP))*AL2(KP) 901 : & +(CYT(KP)*UIY(KP,J)+DYT(KP))*AH2(KP))/12.D0 902 : & *VGGT(J,I)*QUA4 ) 903 : 904 : C &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) ) 905 : C &+ (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4 906 : 907 : ZVGT=AIRE(KP)*( 908 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 909 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 910 : &+ (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0 911 : & *VGGT(J,I)*QUA4 ) 912 : 913 : C &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) ) 914 : C &+ ALFT(KP)/12.D0*VGGT(J,I)*QUA4 915 : 916 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2))*DR(KP,I) 917 : 918 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 919 : 920 : 7015 CONTINUE 921 : 922 : ENDIF 923 : C 924 : C Fin de l'accumulation dans SBF. 925 : C On ajoute ces incréments G. 926 : C 927 : DO 7017 I=1,NP 928 : C*IBMDIR* PREFER VECTOR 929 : DO 7017 K=KDEB,KFIN 930 : KP=K-KDEB+1 931 : NF=IPADL(LE(I,K)) 932 : G(NF) = G(NF)+SBF(KP,I) 933 : 7017 CONTINUE 934 : 935 : 7001 CONTINUE 936 : 937 : C WRITE(6,*)' G DANS YCTSCL ' 938 : C WRITE(6,1984)(M,G(M),M=1,NPTD) 939 : 1984 FORMAT(7(1X,I4,2X,1PE11.4)) 940 : 940 : C CALL ARRET(0) 941 : IPAS=1 942 : C 943 : C 944 : C PRINT *,'HMIN = ', HMIN, 'HMAX = ', HMAX 945 : C 946 : RETURN 947 : 948 : C ******** 949 : C * 3D * 950 : C ******** 951 : 952 : 10 CONTINUE 953 : 953 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 954 : CUB8=0.D0 955 : IF(NP.EQ.8)CUB8=1.D0 956 : 957 : C 958 : C Calcul du nombre de paquets de LRV éléments 959 : C 960 : NNN=MOD(NEL,LRV) 961 : IF(NNN.EQ.0) NPACK=NEL/LRV 962 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 963 : KPACKD=1 964 : KPACKF=NPACK 965 : C 966 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 967 : C 968 : DO 8001 KPACK=KPACKD,KPACKF 969 : C 970 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 971 : C 972 : C 1. Calcul des limites du paquet courant. 973 : KDEB=1+(KPACK-1)*LRV 974 : KFIN=MIN(NEL,KDEB+LRV-1) 975 : C 976 : DO 8002 K=KDEB,KFIN 977 : KP=K-KDEB+1 978 : NK=K+K0 979 : NK1=(1-IND1)*(NK-1)+1 980 : ALFT(KP)=ALFE(NK1)+XPETIT 981 : AIRE(KP)=VOLU(NK) 982 : AL(KP)=COTE(NK,1) 983 : AH(KP)=COTE(NK,2) 984 : AP(KP)=COTE(NK,3) 985 : CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ 986 : & AP(KP)*AH(KP)/AL(KP) 987 : C CF1(KP)=AL(KP)*AH(KP)/AP(KP) 988 : C CF2(KP)=AL(KP)*AP(KP)/AH(KP) 989 : C CF3(KP)=AP(KP)*AH(KP)/AL(KP) 990 : XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0 991 : AL2(KP)=1.D0/AL(KP)/AL(KP) 992 : AH2(KP)=1.D0/AH(KP)/AH(KP) 993 : AP2(KP)=1.D0/AP(KP)/AP(KP) 994 : 8002 CONTINUE 995 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 996 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 997 : DO 8003 K=KDEB,KFIN 998 : KP=K-KDEB+1 999 : NK=K+K0 1000 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 1001 : 8003 CONTINUE 1002 : ENDIF 1003 : 1004 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1005 : 1006 : C Initialisation des UMI avant accumulation 1007 : DO 8005 K=KDEB,KFIN 1008 : KP=K-KDEB+1 1009 : UMI(KP,1)=XPETIT 1010 : UMI(KP,2)=XPETIT 1011 : UMI(KP,3)=XPETIT 1012 : UPI(KP,1)=XPETIT 1013 : UPI(KP,2)=XPETIT 1014 : UPI(KP,3)=XPETIT 1015 : GRADT(KP,1)=XPETIT 1016 : GRADT(KP,2)=XPETIT 1017 : GRADT(KP,3)=XPETIT 1018 : 8005 CONTINUE 1019 : 1020 : DO 8006 I=1,NP 1021 : C*IBMDIR* PREFER VECTOR 1022 : DO 8006 K=KDEB,KFIN 1023 : KP=K-KDEB+1 1024 : NF=IPADL(LE(I,K)) 1025 : NFU=IPADU(LE(I,K)) 1026 : NFU=(1-INDU)*(NFU-1)+1 1027 : DR(KP,I)=DRR(I,K) 1028 : UIX(KP,I)=UN(NFU,1) 1029 : UIY(KP,I)=UN(NFU,2) 1030 : UIZ(KP,I)=UN(NFU,3) 1031 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 1032 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 1033 : UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I) 1034 : 1035 : TETAC(KP,I)=HRN(NF) 1036 : TETAD(KP,I)=TN(NF) 1037 : GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I) 1038 : GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I) 1039 : GRADT(KP,3)=GRADT(KP,3)+HR(K,I,3)*TETAC(KP,I) 1040 : 8006 CONTINUE 1041 : 1042 : C 1043 : C Initialisation de la variable d'accumulation SBF au terme source 1044 : C M 1045 : C write(6,*)' IKomp,ikas=',IKomp,ikas 1046 : C write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU 1047 : 1048 : IF(IKOMP.EQ.0)THEN 1049 : 1050 : DO 80021 K=KDEB,KFIN 1051 : KP=K-KDEB+1 1052 : NK=K+K0 1053 : NKS=(1-IKS)*(NK-1)+1 1054 : QT(KP)=QE(NKS) 1055 : 80021 CONTINUE 1056 : 1057 : DO 80062 I=1,NP 1058 : C*IBMDIR* PREFER VECTOR 1059 : DO 80062 K=KDEB,KFIN 1060 : KP=K-KDEB+1 1061 : SBF(KP,I)=-QT(KP)*DR(KP,I) 1062 : 80062 CONTINUE 1063 : 1064 : ELSEIF(IKOMP.EQ.1)THEN 1065 : 1066 : DO 80023 K=KDEB,KFIN 1067 : KP=K-KDEB+1 1068 : NK=K+K0 1069 : NKS=(1-IKS)*(NK-1)+1 1070 : QT(KP)=QE(NKS) 1071 : 80023 CONTINUE 1072 : 1073 : DO 80064 I=1,NP 1074 : C*IBMDIR* PREFER VECTOR 1075 : DO 80064 K=KDEB,KFIN 1076 : KP=K-KDEB+1 1077 : SBF(KP,I)=-QT(KP)*DR(KP,I) 1078 : 80064 CONTINUE 1079 : 1080 : ENDIF 1081 : 1082 : DO 8007 K=KDEB,KFIN 1083 : KP=K-KDEB+1 1084 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 1085 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 1086 : UMI(KP,3)=UMI(KP,3)/AIRE(KP) 1087 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 1088 : & +UMI(KP,3)*UMI(KP,3) 1089 : UM(KP)=SQRT(UM(KP)) 1090 : A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2) 1091 : & +GRADT(KP,3)*GRADT(KP,3) 1092 : A=SQRT(A)+XPETIT 1093 : GRADT(KP,1)=GRADT(KP,1)/A 1094 : GRADT(KP,2)=GRADT(KP,2)/A 1095 : GRADT(KP,3)=GRADT(KP,3)/A 1096 : UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2) 1097 : & +GRADT(KP,3)*UMI(KP,3) 1098 : UPI(KP,1)=UP(KP)*GRADT(KP,1) 1099 : UPI(KP,2)=UP(KP)*GRADT(KP,2) 1100 : UPI(KP,3)=UP(KP)*GRADT(KP,3) 1101 : 8007 CONTINUE 1102 : 1103 : * DO 80074 K=KDEB,KFIN 1104 : * KP=K-KDEB+1 1105 : * BMX=UMI(KP,1)/AL(KP) 1106 : * BMY=UMI(KP,2)/AH(KP) 1107 : * BMZ=UMI(KP,3)/AP(KP) 1108 : * BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ 1109 : * BM(KP)=SQRT(BM(KP))+XPETIT 1110 : * BPX=UPI(KP,1)/AL(KP) 1111 : * BPY=UPI(KP,2)/AH(KP) 1112 : * BPZ=UPI(KP,3)/AP(KP) 1113 : * BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ 1114 : * BP(KP)=SQRT(BP(KP))+XPETIT 1115 : *80074 CONTINUE 1116 : 1117 : C 1118 : C AVANT DECENTREMENT 1119 : C 1120 : 1121 : IF(IKOMP.EQ.0)THEN 1122 : 1123 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008 1124 : IF(IDCENN.EQ.1)THEN 1125 : DO 80081 K=KDEB,KFIN 1126 : KP=K-KDEB+1 1127 : XMB(KP)=XMH(KP) 1128 : CXT(KP)=0.D0 1129 : CYT(KP)=0.D0 1130 : CZT(KP)=0.D0 1131 : CXY(KP)=0.D0 1132 : CXZ(KP)=0.D0 1133 : CYZ(KP)=0.D0 1134 : 80081 CONTINUE 1135 : 1136 : ELSEIF(IDCENN.EQ.2)THEN 1137 : DO 8008 K=KDEB,KFIN 1138 : KP=K-KDEB+1 1139 : HMK=2.D0*UM(KP)/BM(KP) 1140 : XMB(KP)=HMK 1141 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 1142 : AKSI=ALFA/3.D0 1143 : IF(ALFA.GT.3.D0)AKSI=1.D0 1144 : CCT=AKSI/BM(KP) 1145 : 1146 : HPK=2.D0*UP(KP)/BP(KP) 1147 : ALFA=UP(KP)*HPK/ALFT(KP)/2.D0 1148 : AKSI=ALFA/3.D0 1149 : IF(ALFA.GT.3.D0)AKSI=1.D0 1150 : CCP=AKSI/BP(KP) 1151 : CPT=CCP-CCT 1152 : CC2=0.D0 1153 : IF(CPT.GE.0.D0)CC2=CPT 1154 : 1155 : CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2) 1156 : CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2) 1157 : CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPI(KP,3)*UPI(KP,3)*CC2) 1158 : CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2) 1159 : CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPI(KP,1)*UPI(KP,3)*CC2) 1160 : CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPI(KP,2)*UPI(KP,3)*CC2) 1161 : 1162 : 8008 CONTINUE 1163 : 1164 : ELSEIF(IDCENN.EQ.3)THEN 1165 : DO 8018 K=KDEB,KFIN 1166 : KP=K-KDEB+1 1167 : HMK=2.D0*UM(KP)/BM(KP) 1168 : XMB(KP)=HMK 1169 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 1170 : AKSI=ALFA/3.D0 1171 : IF(ALFA.GT.3.D0)AKSI=1.D0 1172 : CCT=AKSI/BM(KP) 1173 : 1174 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 1175 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 1176 : CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 1177 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 1178 : CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 1179 : CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 1180 : 1181 : 8018 CONTINUE 1182 : 1183 : ELSEIF(IDCENN.EQ.4)THEN 1184 : DT19=DTM1*0.5D0 1185 : DO 8009 K=KDEB,KFIN 1186 : KP=K-KDEB+1 1187 : XMB(KP)=XMH(KP) 1188 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 1189 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 1190 : CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 1191 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 1192 : CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 1193 : CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 1194 : 8009 CONTINUE 1195 : 1196 : ENDIF 1197 : 1198 : ELSEIF(IKOMP.EQ.1)THEN 1199 : 1200 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008 1201 : IF(IDCENN.EQ.1)THEN 1202 : DO 81081 K=KDEB,KFIN 1203 : KP=K-KDEB+1 1204 : XMB(KP)=XMH(KP) 1205 : CXT(KP)=0.D0 1206 : CYT(KP)=0.D0 1207 : CZT(KP)=0.D0 1208 : CXY(KP)=0.D0 1209 : CXZ(KP)=0.D0 1210 : CYZ(KP)=0.D0 1211 : 1212 : DXT(KP)=0.D0 1213 : DYT(KP)=0.D0 1214 : DZT(KP)=0.D0 1215 : DXY(KP)=0.D0 1216 : DXZ(KP)=0.D0 1217 : DYZ(KP)=0.D0 1218 : 1219 : 81081 CONTINUE 1220 : 1221 : ELSEIF(IDCENN.EQ.2)THEN 1222 : DO 81008 K=KDEB,KFIN 1223 : KP=K-KDEB+1 1224 : HMK=2.D0*UM(KP)/BM(KP) 1225 : XMB(KP)=HMK 1226 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 1227 : AKSI=ALFA/3.D0 1228 : IF(ALFA.GT.3.D0)AKSI=1.D0 1229 : CCT=AKSI/BM(KP) 1230 : 1231 : HPK=2.D0*UP(KP)/BP(KP) 1232 : ALFA=UP(KP)*HPK/ALFT(KP)/2.D0 1233 : AKSI=ALFA/3.D0 1234 : IF(ALFA.GT.3.D0)AKSI=1.D0 1235 : CCP=AKSI/BP(KP) 1236 : CPT=CCP-CCT 1237 : CC2=0.D0 1238 : IF(CPT.GE.0.D0)CC2=CPT 1239 : CC2 = CC2*1.3D0 1240 : 1241 : CXT(KP)=UMI(KP,1)*CCT 1242 : CYT(KP)=UMI(KP,2)*CCT 1243 : CZT(KP)=UMI(KP,3)*CCT 1244 : CXY(KP)=CCT 1245 : CXZ(KP)=CCT 1246 : CYZ(KP)=CCT 1247 : 1248 : DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2 1249 : DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2 1250 : DZT(KP)=UPI(KP,3)*UPI(KP,3)*CC2 1251 : DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2 1252 : DXZ(KP)=UPI(KP,1)*UPI(KP,3)*CC2 1253 : DYZ(KP)=UPI(KP,2)*UPI(KP,3)*CC2 1254 : 1255 : 81008 CONTINUE 1256 : 1257 : ELSEIF(IDCENN.EQ.3)THEN 1258 : DO 81018 K=KDEB,KFIN 1259 : KP=K-KDEB+1 1260 : HMK=2.D0*UM(KP)/BM(KP) 1261 : XMB(KP)=HMK 1262 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 1263 : AKSI=ALFA/3.D0 1264 : IF(ALFA.GT.3.D0)AKSI=1.D0 1265 : CCT=AKSI/BM(KP) 1266 : 1267 : CXT(KP)=UMI(KP,1)*CCT 1268 : CYT(KP)=UMI(KP,2)*CCT 1269 : CZT(KP)=UMI(KP,3)*CCT 1270 : CXY(KP)=CCT 1271 : CXZ(KP)=CCT 1272 : CYZ(KP)=CCT 1273 : 1274 : DXT(KP)=0.D0 1275 : DYT(KP)=0.D0 1276 : DZT(KP)=0.D0 1277 : DXY(KP)=0.D0 1278 : DXZ(KP)=0.D0 1279 : DYZ(KP)=0.D0 1280 : 1281 : 81018 CONTINUE 1282 : 1283 : ELSEIF(IDCENN.EQ.4)THEN 1284 : DT19=DTM1*0.5D0 1285 : DO 81009 K=KDEB,KFIN 1286 : KP=K-KDEB+1 1287 : XMB(KP)=XMH(KP) 1288 : 1289 : CXT(KP)=UMI(KP,1)*DT19 1290 : CYT(KP)=UMI(KP,2)*DT19 1291 : CZT(KP)=UMI(KP,3)*DT19 1292 : CXY(KP)=DT19 1293 : CXZ(KP)=DT19 1294 : CYZ(KP)=DT19 1295 : 1296 : DXT(KP)=0.D0 1297 : DYT(KP)=0.D0 1298 : DZT(KP)=0.D0 1299 : DXY(KP)=0.D0 1300 : DXZ(KP)=0.D0 1301 : DYZ(KP)=0.D0 1302 : 1303 : 81009 CONTINUE 1304 : 1305 : ENDIF 1306 : 1307 : ENDIF 1308 : 1309 : C *********************** 1310 : C 1311 : C AVANT CALCUL DT 1312 : C 1313 : 1314 : IF(IKOMP.EQ.0)THEN 1315 : 1316 : C*IBMDIR* PREFER SCALAR 1317 : DO 8010 K=KDEB,KFIN 1318 : KP=K-KDEB+1 1319 : DT0=DT 1320 : 1321 : DT1=2.D0/ 1322 : & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)) 1323 : & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)) 1324 : & + (UMI(KP,3)*UMI(KP,3))/(ALFT(KP)+CZT(KP)) ) 1325 : 1326 : DT2=0.5/ 1327 : & ( (ALFT(KP)+CXT(KP))*AL2(KP) 1328 : & +(ALFT(KP)+CYT(KP))*AH2(KP) 1329 : & +(ALFT(KP)+CZT(KP))*AP2(KP) ) 1330 : 1331 : IF(DT1.LT.DT)DT=DT1 1332 : IF(DT2.LT.DT)DT=DT2 1333 : C IF(DT3.LT.DT)DT=DT3 1334 : IF(DT.NE.DT0) THEN 1335 : DTT1=DT1 1336 : DTT2=DT2 1337 : DTT3=0.D0 1338 : DIAEL=XMB(KP) 1339 : NUEL=K 1340 : END IF 1341 : 8010 CONTINUE 1342 : 1343 : ELSEIF(IKOMP.EQ.1)THEN 1344 : 1345 : C*IBMDIR* PREFER SCALAR 1346 : DO 8110 K=KDEB,KFIN 1347 : KP=K-KDEB+1 1348 : DT0=DT 1349 : 1350 : DT1=2.D0/ 1351 : & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP)) 1352 : & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP)) 1353 : & + (UMI(KP,3)*UMI(KP,3))/(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP)) ) 1354 : 1355 : DT2=0.5/ 1356 : & ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP) 1357 : & +(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP) 1358 : & +(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP))*AP2(KP) ) 1359 : 1360 : IF(DT1.LT.DT)DT=DT1 1361 : IF(DT2.LT.DT)DT=DT2 1362 : C IF(DT3.LT.DT)DT=DT3 1363 : IF(DT.NE.DT0) THEN 1364 : DTT1=DT1 1365 : DTT2=DT2 1366 : DTT3=0.D0 1367 : DIAEL=XMB(KP) 1368 : NUEL=K 1369 : ENDIF 1370 : 8110 CONTINUE 1371 : 1372 : ENDIF 1373 : 1374 : C Le coeur du calcul ... 1375 : 1376 : IF(IKOMP.EQ.0)THEN 1377 : 1378 : DO 8014 I=1,NP 1379 : DO 8014 J= 1,NP 1380 : C*IBMDIR* PREFER VECTOR 1381 : DO 8014 K=KDEB,KFIN 1382 : KP=K-KDEB+1 1383 : C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 1384 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 1385 : ZVGG=AIRE(KP)*( 1386 : & HR(K,I,1)*HR(K,J,1)*CXT(KP) 1387 : &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) 1388 : &+ HR(K,I,3)*HR(K,J,3)*CZT(KP) 1389 : &+ HR(K,I,1)*HR(K,J,2)*CXY(KP) 1390 : &+ HR(K,I,2)*HR(K,J,1)*CXY(KP) 1391 : &+ HR(K,I,1)*HR(K,J,3)*CXZ(KP) 1392 : &+ HR(K,I,3)*HR(K,J,1)*CXZ(KP) 1393 : &+ HR(K,I,2)*HR(K,J,3)*CYZ(KP) 1394 : &+ HR(K,I,3)*HR(K,J,2)*CYZ(KP) ) 1395 : &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1 1396 : C &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1397 : 1398 : ZVGT=AIRE(KP)*( 1399 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 1400 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 1401 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 1402 : &+ ALFT(KP)*XMH(KP)*GEO1 1403 : C &+ ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 1404 : 1405 : V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)*HR(K,J,2) 1406 : & +UMI(KP,3)*DR(KP,I)*HR(K,J,3) 1407 : 1408 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 1409 : 1410 : 8014 CONTINUE 1411 : 1412 : ELSEIF(IKOMP.EQ.1)THEN 1413 : 1414 : DO 8015 I=1,NP 1415 : DO 8015 J= 1,NP 1416 : C*IBMDIR* PREFER VECTOR 1417 : DO 8015 K=KDEB,KFIN 1418 : KP=K-KDEB+1 1419 : C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 1420 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 1421 : ZVGG=AIRE(KP)*( 1422 : & HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP)) 1423 : &+ HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP)) 1424 : &+ HR(K,I,3)*HR(K,J,3)*(CZT(KP)*UIZ(KP,J)+DZT(KP)) 1425 : &+ HR(K,I,1)*HR(K,J,2)*(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP)) 1426 : &+ HR(K,I,2)*HR(K,J,1)*(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP)) 1427 : &+ HR(K,I,1)*HR(K,J,3)*(CXZ(KP)*UMI(KP,1)*UIZ(KP,J)+DXZ(KP)) 1428 : &+ HR(K,I,3)*HR(K,J,1)*(CXZ(KP)*UMI(KP,3)*UIX(KP,J)+DXZ(KP)) 1429 : &+ HR(K,I,2)*HR(K,J,3)*(CYZ(KP)*UMI(KP,2)*UIZ(KP,J)+DYZ(KP)) 1430 : &+ HR(K,I,3)*HR(K,J,2)*(CYZ(KP)*UMI(KP,3)*UIY(KP,J)+DYZ(KP)) 1431 : &+ (CXT(KP)*UIX(KP,J)+DXT(KP)+CYT(KP)*UIY(KP,J)+DYT(KP) 1432 : &+ CZT(KP)*UIZ(KP,J)))/3.D0*GEO1 1433 : C &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1434 : 1435 : ZVGT=AIRE(KP)*( 1436 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 1437 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 1438 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 1439 : &+ ALFT(KP)*XMH(KP)*GEO1 1440 : C &+ ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 1441 : 1442 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+ 1443 : & UIZ(KP,J)*HR(K,J,3))*DR(KP,I) 1444 : 1445 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 1446 : 1447 : 8015 CONTINUE 1448 : 1449 : ENDIF 1450 : 1451 : C 1452 : C Fin de l'accumulation dans SBF. 1453 : C On ajoute ces incréments G. 1454 : C 1455 : DO 8017 I=1,NP 1456 : C*IBMDIR* PREFER VECTOR 1457 : DO 8017 K=KDEB,KFIN 1458 : KP=K-KDEB+1 1459 : NF=IPADL(LE(I,K)) 1460 : G(NF) = G(NF)+SBF(KP,I) 1461 : 8017 CONTINUE 1462 : 1463 : 8001 CONTINUE 1464 : 1465 : C WRITE(6,*)' G DANS YCTSCL ' 1466 : C WRITE(6,1984)(M,G(M),M=1,NPTD) 1467 : 1467 : C CALL ARRET(0) 1468 : IPAS=1 1469 : RETURN 1470 : 1002 FORMAT(10(1X,1PE11.4)) 1471 : END 1472 : 1473 : 1474 : 1475 : 1476 : 1477 : 1478 : 1479 :
© Cast3M 2003 - All rights reserved.
Disclaimer