1 : C ZPSI SOURCE CHAT 05/01/13 04:24:52 5004 2 : SUBROUTINE ZPSI 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 : & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI, 12 : & IPADL,IKOMP,IKAS, 13 : & ALFE,IND1,UN,INDU,NPTU,IPADU, 14 : & TN,QE,IKS, 15 : & HRN,G,ALT,SGT, 16 : & VOLU,COTE,NELZ,ZTE, 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 : -INC CCVQUA4 37 : -INC CCREEL 38 : -INC SMCOORD 39 : -INC CCOPTIO 40 : 41 : C 42 : C Longueur des registres vectoriels de la machine cible 43 : C On prend 64 pour ne pas augmenter la taille des tableaux 44 : C nécessaires à la vectorisation. 45 : C 46 : PARAMETER(LRV=64) 47 : 48 : DIMENSION UN(NPTU,IES),HRN(NPTD),TN(NPTD) 49 : DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*) 50 : DIMENSION ALFE(*),ALT(*) 51 : 52 : DIMENSION IPADL(1),LE(NP,1),IPADU(*) 53 : DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL) 54 : 55 : DIMENSION BF(9,9) 56 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 57 : 58 : DIMENSION COEF(LRV),AIRE(LRV) 59 : DIMENSION AL2(LRV),AH2(LRV),AP2(LRV) 60 : DIMENSION AL(LRV),AH(LRV),AP(LRV) 61 : DIMENSION ALFT(LRV),QT(LRV) 62 : DIMENSION XMB(LRV),XMH(LRV) 63 : DIMENSION CFM(LRV) 64 : DIMENSION CF1(LRV),CF2(LRV),CF3(LRV) 65 : DIMENSION DR(LRV,9) 66 : DIMENSION UM(LRV),UP(LRV) 67 : DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9) 68 : DIMENSION TETAC(LRV,9),TETAD(LRV,9) 69 : DIMENSION UMI(LRV,3) 70 : DIMENSION SBF(LRV,9) 71 : DIMENSION GRADT(LRV,3) 72 : 73 : C? DIMENSION ZTE(LRVH) 74 : DIMENSION ZTE(NPTU) 75 : 76 : INTEGER U(3),D(3),SGN 77 : REAL*8 G(1),ZVGG(4),x(4),y(4),KMAX 78 : REAL*8 KKK(4),PhiSource(LRV),n(3,4) 79 : REAL*8 PhiN3,PhiN4,Phi3NNQ,Phi4NNQ 80 : REAL*8 MINI2,MINI4,MAXI2,MAXI4 81 : 82 : SAVE IPAS,QGGT,Q1,Q2,Q3 83 : 84 : DATA CD/1.D0/ 85 : DATA IPAS/0/ 86 : C********************************************************************* 87 : C 88 : C INITIALISATIONS DIVERSES 89 : C 90 : 91 : NK=K0 92 : C ******** 93 : C * 2D * 94 : C ******** 95 : 96 : IF (IES.EQ.3) GOTO 10 97 : 98 : IAX1=0 99 : IAX2=0 100 : IF (IAXI.EQ.1) IAX2=1 101 : IF (IAXI.EQ.2) IAX1=1 102 : 103 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 104 : 105 : C DIFFERENCES TRIANGLE / QUADRANGLE 106 : QUA4=0.D0 107 : IF (NP.EQ.4) QUA4=1.D0 108 : C 109 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 110 : 111 : C 112 : C Calcul du nombre de paquets de LRV éléments 113 : C 114 : NNN=MOD(NEL,LRV) 115 : IF(NNN.EQ.0) NPACK=NEL/LRV 116 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 117 : KPACKD=1 118 : KPACKF=NPACK 119 : 120 : ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 121 : 122 : DO 7001 KPACK=KPACKD,KPACKF 123 : 124 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 125 : 126 : C 1. Calcul des limites du paquet courant. 127 : 128 : KDEB=1+(KPACK-1)*LRV 129 : KFIN=MIN(NEL,KDEB+LRV-1) 130 : DO 7002 K=KDEB,KFIN 131 : KP=K-KDEB+1 132 : NK=K+K0 133 : NK1=(1-IND1)*(NK-1)+1 134 : ALFT(KP)=ALFE(NK1)+XPETIT 135 : AIRE(KP)=VOLU(NK) 136 : AL(KP)=COTE(NK,1) 137 : AH(KP)=COTE(NK,2) 138 : AL2(KP)=1.D0/AL(KP)/AL(KP) 139 : AH2(KP)=1.D0/AH(KP)/AH(KP) 140 : 7002 CONTINUE 141 : 142 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 143 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 144 : DO 7003 K=KDEB,KFIN 145 : KP=K-KDEB+1 146 : NK=K+K0 147 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 148 : 7003 CONTINUE 149 : ENDIF 150 : 151 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 152 : 153 : C Initialisation des UMI avant accumulation 154 : DO 7005 K=KDEB,KFIN 155 : KP=K-KDEB+1 156 : UMI(KP,1)=XPETIT 157 : UMI(KP,2)=XPETIT 158 : 7005 CONTINUE 159 : 160 : DO 7006 I=1,NP 161 : DO 7006 K=KDEB,KFIN 162 : KP=K-KDEB+1 163 : NF=IPADL(LE(I,K)) 164 : NFU=IPADU(LE(I,K)) 165 : NFU=(1-INDU)*(NFU-1)+1 166 : DR(KP,I)=DRR(I,K) 167 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 168 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 169 : 170 : TETAC(KP,I)=HRN(NF) 171 : TETAD(KP,I)=TN(NF) 172 : 7006 CONTINUE 173 : 174 : * Calcul du FLUX PhiSource du terme source. 175 : 176 : IF (IKOMP.EQ.0) THEN 177 : 178 : DO 70021 K=KDEB,KFIN 179 : KP=K-KDEB+1 180 : NK=K+K0 181 : NKS=(1-IKS)*(NK-1)+1 182 : QT(KP)=QE(NKS) 183 : 70021 CONTINUE 184 : 185 : DO 70062 K=KDEB,KFIN 186 : KP=K-KDEB+1 187 : PhiSource(KP)=-QT(KP)*AIRE(KP) 188 : 70062 CONTINUE 189 : 190 : ELSEIF (IKOMP.EQ.1) THEN 191 : 192 : DO 70023 K=KDEB,KFIN 193 : KP=K-KDEB+1 194 : NK=K+K0 195 : NKS=(1-IKS)*(NK-1)+1 196 : QT(KP)=QE(NKS) 197 : 70023 CONTINUE 198 : 199 : DO 70064 K=KDEB,KFIN 200 : KP=K-KDEB+1 201 : PhiSource(KP)=-QT(KP)*AIRE(KP) 202 : 70064 CONTINUE 203 : 204 : ENDIF 205 : 206 : 207 : DO 7007 K=KDEB,KFIN 208 : KP=K-KDEB+1 209 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 210 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 211 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 212 : UM(KP)=SQRT(UM(KP)) + XPETIT 213 : 214 : 7007 CONTINUE 215 : 216 : ************************************************** 217 : 217 : CALL INITD(ZTE,NPTU,XPETIT) 218 : 219 : IF (IKOMP.EQ.0) THEN 220 : 221 : DO 7008 K=KDEB,KFIN 222 : KP=K-KDEB+1 223 : 224 : DO 23 I=1,NP 225 : 23 ZVGG(I) =0.D0 226 : 227 : ************************************************** 228 : 229 : * NP=3 correspond à PSI sur triangles 230 : * NP=4 correspond à PSI sur quadrangles ,cad NNQ 231 : 232 : * Calcul des Ki:********************************** 233 : 234 : IF (NP.EQ.3) THEN 235 : 236 : DO 1 I=1,NP 237 : 1 KKK(I)=AIRE(KP)*(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2)) 238 : 239 : ELSEIF (NP.EQ.4) THEN 240 : 241 : * Calcul des normales: 242 : 243 : DO 36 I=1,NP 244 : 245 : x(I)=XCOOR((LE(I,K)-1)*(IES+1)+1) 246 : 36 y(I)=XCOOR((LE(I,K)-1)*(IES+1)+2) 247 : 248 : DO 39 I=1,NP 249 : 250 : IF (I.EQ.4) THEN 251 : J1=1 252 : ELSE 253 : J1=I+1 254 : ENDIF 255 : 256 : n(1,I)= (y(I)-y(J1)) 257 : n(2,I)= (x(J1)-x(I)) 258 : 259 : 39 CONTINUE 260 : 261 : * Calcul des Ki : 262 : 263 : DO 17 I=1,NP 264 : 17 KKK(I)=0.5D0*(UMI(KP,1)*n(1,I)+UMI(KP,2)*n(2,I)) 265 : 266 : ENDIF 267 : 268 : 269 : **************************************************** 270 : * Calcul du DT: 271 : **************************************************** 272 : 273 : KMAX=KKK(1) 274 : DO 29 I=2,NP 275 : IF (KKK(I).GT.KMAX) KMAX=KKK(I) 276 : 29 CONTINUE 277 : 278 : DO 27 J=1,NP 279 : IP=IPADU(LE(J,K)) 280 : 281 : ZTE(IP) = ZTE(IP) + KMAX/(2.D0*AIRE(KP)+XPETIT) 282 : 27 CONTINUE 283 : 284 : **************************************************** 285 : 286 : 287 : IF (NP.EQ.4) THEN 288 : 289 : ************************************************** 290 : * Schéma NNQ (Pour des Quadrangles) . 291 : ************************************************** 292 : 293 : * Tests: 294 : ******** 295 : Nd=0 296 : 297 : rrt =KKK(2)+KKK(3) 298 : 299 : IF (rrt.GT.0.D0) THEN 300 : Nd=Nd+1 301 : D(Nd)=1 302 : ENDIF 303 : 304 : rrt =KKK(3)+KKK(4) 305 : 306 : IF (rrt.GT.0.D0) THEN 307 : Nd=Nd+1 308 : D(Nd)=2 309 : ENDIF 310 : 311 : rrt =KKK(1)+KKK(4) 312 : 313 : IF (rrt.GT.0.D0) THEN 314 : Nd=Nd+1 315 : D(Nd)=3 316 : ENDIF 317 : 318 : rrt =KKK(1)+KKK(2) 319 : 320 : IF (rrt.GT.0.D0) THEN 321 : Nd=Nd+1 322 : D(Nd)=4 323 : ENDIF 324 : 325 : IF (Nd.EQ.2) THEN 326 : IF ((D(1)+1).NE.D(2)) THEN 327 : jj = D(1) 328 : D(1) = D(2) 329 : D(2) = jj 330 : ENDIF 331 : ENDIF 332 : 333 : Nu=0 334 : DO 66 I=1,NP 335 : IF (Nd.EQ.2) THEN 336 : IF ((I.NE.D(1)).AND.(I.NE.D(2))) THEN 337 : Nu=Nu+1 338 : U(Nu)=I 339 : ENDIF 340 : ELSEIF (Nd.EQ.1) THEN 341 : IF (I.NE.D(1)) THEN 342 : Nu=Nu+1 343 : U(Nu)=I 344 : ENDIF 345 : ENDIF 346 : 66 CONTINUE 347 : 348 : IF (Nu+Nd.NE.4) THEN 349 : Print *,'Nd=,Nu=',Nd,Nu 350 : ENDIF 351 : 352 : IF (Nu.EQ.2) THEN 353 : IF ((U(1)+1).NE.U(2)) THEN 354 : jj=U(1) 355 : U(1)=U(2) 356 : U(2)=jj 357 : ENDIF 358 : ENDIF 359 : 360 : N1 =U(1) 361 : N2 =U(2) 362 : N3 =D(1) 363 : N4 =D(2) 364 : 365 : IF (Nd.EQ.2) THEN 366 : 367 : MINI2=0.D0 368 : IF (0.D0.GE.KKK(N2)) MINI2=KKK(N2) 369 : 370 : MINI4=0.D0 371 : IF (0.D0.GE.KKK(N4)) MINI4=KKK(N4) 372 : 373 : MAXI2=0.D0 374 : IF (0.D0.LE.KKK(N2)) MAXI2=KKK(N2) 375 : 376 : MAXI4=0.D0 377 : IF (0.D0.LE.KKK(N4)) MAXI4=KKK(N4) 378 : 379 : PhiN3 =( KKK(N1) + MINI2 + MINI4 ) 380 : * *(TETAC(KP,N3)-TETAC(KP,N2)) 381 : * +( MAXI4-MINI2 )*(TETAC(KP,N3)-TETAC(KP,N1)) 382 : 383 : PhiN4 =( KKK(N1) + MINI2 + MINI4 ) 384 : * *(TETAC(KP,N4)-TETAC(KP,N1)) 385 : * +( MAXI2-MINI4 )*(TETAC(KP,N4)-TETAC(KP,N2)) 386 : 387 : ZVGG(N3)=PhiN3 388 : ZVGG(N4)=PhiN4 389 : 390 : r= -ZVGG(N3) / (ZVGG(N4)+XPETIT) 391 : IF (r.GT.0.D0) THEN 392 : IF (r.GT.1.D0) THEN 393 : ZVGG(N3)=ZVGG(N3)+ZVGG(N4) 394 : ZVGG(N4)=0.D0 395 : ELSE 396 : ZVGG(N4)=ZVGG(N3)+ZVGG(N4) 397 : ZVGG(N3)=0.D0 398 : ENDIF 399 : ENDIF 400 : 401 : ELSEIF (ND.EQ.1) THEN 402 : 403 : PhiQ =0.5D0*(UMI(KP,1)*(n(1,1)+n(1,4)) + 404 : * UMI(KP,2)*(n(2,1)+n(2,4)))*(TETAC(KP,3)-TETAC(KP,1))+ 405 : * 0.5D0*(UMI(KP,1)*(n(1,1)+n(1,2)) + 406 : * UMI(KP,2)*(n(2,1)+n(2,2)))*(TETAC(KP,4)-TETAC(KP,2)) 407 : 408 : ZVGG(N3) = PhiQ 409 : ENDIF 410 : 411 : ELSEIF (NP.EQ.3) THEN 412 : 413 : ************************************************** 414 : * Schéma PSI (Min-Mod) Pour des Triangles : 415 : ************************************************** 416 : 417 : * U Noeuds amonts 418 : * D Noeuds avals . 419 : ********************** 420 : 421 : Nd=0 422 : Nu=0 423 : 424 : DO 13 I=1,NP 425 : 426 : IF (KKK(I).GT.0.D0) THEN 427 : Nd=Nd+1 428 : D(Nd)=I 429 : ELSE 430 : Nu=Nu+1 431 : U(Nu)=I 432 : ENDIF 433 : 434 : 13 CONTINUE 435 : 436 : IF (Nd.EQ.1) THEN 437 : 438 : *************************** 439 : * 1 Target : 440 : *************************** 441 : 442 : N1 = D(1) 443 : N2 = U(1) 444 : N3 = U(2) 445 : 446 : ZVGG(N1) = PhiSource(KP) + KKK(1)*TETAC(KP,1)+ 447 : * KKK(2)*TETAC(KP,2) + KKK(3)*TETAC(KP,3) 448 : ZVGG(N2) = 0.D0 449 : ZVGG(N3) = 0.D0 450 : 451 : ELSEIF (Nd.EQ.2) THEN 452 : 453 : *************************** 454 : * 2 Targets : 455 : *************************** 456 : 457 : N1 = D(1) 458 : N2 = D(2) 459 : N3 = U(1) 460 : 461 : ZVGG(N1) = KKK(N1)*(TETAC(KP,N1)-TETAC(KP,N3))-PhiSource(KP)* 462 : * KKK(N1)/(KKK(N3)+XPETIT) 463 : 464 : ZVGG(N2) = KKK(N2)*(TETAC(KP,N2)-TETAC(KP,N3))-PhiSource(KP)* 465 : * KKK(N2)/(KKK(N3)+XPETIT) 466 : 467 : ZVGG(N3) = 0.D0 468 : 469 : r= -ZVGG(N1) / (ZVGG(N2)+XPETIT) 470 : IF (r.GT.0.D0) THEN 471 : IF (r.GT.1.D0) THEN 472 : ZVGG(N1)=ZVGG(N1)+ZVGG(N2) 473 : ZVGG(N2)=0.D0 474 : ELSE 475 : ZVGG(N2)=ZVGG(N1)+ZVGG(N2) 476 : ZVGG(N1)=0.D0 477 : ENDIF 478 : ENDIF 479 : 480 : ENDIF 481 : 482 : ENDIF 483 : 484 : DO 7 I=1,NP 485 : SBF(KP,I)=ZVGG(I) 486 : 7 CONTINUE 487 : 488 : 7008 CONTINUE 489 : ************************************************** 490 : 491 : ELSEIF (IKOMP.EQ.1) THEN 492 : 493 : * Prog à compléter 494 : RETURN 495 : 496 : ENDIF 497 : 498 : * Pas de temps DT : 499 : 500 : DT=100. 501 : 502 : DO 32 K1=KDEB,KFIN 503 : DO 32 J=1,NP 504 : IP=IPADU(LE(J,K1)) 505 : DT1=1.D0 /( ZTE(IP)+XPETIT) 506 : 507 : c IF (Abs(DT1).GE.1000) THEN 508 : c Print *,'DT1 trop grand =',DT1 509 : c ENDIF 510 : 511 : IF (DT1.LE.DT) DT=DT1 512 : 32 CONTINUE 513 : 514 : 515 : DTT1=DT 516 : DTT2=DT 517 : DTT3=0.D0 518 : DIAEL=0.D0 519 : NUEL=0 520 : 521 : * PRINT *, 'PAS DE TEMPS = ', DT 522 : 523 : 524 : **************************************************** 525 : * Contribution du terme diffusif: 526 : 527 : DO 8 I=1,NP 528 : DO 8 J=1,NP 529 : DO 8 K=KDEB,KFIN 530 : KP=K-KDEB+1 531 : 532 : ZVGT=AIRE(KP)*( 533 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 534 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 535 : &+ (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0 536 : & *VGGT(J,I)*QUA4 ) 537 : 538 : 539 : SBF(KP,I)=SBF(KP,I)+ TETAD(KP,J)*ZVGT 540 : 541 : 8 CONTINUE 542 : **************************************************** 543 : 544 : C Fin de l'accumulation dans SBF. 545 : C On ajoute ces incréments G. 546 : 547 : 548 : DO 7017 I=1,NP 549 : DO 7017 K=KDEB,KFIN 550 : KP=K-KDEB+1 551 : NF=IPADL(LE(I,K)) 552 : 553 : G(NF) = G(NF)+SBF(KP,I) 554 : 555 : 7017 CONTINUE 556 : 557 : 558 : 7001 CONTINUE 559 : IPAS=1 560 : 561 : RETURN 562 : 563 : 564 : C ******** 565 : C * 3D * 566 : C ******** 567 : 568 : 10 CONTINUE 569 : 569 : IF (IPAS.EQ.0) CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 570 : CUB8 = 0.D0 571 : IF (NP.EQ.8) CUB8=1.D0 572 : 573 : C 574 : C Calcul du nombre de paquets de LRV éléments 575 : C 576 : NNN=MOD(NEL,LRV) 577 : IF(NNN.EQ.0) NPACK=NEL/LRV 578 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 579 : KPACKD=1 580 : KPACKF=NPACK 581 : C 582 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 583 : C 584 : DO 8001 KPACK=KPACKD,KPACKF 585 : C 586 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 587 : C 588 : C 1. Calcul des limites du paquet courant. 589 : KDEB=1+(KPACK-1)*LRV 590 : KFIN=MIN(NEL,KDEB+LRV-1) 591 : C 592 : DO 8002 K=KDEB,KFIN 593 : KP=K-KDEB+1 594 : NK=K+K0 595 : NK1=(1-IND1)*(NK-1)+1 596 : ALFT(KP)=ALFE(NK1)+XPETIT 597 : AIRE(KP)=VOLU(NK) 598 : AL(KP)=COTE(NK,1) 599 : AH(KP)=COTE(NK,2) 600 : AP(KP)=COTE(NK,3) 601 : 602 : CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ 603 : & AP(KP)*AH(KP)/AL(KP) 604 : C CF1(KP)=AL(KP)*AH(KP)/AP(KP) 605 : C CF2(KP)=AL(KP)*AP(KP)/AH(KP) 606 : C CF3(KP)=AP(KP)*AH(KP)/AL(KP) 607 : XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0 608 : AL2(KP)=1.D0/AL(KP)/AL(KP) 609 : AH2(KP)=1.D0/AH(KP)/AH(KP) 610 : AP2(KP)=1.D0/AP(KP)/AP(KP) 611 : 612 : 8002 CONTINUE 613 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 614 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 615 : DO 8003 K=KDEB,KFIN 616 : KP=K-KDEB+1 617 : NK=K+K0 618 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 619 : 8003 CONTINUE 620 : ENDIF 621 : 622 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 623 : 624 : C Initialisation des UMI avant accumulation 625 : DO 8005 K=KDEB,KFIN 626 : KP=K-KDEB+1 627 : UMI(KP,1)=XPETIT 628 : UMI(KP,2)=XPETIT 629 : UMI(KP,3)=XPETIT 630 : 8005 CONTINUE 631 : 632 : DO 8006 I=1,NP 633 : DO 8006 K=KDEB,KFIN 634 : KP=K-KDEB+1 635 : NF=IPADL(LE(I,K)) 636 : NFU=IPADU(LE(I,K)) 637 : NFU=(1-INDU)*(NFU-1)+1 638 : DR(KP,I)=DRR(I,K) 639 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 640 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 641 : UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I) 642 : 643 : TETAC(KP,I)=HRN(NF) 644 : TETAD(KP,I)=TN(NF) 645 : 8006 CONTINUE 646 : 647 : C 648 : C Initialisation de la variable d'accumulation SBF au terme source 649 : C 650 : C write(6,*)' IKomp,ikas=',IKomp,ikas 651 : C write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU 652 : 653 : IF(IKOMP.EQ.0)THEN 654 : 655 : DO 80021 K=KDEB,KFIN 656 : KP=K-KDEB+1 657 : NK=K+K0 658 : NKS=(1-IKS)*(NK-1)+1 659 : QT(KP)=QE(NKS) 660 : 80021 CONTINUE 661 : 662 : DO 80062 I=1,NP 663 : DO 80062 K=KDEB,KFIN 664 : KP=K-KDEB+1 665 : SBF(KP,I)=-QT(KP)*DR(KP,I) 666 : 80062 CONTINUE 667 : 668 : ELSEIF(IKOMP.EQ.1)THEN 669 : 670 : DO 80023 K=KDEB,KFIN 671 : KP=K-KDEB+1 672 : NK=K+K0 673 : NKS=(1-IKS)*(NK-1)+1 674 : QT(KP)=QE(NKS) 675 : 80023 CONTINUE 676 : 677 : DO 80064 I=1,NP 678 : DO 80064 K=KDEB,KFIN 679 : KP=K-KDEB+1 680 : SBF(KP,I)=-QT(KP)*DR(KP,I) 681 : 80064 CONTINUE 682 : 683 : ENDIF 684 : 685 : DO 8007 K=KDEB,KFIN 686 : KP=K-KDEB+1 687 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 688 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 689 : UMI(KP,3)=UMI(KP,3)/AIRE(KP) 690 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 691 : & +UMI(KP,3)*UMI(KP,3) 692 : UM(KP)=SQRT(UM(KP)) 693 : 8007 CONTINUE 694 : 695 : 696 : ******************************************************** 697 : * Debut de PSI: 698 : 699 : 700 : 701 : 702 : 703 : 704 : 705 : 706 : ******************************************************** 707 : 708 : 709 : C Le coeur du calcul ... 710 : 711 : * IF(IKOMP.EQ.0)THEN 712 : 713 : C DO 8014 I=1,NP 714 : C DO 8014 K=KDEB,KFIN 715 : C KP=K-KDEB+1 716 : 717 : * ELSEIF(IKOMP.EQ.1)THEN 718 : * RETURN 719 : 720 : * ENDIF 721 : 722 : C 723 : C Fin de l'accumulation dans SBF. 724 : C On ajoute ces incréments G. 725 : C 726 : DO 8017 I=1,NP 727 : DO 8017 K=KDEB,KFIN 728 : 729 : KP=K-KDEB+1 730 : NF=IPADL(LE(I,K)) 731 : G(NF) = G(NF)+SBF(KP,I) 732 : 8017 CONTINUE 733 : 734 : 8001 CONTINUE 735 : 736 : 737 : IPAS=1 738 : RETURN 739 : 1002 FORMAT(10(1X,1PE11.4)) 740 : 741 : END 742 : 743 : 744 : 745 : 746 :
© Cast3M 2003 - All rights reserved.
Disclaimer