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