1 : C ZCTSCA SOURCE CHAT 05/01/13 04:21:34 5004 2 : SUBROUTINE ZCTSCA 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 : C! Attention : IAXI et RPG ne sont pas utilisées... 12 : & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI, 13 : & IPADL,IKOMP,IKAS,IPADU,IPADD, 14 : & ALFE,IND1,UN,INDU,NPTU,TN,QE,IKS, 15 : & HRN,G,ALT,SGT, 16 : & VOLU,COTE,NELZ,IDCEN, 17 : & DTM1,DT,DTT1,DTT2,NUEL,DIAEL) 18 : 19 : IMPLICIT INTEGER(I-N) 20 : IMPLICIT REAL*8 (A-H,O-Z) 21 : 22 : C*********************************************************************** 23 : C 24 : C CE SP DISCRETISE UNE EQUATION GENERALE DE TRANSPORT-DIFFUSION AVEC 25 : C SOURCE. 26 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI 27 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 28 : C LES OPERATEURS SONT "SOUS-INTEGRES" 29 : C 30 : C 31 : C APPELE PAR YTSCAL 32 : C 33 : C 34 : C*********************************************************************** 35 : 36 : C-INC CCOPTIO 37 : -INC CCVQUA4 38 : -INC CCREEL 39 : C 40 : C Longueur des registres vectoriels de la machine cible 41 : C On prend 64 pour ne pas augmenter la taille des tableaux 42 : C nécessaires à la vectorisation. 43 : C 44 : PARAMETER(LRV=64) 45 : 46 : DIMENSION UN(NPTU,IES),HRN(NPTD),TN(NPTD) 47 : DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*) 48 : DIMENSION ALFE(*),ALT(*) 49 : 50 : DIMENSION IPADL(*),LE(NP,*),IPADU(*),IPADD(*) 51 : DIMENSION HR(NEL,NP,IES),RPG(*),DRR(NP,NEL) 52 : 53 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 54 : 55 : DIMENSION AIRE(LRV) 56 : DIMENSION AL2(LRV),AH2(LRV),AP2(LRV) 57 : DIMENSION AL(LRV),AH(LRV),AP(LRV) 58 : DIMENSION ALFT(LRV),QT(LRV) 59 : DIMENSION XMB(LRV),XMH(LRV) 60 : DIMENSION CFM(LRV) 61 : C DIMENSION CF1(LRV),CF2(LRV),CF3(LRV) 62 : DIMENSION DR(LRV,9) 63 : DIMENSION UM(LRV),UP(LRV) 64 : DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9) 65 : DIMENSION TETAC(LRV,9),TETAD(LRV,9) 66 : DIMENSION UMI(LRV,3),UPI(LRV,3) 67 : DIMENSION CXT(LRV),CYT(LRV),CXY(LRV) 68 : DIMENSION DXT(LRV),DYT(LRV),DXY(LRV) 69 : DIMENSION CZT(LRV),CXZ(LRV),CYZ(LRV) 70 : DIMENSION DZT(LRV),DXZ(LRV),DYZ(LRV) 71 : DIMENSION SBF(LRV,9) 72 : DIMENSION GRADT(LRV,3) 73 : DIMENSION BM(LRV),BP(LRV) 74 : 75 : REAL*8 G(*) 76 : 77 : SAVE IPAS,QGGT,Q1,Q2,Q3 78 : 79 : DATA CD/1.D0/ 80 : 81 : DATA IPAS/0/ 82 : C************************************************************************ 83 : C 84 : C INITIALISATIONS DIVERSES 85 : C 86 : C WRITE(IOIMP,*)' debut YCTSCL',IES 87 : ZERMA=1.D-20 88 : 89 : NK=K0 90 : C ******** 91 : C * 2D * 92 : C ******** 93 : 94 : IF(IES.EQ.3)GO TO 10 95 : 96 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 97 : 98 : C DIFFERENCES TRIANGLE / QUADRANGLE 99 : QUA4=0.D0 100 : IF(NP.EQ.4)QUA4=1.D0 101 : C 102 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 103 : 104 : C 105 : C Calcul du nombre de paquets de LRV éléments 106 : C 107 : NNN=MOD(NEL,LRV) 108 : IF(NNN.EQ.0) NPACK=NEL/LRV 109 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 110 : KPACKD=1 111 : KPACKF=NPACK 112 : C 113 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 114 : C 115 : DO 7001 KPACK=KPACKD,KPACKF 116 : C 117 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 118 : C 119 : C 1.D0 Calcul des limites du paquet courant. 120 : KDEB=1+(KPACK-1)*LRV 121 : KFIN=MIN(NEL,KDEB+LRV-1) 122 : C 123 : DO 7002 K=KDEB,KFIN 124 : KP=K-KDEB+1 125 : NK=K+K0 126 : NK1=(1-IND1)*(NK-1)+1 127 : ALFT(KP)=ALFE(NK1)+ZERMA 128 : AIRE(KP)=VOLU(NK) 129 : AL(KP)=COTE(NK,1) 130 : AH(KP)=COTE(NK,2) 131 : AL2(KP)=1.D0/AL(KP)/AL(KP) 132 : AH2(KP)=1.D0/AH(KP)/AH(KP) 133 : XMH(KP)=(AL(KP)+AH(KP))/2.D0 134 : 7002 CONTINUE 135 : 136 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 137 : & (IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 138 : DO 7003 K=KDEB,KFIN 139 : KP=K-KDEB+1 140 : NK=K+K0 141 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 142 : 7003 CONTINUE 143 : ENDIF 144 : 145 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 146 : 147 : C Initialisation des UMI avant accumulation 148 : DO 7005 K=KDEB,KFIN 149 : KP=K-KDEB+1 150 : UMI(KP,1)=ZERMA 151 : UMI(KP,2)=ZERMA 152 : UPI(KP,1)=ZERMA 153 : UPI(KP,2)=ZERMA 154 : GRADT(KP,1)=ZERMA 155 : GRADT(KP,2)=ZERMA 156 : 7005 CONTINUE 157 : 158 : DO 7006 I=1,NP 159 : C*IBMDIR* PREFER VECTOR 160 : DO 70061 K=KDEB,KFIN 161 : KP=K-KDEB+1 162 : NF=IPADL(LE(I,K)) 163 : NFD=IPADD(LE(I,K)) 164 : NFU=IPADU(LE(I,K)) 165 : NFU=(1-INDU)*(NFU-1)+1 166 : DR(KP,I)=DRR(I,K) 167 : UIX(KP,I)=UN(NFU,1) 168 : UIY(KP,I)=UN(NFU,2) 169 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 170 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 171 : 172 : TETAC(KP,I)=HRN(NF) 173 : TETAD(KP,I)=TN(NFD) 174 : GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I) 175 : GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I) 176 : 70061 CONTINUE 177 : 7006 CONTINUE 178 : 179 : C 180 : C Initialisation de la variable d'accumulation SBF au terme source 181 : C 182 : C WRITE(IOIMP,*)' IKomp,ikas=',IKomp,ikas 183 : C WRITE(IOIMP,*)' IKS,IND1,INDU=',IKS,IND1,INDU 184 : 185 : IF(IKOMP.EQ.0)THEN 186 : 187 : DO 70021 K=KDEB,KFIN 188 : KP=K-KDEB+1 189 : NK=K+K0 190 : NKS=(1-IKS)*(NK-1)+1 191 : QT(KP)=QE(NKS) 192 : 70021 CONTINUE 193 : 194 : DO 70062 I=1,NP 195 : C*IBMDIR* PREFER VECTOR 196 : DO 70162 K=KDEB,KFIN 197 : KP=K-KDEB+1 198 : SBF(KP,I)=-QT(KP)*DR(KP,I) 199 : 70162 CONTINUE 200 : 70062 CONTINUE 201 : 202 : ELSEIF(IKOMP.EQ.1)THEN 203 : 204 : DO 70023 K=KDEB,KFIN 205 : KP=K-KDEB+1 206 : NK=K+K0 207 : NKS=(1-IKS)*(NK-1)+1 208 : QT(KP)=QE(NKS) 209 : 70023 CONTINUE 210 : 211 : DO 70064 I=1,NP 212 : C*IBMDIR* PREFER VECTOR 213 : DO 70164 K=KDEB,KFIN 214 : KP=K-KDEB+1 215 : SBF(KP,I)=-QT(KP)*DR(KP,I) 216 : 70164 CONTINUE 217 : 70064 CONTINUE 218 : 219 : ENDIF 220 : 221 : DO 7007 K=KDEB,KFIN 222 : KP=K-KDEB+1 223 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 224 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 225 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 226 : UM(KP)=SQRT(UM(KP)) 227 : A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2) 228 : A=SQRT(A)+XPETIT 229 : GRADT(KP,1)=GRADT(KP,1)/A 230 : GRADT(KP,2)=GRADT(KP,2)/A 231 : UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2) 232 : UPI(KP,1)=UP(KP)*GRADT(KP,1) 233 : UPI(KP,2)=UP(KP)*GRADT(KP,2) 234 : 7007 CONTINUE 235 : 236 : DO 70074 K=KDEB,KFIN 237 : KP=K-KDEB+1 238 : BMX=UMI(KP,1)/AL(KP) 239 : BMY=UMI(KP,2)/AH(KP) 240 : BM(KP)=BMX*BMX+BMY*BMY 241 : BM(KP)=SQRT(BM(KP))+XPETIT 242 : BPX=UPI(KP,1)/AL(KP) 243 : BPY=UPI(KP,2)/AH(KP) 244 : BP(KP)=BPX*BPX+BPY*BPY 245 : BP(KP)=SQRT(BP(KP))+XPETIT 246 : 70074 CONTINUE 247 : 248 : C 249 : C AVANT DECENTREMENT 250 : C 251 : 252 : IF(IKOMP.EQ.0)THEN 253 : 254 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008 255 : IF(IDCEN.EQ.1)THEN 256 : DO 70081 K=KDEB,KFIN 257 : KP=K-KDEB+1 258 : XMB(KP)=XMH(KP) 259 : CXT(KP)=0.D0 260 : CYT(KP)=0.D0 261 : CXY(KP)=0.D0 262 : 70081 CONTINUE 263 : 264 : ELSEIF(IDCEN.EQ.2)THEN 265 : DO 7008 K=KDEB,KFIN 266 : KP=K-KDEB+1 267 : HMK=2.D0*UM(KP)/BM(KP) 268 : XMB(KP)=HMK 269 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 270 : AKSI=ALFA/3.D0 271 : IF(ALFA.GT.3.D0)AKSI=1.D0 272 : CCT=AKSI/BM(KP) 273 : 274 : HPK=2.D0*UP(KP)/BP(KP) 275 : ALFA=UP(KP)*HPK/ALFT(KP)/2.D0 276 : AKSI=ALFA/3.D0 277 : IF(ALFA.GT.3.D0)AKSI=1.D0 278 : CCP=AKSI/BP(KP) 279 : CPT=CCP-CCT 280 : CC2=0.D0 281 : IF(CPT.GE.0.D0)CC2=CPT 282 : 283 : CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT 284 : $ +UPI(KP,1)*UPI(KP,1)*CC2) 285 : CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT 286 : $ +UPI(KP,2)*UPI(KP,2)*CC2) 287 : CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT 288 : $ +UPI(KP,1)*UPI(KP,2)*CC2) 289 : 290 : 7008 CONTINUE 291 : 292 : ELSEIF(IDCEN.EQ.3)THEN 293 : DO 7018 K=KDEB,KFIN 294 : KP=K-KDEB+1 295 : HMK=2.D0*UM(KP)/BM(KP) 296 : XMB(KP)=HMK 297 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 298 : AKSI=ALFA/3.D0 299 : IF(ALFA.GT.3.D0)AKSI=1.D0 300 : CCT=AKSI/BM(KP) 301 : 302 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 303 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 304 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 305 : 306 : 7018 CONTINUE 307 : 308 : ELSEIF(IDCEN.EQ.4)THEN 309 : DT19=DTM1*0.5D0 310 : DO 7009 K=KDEB,KFIN 311 : KP=K-KDEB+1 312 : XMB(KP)=XMH(KP) 313 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 314 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 315 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 316 : 7009 CONTINUE 317 : 318 : ENDIF 319 : 320 : ELSEIF(IKOMP.EQ.1)THEN 321 : 322 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008 323 : IF(IDCEN.EQ.1)THEN 324 : DO 71081 K=KDEB,KFIN 325 : KP=K-KDEB+1 326 : XMB(KP)=XMH(KP) 327 : 328 : CXT(KP)=0.D0 329 : CYT(KP)=0.D0 330 : CXY(KP)=0.D0 331 : 332 : DXT(KP)=0.D0 333 : DYT(KP)=0.D0 334 : DXY(KP)=0.D0 335 : 336 : 71081 CONTINUE 337 : 338 : ELSEIF(IDCEN.EQ.2)THEN 339 : DO 7108 K=KDEB,KFIN 340 : KP=K-KDEB+1 341 : HMK=2.D0*UM(KP)/BM(KP) 342 : XMB(KP)=HMK 343 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 344 : AKSI=ALFA/3.D0 345 : IF(ALFA.GT.3.D0)AKSI=1.D0 346 : CCT=AKSI/BM(KP) 347 : 348 : HPK=2.D0*UP(KP)/BP(KP) 349 : ALFA=UP(KP)*HPK/ALFT(KP)/2.D0 350 : AKSI=ALFA/3.D0 351 : IF(ALFA.GT.3.D0)AKSI=1.D0 352 : CCP=AKSI/BP(KP) 353 : CPT=CCP-CCT 354 : CC2=0.D0 355 : IF(CPT.GE.0.D0)CC2=CPT 356 : CC2 = CC2*1.3D0 357 : 358 : CXT(KP)=UMI(KP,1)*CCT 359 : CYT(KP)=UMI(KP,2)*CCT 360 : CXY(KP)=CCT 361 : 362 : DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2 363 : DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2 364 : DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2 365 : 366 : 7108 CONTINUE 367 : 368 : ELSEIF(IDCEN.EQ.3)THEN 369 : DO 71018 K=KDEB,KFIN 370 : KP=K-KDEB+1 371 : HMK=2.D0*UM(KP)/BM(KP) 372 : XMB(KP)=HMK 373 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 374 : AKSI=ALFA/3.D0 375 : IF(ALFA.GT.3.D0)AKSI=1.D0 376 : CCT=AKSI/BM(KP) 377 : 378 : CXT(KP)=UMI(KP,1)*CCT 379 : CYT(KP)=UMI(KP,2)*CCT 380 : CXY(KP)=CCT 381 : 382 : DXT(KP)=0.D0 383 : DYT(KP)=0.D0 384 : DXY(KP)=0.D0 385 : 386 : 71018 CONTINUE 387 : 388 : ELSEIF(IDCEN.EQ.4)THEN 389 : DT19=DTM1*0.5D0 390 : DO 71009 K=KDEB,KFIN 391 : KP=K-KDEB+1 392 : XMB(KP)=XMH(KP) 393 : CXT(KP)=UMI(KP,1)*DT19 394 : CYT(KP)=UMI(KP,2)*DT19 395 : CXY(KP)=DT19 396 : 397 : DXT(KP)=0.D0 398 : DYT(KP)=0.D0 399 : DXY(KP)=0.D0 400 : 401 : 71009 CONTINUE 402 : 403 : ENDIF 404 : 405 : ENDIF 406 : C *********************** 407 : C 408 : C AVANT CALCUL DT 409 : C 410 : IF(IKOMP.EQ.0)THEN 411 : 412 : C*IBMDIR* PREFER SCALAR 413 : DO 7010 K=KDEB,KFIN 414 : KP=K-KDEB+1 415 : DT0=DT 416 : 417 : DT1=2.D0/ 418 : & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)) 419 : & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)) ) 420 : 421 : DT2=0.5D0/ 422 : & ( (ALFT(KP)+CXT(KP))*AL2(KP) 423 : & +(ALFT(KP)+CYT(KP))*AH2(KP) ) 424 : 425 : IF(DT1.LT.DT)DT=DT1 426 : IF(DT2.LT.DT)DT=DT2 427 : C IF(DT3.LT.DT)DT=DT3 428 : IF(DT.NE.DT0) THEN 429 : DTT1=DT1 430 : DTT2=DT2 431 : C DTT3=0.D0 432 : DIAEL=XMB(KP) 433 : NUEL=K 434 : END IF 435 : CXY(KP)=CXY(KP)*CD 436 : 7010 CONTINUE 437 : 438 : ELSEIF(IKOMP.EQ.1)THEN 439 : 440 : C*IBMDIR* PREFER SCALAR 441 : DO 7110 K=KDEB,KFIN 442 : KP=K-KDEB+1 443 : DT0=DT 444 : 445 : DT1=2.D0/ 446 : & ( (UMI(KP,1)*UMI(KP,1)) 447 : $ /(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP)) 448 : $ + (UMI(KP,2)*UMI(KP,2)) 449 : $ /(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP)) ) 450 : 451 : DT2=0.5D0/ 452 : & ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP) 453 : & + (ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP) ) 454 : 455 : IF(DT1.LT.DT)DT=DT1 456 : IF(DT2.LT.DT)DT=DT2 457 : C IF(DT3.LT.DT)DT=DT3 458 : IF(DT.NE.DT0) THEN 459 : DTT1=DT1 460 : DTT2=DT2 461 : C DTT3=0.D0 462 : DIAEL=XMB(KP) 463 : NUEL=K 464 : END IF 465 : CXY(KP)=CXY(KP)*CD 466 : 7110 CONTINUE 467 : 468 : ENDIF 469 : C Le coeur du calcul ... 470 : 471 : IF(IKOMP.EQ.0)THEN 472 : 473 : DO 7014 I=1,NP 474 : DO 70141 J= 1,NP 475 : C*IBMDIR* PREFER VECTOR 476 : DO 70142 K=KDEB,KFIN 477 : KP=K-KDEB+1 478 : ZVGG=AIRE(KP)*( 479 : & HR(K,I,1)*HR(K,J,1)*CXT(KP) 480 : & + HR(K,I,1)*HR(K,J,2)*CXY(KP) 481 : & + HR(K,I,2)*HR(K,J,1)*CXY(KP) 482 : & + HR(K,I,2)*HR(K,J,2)*CYT(KP) 483 : & + (CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0 484 : & *VGGT(J,I)*QUA4 ) 485 : 486 : C &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) ) 487 : C &+ (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4 488 : 489 : ZVGT=AIRE(KP)*( 490 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 491 : & + HR(K,I,2)*HR(K,J,2)*ALFT(KP) 492 : & + (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0 493 : & *VGGT(J,I)*QUA4 ) 494 : 495 : C &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) ) 496 : C &+ ALFT(KP)/12.D0*VGGT(J,I)*QUA4 497 : 498 : C? V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I) 499 : V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1) 500 : $ +UMI(KP,2)*DR(KP,I)*HR(K,J,2) 501 : 502 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2) 503 : $ +TETAD(KP,J)*ZVGT 504 : 505 : 70142 CONTINUE 506 : 70141 CONTINUE 507 : 7014 CONTINUE 508 : 509 : ELSEIF(IKOMP.EQ.1)THEN 510 : 511 : DO 7015 I=1,NP 512 : DO 70151 J= 1,NP 513 : C*IBMDIR* PREFER VECTOR 514 : DO 70152 K=KDEB,KFIN 515 : KP=K-KDEB+1 516 : ZVGG=AIRE(KP)*( 517 : & HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP)) 518 : &+ HR(K,I,1)*HR(K,J,2) 519 : & *(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP)) 520 : &+ HR(K,I,2)*HR(K,J,1) 521 : & *(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP)) 522 : &+ HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP)) 523 : &+ ((CXT(KP)*UIX(KP,J)+DXT(KP))*AL2(KP) 524 : & +(CYT(KP)*UIY(KP,J)+DYT(KP))*AH2(KP))/12.D0 525 : & *VGGT(J,I)*QUA4 ) 526 : 527 : C &+ HR(K,I,2)*HR(K,J,2)*CYT(KP) ) 528 : C &+ (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4 529 : 530 : ZVGT=AIRE(KP)*( 531 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 532 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 533 : &+ (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0 534 : & *VGGT(J,I)*QUA4 ) 535 : 536 : C &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) ) 537 : C &+ ALFT(KP)/12.D0*VGGT(J,I)*QUA4 538 : 539 : V2=(UIX(KP,J)*HR(K,J,1) 540 : $ +UIY(KP,J)*HR(K,J,2))*DR(KP,I) 541 : 542 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2) 543 : $ +TETAD(KP,J)*ZVGT 544 : 545 : 70152 CONTINUE 546 : 70151 CONTINUE 547 : 7015 CONTINUE 548 : 549 : ENDIF 550 : C 551 : C Fin de l'accumulation dans SBF. 552 : C On ajoute ces incréments G. 553 : C 554 : DO 7017 I=1,NP 555 : C*IBMDIR* PREFER VECTOR 556 : DO 70171 K=KDEB,KFIN 557 : KP=K-KDEB+1 558 : NF=IPADL(LE(I,K)) 559 : G(NF) = G(NF)+SBF(KP,I) 560 : 70171 CONTINUE 561 : 7017 CONTINUE 562 : 563 : 7001 CONTINUE 564 : 565 : C WRITE(IOIMP,*)' G DANS YCTSCL ' 566 : C WRITE(IOIMP,1984)(M,G(M),M=1,NPTD) 567 : 1984 FORMAT(7(1X,I4,2X,1PE11.4)) 568 : 568 : C CALL ARRET(0) 569 : IPAS=1 570 : RETURN 571 : 572 : C ******** 573 : C * 3D * 574 : C ******** 575 : 576 : 10 CONTINUE 577 : 577 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 578 : CUB8=0.D0 579 : IF(NP.EQ.8)CUB8=1.D0 580 : 581 : C 582 : C Calcul du nombre de paquets de LRV éléments 583 : C 584 : NNN=MOD(NEL,LRV) 585 : IF(NNN.EQ.0) NPACK=NEL/LRV 586 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 587 : KPACKD=1 588 : KPACKF=NPACK 589 : C 590 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 591 : C 592 : DO 8001 KPACK=KPACKD,KPACKF 593 : C 594 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 595 : C 596 : C 1.D0 Calcul des limites du paquet courant. 597 : KDEB=1+(KPACK-1)*LRV 598 : KFIN=MIN(NEL,KDEB+LRV-1) 599 : C 600 : DO 8002 K=KDEB,KFIN 601 : KP=K-KDEB+1 602 : NK=K+K0 603 : NK1=(1-IND1)*(NK-1)+1 604 : ALFT(KP)=ALFE(NK1)+ZERMA 605 : AIRE(KP)=VOLU(NK) 606 : AL(KP)=COTE(NK,1) 607 : AH(KP)=COTE(NK,2) 608 : AP(KP)=COTE(NK,3) 609 : CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ 610 : & AP(KP)*AH(KP)/AL(KP) 611 : C CF1(KP)=AL(KP)*AH(KP)/AP(KP) 612 : C CF2(KP)=AL(KP)*AP(KP)/AH(KP) 613 : C CF3(KP)=AP(KP)*AH(KP)/AL(KP) 614 : XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0 615 : AL2(KP)=1.D0/AL(KP)/AL(KP) 616 : AH2(KP)=1.D0/AH(KP)/AH(KP) 617 : AP2(KP)=1.D0/AP(KP)/AP(KP) 618 : 8002 CONTINUE 619 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 620 : & (IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 621 : DO 8003 K=KDEB,KFIN 622 : KP=K-KDEB+1 623 : NK=K+K0 624 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 625 : 8003 CONTINUE 626 : ENDIF 627 : 628 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 629 : 630 : C Initialisation des UMI avant accumulation 631 : DO 8005 K=KDEB,KFIN 632 : KP=K-KDEB+1 633 : UMI(KP,1)=ZERMA 634 : UMI(KP,2)=ZERMA 635 : UMI(KP,3)=ZERMA 636 : UPI(KP,1)=ZERMA 637 : UPI(KP,2)=ZERMA 638 : UPI(KP,3)=ZERMA 639 : GRADT(KP,1)=ZERMA 640 : GRADT(KP,2)=ZERMA 641 : GRADT(KP,3)=ZERMA 642 : 8005 CONTINUE 643 : 644 : DO 8006 I=1,NP 645 : C*IBMDIR* PREFER VECTOR 646 : DO 80061 K=KDEB,KFIN 647 : KP=K-KDEB+1 648 : NF=IPADL(LE(I,K)) 649 : NFD=IPADD(LE(I,K)) 650 : NFU=IPADU(LE(I,K)) 651 : NFU=(1-INDU)*(NFU-1)+1 652 : DR(KP,I)=DRR(I,K) 653 : UIX(KP,I)=UN(NFU,1) 654 : UIY(KP,I)=UN(NFU,2) 655 : UIZ(KP,I)=UN(NFU,3) 656 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 657 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 658 : UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I) 659 : 660 : TETAC(KP,I)=HRN(NF) 661 : TETAD(KP,I)=TN(NFD) 662 : GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I) 663 : GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I) 664 : GRADT(KP,3)=GRADT(KP,3)+HR(K,I,3)*TETAC(KP,I) 665 : 80061 CONTINUE 666 : 8006 CONTINUE 667 : 668 : C 669 : C Initialisation de la variable d'accumulation SBF au terme source 670 : C M 671 : C WRITE(IOIMP,*)' IKomp,ikas=',IKomp,ikas 672 : C WRITE(IOIMP,*)' IKS,IND1,INDU=',IKS,IND1,INDU 673 : 674 : IF(IKOMP.EQ.0)THEN 675 : 676 : DO 80021 K=KDEB,KFIN 677 : KP=K-KDEB+1 678 : NK=K+K0 679 : NKS=(1-IKS)*(NK-1)+1 680 : QT(KP)=QE(NKS) 681 : 80021 CONTINUE 682 : 683 : DO 80062 I=1,NP 684 : C*IBMDIR* PREFER VECTOR 685 : DO 80162 K=KDEB,KFIN 686 : KP=K-KDEB+1 687 : SBF(KP,I)=-QT(KP)*DR(KP,I) 688 : 80162 CONTINUE 689 : 80062 CONTINUE 690 : 691 : ELSEIF(IKOMP.EQ.1)THEN 692 : 693 : DO 80023 K=KDEB,KFIN 694 : KP=K-KDEB+1 695 : NK=K+K0 696 : NKS=(1-IKS)*(NK-1)+1 697 : QT(KP)=QE(NKS) 698 : 80023 CONTINUE 699 : 700 : DO 80064 I=1,NP 701 : C*IBMDIR* PREFER VECTOR 702 : DO 80164 K=KDEB,KFIN 703 : KP=K-KDEB+1 704 : SBF(KP,I)=-QT(KP)*DR(KP,I) 705 : 80164 CONTINUE 706 : 80064 CONTINUE 707 : 708 : ENDIF 709 : 710 : DO 8007 K=KDEB,KFIN 711 : KP=K-KDEB+1 712 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 713 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 714 : UMI(KP,3)=UMI(KP,3)/AIRE(KP) 715 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 716 : & +UMI(KP,3)*UMI(KP,3) 717 : UM(KP)=SQRT(UM(KP)) 718 : A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2) 719 : & +GRADT(KP,3)*GRADT(KP,3) 720 : A=SQRT(A)+XPETIT 721 : GRADT(KP,1)=GRADT(KP,1)/A 722 : GRADT(KP,2)=GRADT(KP,2)/A 723 : GRADT(KP,3)=GRADT(KP,3)/A 724 : UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2) 725 : & +GRADT(KP,3)*UMI(KP,3) 726 : UPI(KP,1)=UP(KP)*GRADT(KP,1) 727 : UPI(KP,2)=UP(KP)*GRADT(KP,2) 728 : UPI(KP,3)=UP(KP)*GRADT(KP,3) 729 : 8007 CONTINUE 730 : 731 : DO 80074 K=KDEB,KFIN 732 : KP=K-KDEB+1 733 : BMX=UMI(KP,1)/AL(KP) 734 : BMY=UMI(KP,2)/AH(KP) 735 : BMZ=UMI(KP,3)/AP(KP) 736 : BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ 737 : BM(KP)=SQRT(BM(KP))+XPETIT 738 : BPX=UPI(KP,1)/AL(KP) 739 : BPY=UPI(KP,2)/AH(KP) 740 : BPZ=UPI(KP,3)/AP(KP) 741 : BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ 742 : BP(KP)=SQRT(BP(KP))+XPETIT 743 : 80074 CONTINUE 744 : 745 : C 746 : C AVANT DECENTREMENT 747 : C 748 : 749 : IF(IKOMP.EQ.0)THEN 750 : 751 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008 752 : IF(IDCEN.EQ.1)THEN 753 : DO 80081 K=KDEB,KFIN 754 : KP=K-KDEB+1 755 : XMB(KP)=XMH(KP) 756 : CXT(KP)=0.D0 757 : CYT(KP)=0.D0 758 : CZT(KP)=0.D0 759 : CXY(KP)=0.D0 760 : CXZ(KP)=0.D0 761 : CYZ(KP)=0.D0 762 : 80081 CONTINUE 763 : 764 : ELSEIF(IDCEN.EQ.2)THEN 765 : DO 8008 K=KDEB,KFIN 766 : KP=K-KDEB+1 767 : HMK=2.D0*UM(KP)/BM(KP) 768 : XMB(KP)=HMK 769 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 770 : AKSI=ALFA/3.D0 771 : IF(ALFA.GT.3.D0)AKSI=1.D0 772 : CCT=AKSI/BM(KP) 773 : 774 : HPK=2.D0*UP(KP)/BP(KP) 775 : ALFA=UP(KP)*HPK/ALFT(KP)/2.D0 776 : AKSI=ALFA/3.D0 777 : IF(ALFA.GT.3.D0)AKSI=1.D0 778 : CCP=AKSI/BP(KP) 779 : CPT=CCP-CCT 780 : CC2=0.D0 781 : IF(CPT.GE.0.D0)CC2=CPT 782 : 783 : CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT 784 : $ +UPI(KP,1)*UPI(KP,1)*CC2) 785 : CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT 786 : $ +UPI(KP,2)*UPI(KP,2)*CC2) 787 : CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT 788 : $ +UPI(KP,3)*UPI(KP,3)*CC2) 789 : CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT 790 : $ +UPI(KP,1)*UPI(KP,2)*CC2) 791 : CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT 792 : $ +UPI(KP,1)*UPI(KP,3)*CC2) 793 : CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT 794 : $ +UPI(KP,2)*UPI(KP,3)*CC2) 795 : 796 : 8008 CONTINUE 797 : 798 : ELSEIF(IDCEN.EQ.3)THEN 799 : DO 8018 K=KDEB,KFIN 800 : KP=K-KDEB+1 801 : HMK=2.D0*UM(KP)/BM(KP) 802 : XMB(KP)=HMK 803 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 804 : AKSI=ALFA/3.D0 805 : IF(ALFA.GT.3.D0)AKSI=1.D0 806 : CCT=AKSI/BM(KP) 807 : 808 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 809 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 810 : CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 811 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 812 : CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 813 : CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 814 : 815 : 8018 CONTINUE 816 : 817 : ELSEIF(IDCEN.EQ.4)THEN 818 : DT19=DTM1*0.5D0 819 : DO 8009 K=KDEB,KFIN 820 : KP=K-KDEB+1 821 : XMB(KP)=XMH(KP) 822 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 823 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 824 : CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 825 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 826 : CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 827 : CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 828 : 8009 CONTINUE 829 : 830 : ENDIF 831 : 832 : ELSEIF(IKOMP.EQ.1)THEN 833 : 834 : C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008 835 : IF(IDCEN.EQ.1)THEN 836 : DO 81081 K=KDEB,KFIN 837 : KP=K-KDEB+1 838 : XMB(KP)=XMH(KP) 839 : CXT(KP)=0.D0 840 : CYT(KP)=0.D0 841 : CZT(KP)=0.D0 842 : CXY(KP)=0.D0 843 : CXZ(KP)=0.D0 844 : CYZ(KP)=0.D0 845 : 846 : DXT(KP)=0.D0 847 : DYT(KP)=0.D0 848 : DZT(KP)=0.D0 849 : DXY(KP)=0.D0 850 : DXZ(KP)=0.D0 851 : DYZ(KP)=0.D0 852 : 853 : 81081 CONTINUE 854 : 855 : ELSEIF(IDCEN.EQ.2)THEN 856 : DO 81008 K=KDEB,KFIN 857 : KP=K-KDEB+1 858 : HMK=2.D0*UM(KP)/BM(KP) 859 : XMB(KP)=HMK 860 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 861 : AKSI=ALFA/3.D0 862 : IF(ALFA.GT.3.D0)AKSI=1.D0 863 : CCT=AKSI/BM(KP) 864 : 865 : HPK=2.D0*UP(KP)/BP(KP) 866 : ALFA=UP(KP)*HPK/ALFT(KP)/2.D0 867 : AKSI=ALFA/3.D0 868 : IF(ALFA.GT.3.D0)AKSI=1.D0 869 : CCP=AKSI/BP(KP) 870 : CPT=CCP-CCT 871 : CC2=0.D0 872 : IF(CPT.GE.0.D0)CC2=CPT 873 : CC2 = CC2*1.3D0 874 : 875 : CXT(KP)=UMI(KP,1)*CCT 876 : CYT(KP)=UMI(KP,2)*CCT 877 : CZT(KP)=UMI(KP,3)*CCT 878 : CXY(KP)=CCT 879 : CXZ(KP)=CCT 880 : CYZ(KP)=CCT 881 : 882 : DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2 883 : DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2 884 : DZT(KP)=UPI(KP,3)*UPI(KP,3)*CC2 885 : DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2 886 : DXZ(KP)=UPI(KP,1)*UPI(KP,3)*CC2 887 : DYZ(KP)=UPI(KP,2)*UPI(KP,3)*CC2 888 : 889 : 81008 CONTINUE 890 : 891 : ELSEIF(IDCEN.EQ.3)THEN 892 : DO 81018 K=KDEB,KFIN 893 : KP=K-KDEB+1 894 : HMK=2.D0*UM(KP)/BM(KP) 895 : XMB(KP)=HMK 896 : ALFA=UM(KP)*HMK/ALFT(KP)/2.D0 897 : AKSI=ALFA/3.D0 898 : IF(ALFA.GT.3.D0)AKSI=1.D0 899 : CCT=AKSI/BM(KP) 900 : 901 : CXT(KP)=UMI(KP,1)*CCT 902 : CYT(KP)=UMI(KP,2)*CCT 903 : CZT(KP)=UMI(KP,3)*CCT 904 : CXY(KP)=CCT 905 : CXZ(KP)=CCT 906 : CYZ(KP)=CCT 907 : 908 : DXT(KP)=0.D0 909 : DYT(KP)=0.D0 910 : DZT(KP)=0.D0 911 : DXY(KP)=0.D0 912 : DXZ(KP)=0.D0 913 : DYZ(KP)=0.D0 914 : 915 : 81018 CONTINUE 916 : 917 : ELSEIF(IDCEN.EQ.4)THEN 918 : DT19=DTM1*0.5D0 919 : DO 81009 K=KDEB,KFIN 920 : KP=K-KDEB+1 921 : XMB(KP)=XMH(KP) 922 : 923 : CXT(KP)=UMI(KP,1)*DT19 924 : CYT(KP)=UMI(KP,2)*DT19 925 : CZT(KP)=UMI(KP,3)*DT19 926 : CXY(KP)=DT19 927 : CXZ(KP)=DT19 928 : CYZ(KP)=DT19 929 : 930 : DXT(KP)=0.D0 931 : DYT(KP)=0.D0 932 : DZT(KP)=0.D0 933 : DXY(KP)=0.D0 934 : DXZ(KP)=0.D0 935 : DYZ(KP)=0.D0 936 : 937 : 81009 CONTINUE 938 : 939 : ENDIF 940 : 941 : ENDIF 942 : 943 : C *********************** 944 : C 945 : C AVANT CALCUL DT 946 : C 947 : 948 : IF(IKOMP.EQ.0)THEN 949 : 950 : C*IBMDIR* PREFER SCALAR 951 : DO 8010 K=KDEB,KFIN 952 : KP=K-KDEB+1 953 : DT0=DT 954 : 955 : DT1=2.D0/ 956 : & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)) 957 : & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)) 958 : & + (UMI(KP,3)*UMI(KP,3))/(ALFT(KP)+CZT(KP)) ) 959 : 960 : DT2=0.5D0/ 961 : & ( (ALFT(KP)+CXT(KP))*AL2(KP) 962 : & +(ALFT(KP)+CYT(KP))*AH2(KP) 963 : & +(ALFT(KP)+CZT(KP))*AP2(KP) ) 964 : 965 : IF(DT1.LT.DT)DT=DT1 966 : IF(DT2.LT.DT)DT=DT2 967 : C IF(DT3.LT.DT)DT=DT3 968 : IF(DT.NE.DT0) THEN 969 : DTT1=DT1 970 : DTT2=DT2 971 : C DTT3=0.D0 972 : DIAEL=XMB(KP) 973 : NUEL=K 974 : END IF 975 : 8010 CONTINUE 976 : 977 : ELSEIF(IKOMP.EQ.1)THEN 978 : 979 : C*IBMDIR* PREFER SCALAR 980 : DO 8110 K=KDEB,KFIN 981 : KP=K-KDEB+1 982 : DT0=DT 983 : 984 : DT1=2.D0/ 985 : & ( (UMI(KP,1)*UMI(KP,1)) 986 : $ /(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP)) 987 : $ + (UMI(KP,2)*UMI(KP,2)) 988 : $ /(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP)) 989 : $ + (UMI(KP,3)*UMI(KP,3)) 990 : $ /(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP)) ) 991 : 992 : DT2=0.5D0/ 993 : & ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP) 994 : & +(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP) 995 : & +(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP))*AP2(KP) ) 996 : 997 : IF(DT1.LT.DT)DT=DT1 998 : IF(DT2.LT.DT)DT=DT2 999 : C IF(DT3.LT.DT)DT=DT3 1000 : IF(DT.NE.DT0) THEN 1001 : DTT1=DT1 1002 : DTT2=DT2 1003 : C DTT3=0.D0 1004 : DIAEL=XMB(KP) 1005 : NUEL=K 1006 : ENDIF 1007 : 8110 CONTINUE 1008 : 1009 : ENDIF 1010 : 1011 : C Le coeur du calcul ... 1012 : 1013 : IF(IKOMP.EQ.0)THEN 1014 : 1015 : DO 8014 I=1,NP 1016 : DO 80141 J= 1,NP 1017 : C*IBMDIR* PREFER VECTOR 1018 : DO 80142 K=KDEB,KFIN 1019 : KP=K-KDEB+1 1020 : C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 1021 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 1022 : ZVGG=AIRE(KP)*( 1023 : & HR(K,I,1)*HR(K,J,1)*CXT(KP) 1024 : & + HR(K,I,2)*HR(K,J,2)*CYT(KP) 1025 : & + HR(K,I,3)*HR(K,J,3)*CZT(KP) 1026 : & + HR(K,I,1)*HR(K,J,2)*CXY(KP) 1027 : & + HR(K,I,2)*HR(K,J,1)*CXY(KP) 1028 : & + HR(K,I,1)*HR(K,J,3)*CXZ(KP) 1029 : & + HR(K,I,3)*HR(K,J,1)*CXZ(KP) 1030 : & + HR(K,I,2)*HR(K,J,3)*CYZ(KP) 1031 : & + HR(K,I,3)*HR(K,J,2)*CYZ(KP) ) 1032 : & + (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1 1033 : C &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1034 : 1035 : ZVGT=AIRE(KP)*( 1036 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 1037 : & + HR(K,I,2)*HR(K,J,2)*ALFT(KP) 1038 : & + HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 1039 : & + ALFT(KP)*XMH(KP)*GEO1 1040 : C &+ ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 1041 : 1042 : V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1) 1043 : $ +UMI(KP,2)*DR(KP,I)*HR(K,J,2) 1044 : $ +UMI(KP,3)*DR(KP,I)*HR(K,J,3) 1045 : 1046 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2) 1047 : $ +TETAD(KP,J)*ZVGT 1048 : 1049 : 80142 CONTINUE 1050 : 80141 CONTINUE 1051 : 8014 CONTINUE 1052 : 1053 : ELSEIF(IKOMP.EQ.1)THEN 1054 : 1055 : DO 8015 I=1,NP 1056 : DO 80151 J= 1,NP 1057 : C*IBMDIR* PREFER VECTOR 1058 : DO 80152 K=KDEB,KFIN 1059 : KP=K-KDEB+1 1060 : C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 1061 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 1062 : ZVGG=AIRE(KP)*( 1063 : & HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP)) 1064 : &+ HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP)) 1065 : &+ HR(K,I,3)*HR(K,J,3)*(CZT(KP)*UIZ(KP,J)+DZT(KP)) 1066 : &+ HR(K,I,1)*HR(K,J,2)*(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP)) 1067 : &+ HR(K,I,2)*HR(K,J,1)*(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP)) 1068 : &+ HR(K,I,1)*HR(K,J,3)*(CXZ(KP)*UMI(KP,1)*UIZ(KP,J)+DXZ(KP)) 1069 : &+ HR(K,I,3)*HR(K,J,1)*(CXZ(KP)*UMI(KP,3)*UIX(KP,J)+DXZ(KP)) 1070 : &+ HR(K,I,2)*HR(K,J,3)*(CYZ(KP)*UMI(KP,2)*UIZ(KP,J)+DYZ(KP)) 1071 : &+ HR(K,I,3)*HR(K,J,2)*(CYZ(KP)*UMI(KP,3)*UIY(KP,J)+DYZ(KP)) 1072 : &+ (CXT(KP)*UIX(KP,J)+DXT(KP)+CYT(KP)*UIY(KP,J)+DYT(KP) 1073 : &+ CZT(KP)*UIZ(KP,J)))/3.D0*GEO1 1074 : C &+ (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1075 : 1076 : ZVGT=AIRE(KP)*( 1077 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 1078 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 1079 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 1080 : &+ ALFT(KP)*XMH(KP)*GEO1 1081 : C &+ ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 1082 : 1083 : V2=(UIX(KP,J)*HR(K,J,1) 1084 : $ +UIY(KP,J)*HR(K,J,2) 1085 : $ +UIZ(KP,J)*HR(K,J,3))*DR(KP,I) 1086 : 1087 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2) 1088 : $ +TETAD(KP,J)*ZVGT 1089 : 1090 : 80152 CONTINUE 1091 : 80151 CONTINUE 1092 : 8015 CONTINUE 1093 : 1094 : ENDIF 1095 : 1096 : C 1097 : C Fin de l'accumulation dans SBF. 1098 : C On ajoute ces incréments G. 1099 : C 1100 : DO 8017 I=1,NP 1101 : C*IBMDIR* PREFER VECTOR 1102 : DO 80171 K=KDEB,KFIN 1103 : KP=K-KDEB+1 1104 : NF=IPADL(LE(I,K)) 1105 : G(NF) = G(NF)+SBF(KP,I) 1106 : 80171 CONTINUE 1107 : 8017 CONTINUE 1108 : 1109 : 8001 CONTINUE 1110 : 1111 : C WRITE(IOIMP,*)' G DANS YCTSCL ' 1112 : C WRITE(IOIMP,1984)(M,G(M),M=1,NPTD) 1113 : 1113 : C CALL ARRET(0) 1114 : IPAS=1 1115 : RETURN 1116 : 1002 FORMAT(10(1X,1PE11.4)) 1117 : END 1118 : 1119 : 1120 : 1121 : 1122 : 1123 :
© Cast3M 2003 - All rights reserved.
Disclaimer