1 : C YCTSCL SOURCE CHAT 05/01/13 04:16:35 5004 2 : SUBROUTINE YCTSCL 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,IES,NP,IAXI, 12 : & IPADI,IPADS,IPADF,IKOMP,IKAS, 13 : & ALFE,IND1,UN,INDU,NPTS,TN,NPTD,QE,IKS, 14 : & HRN,G,NPTI, 15 : & ALT,SGT, 16 : & VOLU,COTE,NELZ,IDCEN,IPG, 17 : & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,FN) 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 : C 39 : C Longueur des registres vectoriels de la machine cible 40 : C On prend 64 pour ne pas augmenter la taille des tableaux 41 : C nécessaires à la vectorisation. 42 : C 43 : PARAMETER(LRV=64) 44 : 45 : DIMENSION UN(NPTS,IES),HRN(NPTI),TN(NPTD) 46 : DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*) 47 : DIMENSION ALFE(*),ALT(*) 48 : 49 : DIMENSION IPADI(*),IPADS(*),IPADF(*),LE(NP,1) 50 : DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL) 51 : 52 : DIMENSION BF(9,9) 53 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 54 : 55 : DIMENSION AIRE(LRV) 56 : DIMENSION AL(LRV),AH(LRV),AP(LRV) 57 : DIMENSION ALFT(LRV),QT(LRV) 58 : DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9) 59 : DIMENSION TETAC(LRV,9),TETAD(LRV,9),TETA(LRV,9) 60 : DIMENSION UMI(LRV,3) 61 : DIMENSION SBF(LRV,9) 62 : DIMENSION WT(LRV,9),CHGLD(LRV),CHGLP(LRV) 63 : 64 : REAL*8 G(NPTS),FN(NP,*) 65 : 66 : SAVE IPAS,QGGT,Q1,Q2,Q3 67 : 68 : DATA CD/1.D0/ 69 : 70 : DATA IPAS/0/ 71 : C************************************************************************ 72 : C 73 : C INITIALISATIONS DIVERSES 74 : C 75 : ZERMA=XPETIT 76 : 77 : NK=K0 78 : C ******** 79 : C * 2D * 80 : C ******** 81 : 82 : 83 : IF(IES.EQ.3)GO TO 10 84 : 85 : IAX1=0 86 : IAX2=0 87 : IF(IAXI.EQ.1)IAX2=1 88 : IF(IAXI.EQ.2)IAX1=1 89 : 90 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 91 : 92 : C DIFFERENCES TRIANGLE / QUADRANGLE 93 : QUA4=0.D0 94 : IF(NP.EQ.4)QUA4=1.D0 95 : C 96 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 97 : 98 : C 99 : C Calcul du nombre de paquets de LRV éléments 100 : C 101 : NNN=MOD(NEL,LRV) 102 : IF(NNN.EQ.0) NPACK=NEL/LRV 103 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 104 : KPACKD=1 105 : KPACKF=NPACK 106 : C 107 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 108 : C 109 : DO 7001 KPACK=KPACKD,KPACKF 110 : C 111 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 112 : C 113 : C 1. Calcul des limites du paquet courant. 114 : KDEB=1+(KPACK-1)*LRV 115 : KFIN=MIN(NEL,KDEB+LRV-1) 116 : C 117 : DO 7002 K=KDEB,KFIN 118 : KP=K-KDEB+1 119 : NK=K+K0 120 : NK1=(1-IND1)*(NK-1)+1 121 : ALFT(KP)=ALFE(NK1)+ZERMA 122 : AIRE(KP)=VOLU(NK) 123 : AL(KP)=COTE(NK,1) 124 : AH(KP)=COTE(NK,2) 125 : 7002 CONTINUE 126 : 127 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 128 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 129 : DO 7003 K=KDEB,KFIN 130 : KP=K-KDEB+1 131 : NK=K+K0 132 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 133 : 7003 CONTINUE 134 : ENDIF 135 : 136 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 137 : 138 : DO 7006 I=1,NP 139 : DO 7006 K=KDEB,KFIN 140 : KP=K-KDEB+1 141 : NF=IPADF(LE(I,K)) 142 : NU=IPADS(LE(I,K)) 143 : NI=IPADI(LE(I,K)) 144 : NFU=(1-INDU)*(NU-1)+1 145 : UIX(KP,I)=UN(NFU,1) 146 : UIY(KP,I)=UN(NFU,2) 147 : TETAC(KP,I)=HRN(NI) 148 : TETAD(KP,I)=TN(NF) 149 : 7006 CONTINUE 150 : 150 : CALL KSUPG1(WT,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV, 151 : &HRN,IPADI,UN,ALFT,NPTS,NEL,NP,DRR,HR,FN, 152 : &AIRE,AL,AH,AP,IDCEN,IPADS,LE,QUA4,IKOMP, 153 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 154 : C 155 : C Initialisation de la variable d'accumulation SBF au terme source 156 : C 157 : 158 : IF(IKOMP.EQ.0)THEN 159 : 160 : DO 70021 K=KDEB,KFIN 161 : KP=K-KDEB+1 162 : NK=K+K0 163 : NKS=(1-IKS)*(NK-1)+1 164 : QT(KP)=QE(NKS) 165 : 70021 CONTINUE 166 : 167 : IF(IPG.EQ.0)THEN 168 : DO 70062 I=1,NP 169 : DO 70062 K=KDEB,KFIN 170 : KP=K-KDEB+1 171 : SBF(KP,I)=-QT(KP)*DRR(I,K) 172 : 70062 CONTINUE 173 : ELSE 174 : DO 71062 I=1,NP 175 : DO 71062 K=KDEB,KFIN 176 : KP=K-KDEB+1 177 : SBF(KP,I)=-QT(KP)*WT(KP,I) 178 : 71062 CONTINUE 179 : ENDIF 180 : 181 : ELSEIF(IKOMP.EQ.1)THEN 182 : 183 : DO 70023 K=KDEB,KFIN 184 : KP=K-KDEB+1 185 : NK=K+K0 186 : NKS=(1-IKS)*(NK-1)+1 187 : QT(KP)=QE(NKS) 188 : 70023 CONTINUE 189 : 190 : IF(IPG.EQ.0)THEN 191 : DO 70064 I=1,NP 192 : DO 70064 K=KDEB,KFIN 193 : KP=K-KDEB+1 194 : SBF(KP,I)=-QT(KP)*DRR(I,K) 195 : 70064 CONTINUE 196 : ELSE 197 : DO 71064 I=1,NP 198 : DO 71064 K=KDEB,KFIN 199 : KP=K-KDEB+1 200 : SBF(KP,I)=-QT(KP)*WT(KP,I) 201 : 71064 CONTINUE 202 : ENDIF 203 : 204 : ENDIF 205 : 206 : C Le coeur du calcul ... 207 : 208 : IF(IKOMP.EQ.0)THEN 209 : 210 : DO 7014 I=1,NP 211 : DO 7014 J= 1,NP 212 : DO 7014 K=KDEB,KFIN 213 : KP=K-KDEB+1 214 : ZVGG=AIRE(KP)*CHGLP(KP)*VGGT(J,I) 215 : 216 : ZVGT=AIRE(KP)*( 217 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 218 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 219 : &+ CHGLD(KP)*VGGT(J,I) ) 220 : 221 : V2=(UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2))*WT(KP,I) 222 : 223 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 224 : 225 : 7014 CONTINUE 226 : 227 : ELSEIF(IKOMP.EQ.1)THEN 228 : 229 : DO 7015 I=1,NP 230 : DO 7015 J= 1,NP 231 : DO 7015 K=KDEB,KFIN 232 : KP=K-KDEB+1 233 : ZVGG=AIRE(KP)*CHGLP(KP)*VGGT(J,I) 234 : 235 : ZVGT=AIRE(KP)*( 236 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 237 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 238 : &+ CHGLD(KP)*VGGT(J,I) ) 239 : 240 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2))*WT(KP,I) 241 : 242 : SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT 243 : 244 : 7015 CONTINUE 245 : 246 : ENDIF 247 : C 248 : C Fin de l'accumulation dans SBF. 249 : C On ajoute ces incréments G. 250 : C 251 : DO 7017 I=1,NP 252 : DO 7017 K=KDEB,KFIN 253 : KP=K-KDEB+1 254 : NF=IPADS(LE(I,K)) 255 : G(NF) = G(NF)-SBF(KP,I) 256 : 7017 CONTINUE 257 : 258 : 7001 CONTINUE 259 : 260 : C WRITE(6,*)' G DANS YCTSCL ' 261 : C WRITE(6,1984)(M,G(M),M=1,NPTS) 262 : 1984 FORMAT(7(1X,I4,2X,1PE11.4)) 263 : 263 : C CALL ARRET(0) 264 : IPAS=1 265 : RETURN 266 : 267 : C ******** 268 : C * 3D * 269 : C ******** 270 : 271 : 10 CONTINUE 272 : 272 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 273 : CUB8=0.D0 274 : IF(NP.EQ.8)CUB8=1.D0 275 : 276 : C 277 : C Calcul du nombre de paquets de LRV éléments 278 : C 279 : NNN=MOD(NEL,LRV) 280 : IF(NNN.EQ.0) NPACK=NEL/LRV 281 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 282 : KPACKD=1 283 : KPACKF=NPACK 284 : C 285 : C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS ********** 286 : C 287 : DO 8001 KPACK=KPACKD,KPACKF 288 : C 289 : C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS ======= 290 : C 291 : C 1. Calcul des limites du paquet courant. 292 : KDEB=1+(KPACK-1)*LRV 293 : KFIN=MIN(NEL,KDEB+LRV-1) 294 : C 295 : DO 8002 K=KDEB,KFIN 296 : KP=K-KDEB+1 297 : NK=K+K0 298 : NK1=(1-IND1)*(NK-1)+1 299 : ALFT(KP)=ALFE(NK1)+ZERMA 300 : AIRE(KP)=VOLU(NK) 301 : AL(KP)=COTE(NK,1) 302 : AH(KP)=COTE(NK,2) 303 : AP(KP)=COTE(NK,3) 304 : 8002 CONTINUE 305 : IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR. 306 : &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN 307 : DO 8003 K=KDEB,KFIN 308 : KP=K-KDEB+1 309 : NK=K+K0 310 : ALFT(KP)=ALFT(KP)+ALT(NK)/SGT 311 : 8003 CONTINUE 312 : ENDIF 313 : 314 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 315 : 316 : DO 8006 I=1,NP 317 : DO 8006 K=KDEB,KFIN 318 : KP=K-KDEB+1 319 : NF=IPADF(LE(I,K)) 320 : NU=IPADS(LE(I,K)) 321 : NI=IPADI(LE(I,K)) 322 : NFU=(1-INDU)*(NU-1)+1 323 : UIX(KP,I)=UN(NFU,1) 324 : UIY(KP,I)=UN(NFU,2) 325 : UIZ(KP,I)=UN(NFU,3) 326 : TETAC(KP,I)=HRN(NI) 327 : TETAD(KP,I)=TN(NF) 328 : 8006 CONTINUE 329 : 329 : CALL KSUPG1(WT,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV, 330 : &HRN,IPADI,UN,ALFT,NPTS,NEL,NP,DRR,HR,FN, 331 : &AIRE,AL,AH,AP,IDCEN,IPADS,LE,CUB8,IKOMP, 332 : &DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 333 : 334 : C 335 : C Initialisation de la variable d'accumulation SBF au terme source 336 : C M 337 : IF(IKOMP.EQ.0)THEN 338 : 339 : DO 80021 K=KDEB,KFIN 340 : KP=K-KDEB+1 341 : NK=K+K0 342 : NKS=(1-IKS)*(NK-1)+1 343 : QT(KP)=QE(NKS) 344 : 80021 CONTINUE 345 : 346 : IF(IPG.EQ.0)THEN 347 : DO 80062 I=1,NP 348 : DO 80062 K=KDEB,KFIN 349 : KP=K-KDEB+1 350 : SBF(KP,I)=-QT(KP)*DRR(I,K) 351 : 80062 CONTINUE 352 : ELSE 353 : DO 81062 I=1,NP 354 : DO 81062 K=KDEB,KFIN 355 : KP=K-KDEB+1 356 : SBF(KP,I)=-QT(KP)*WT(KP,I) 357 : 81062 CONTINUE 358 : ENDIF 359 : 360 : ELSEIF(IKOMP.EQ.1)THEN 361 : 362 : DO 80023 K=KDEB,KFIN 363 : KP=K-KDEB+1 364 : NK=K+K0 365 : NKS=(1-IKS)*(NK-1)+1 366 : QT(KP)=QE(NKS) 367 : 80023 CONTINUE 368 : 369 : IF(IPG.EQ.0)THEN 370 : DO 80064 I=1,NP 371 : DO 80064 K=KDEB,KFIN 372 : KP=K-KDEB+1 373 : SBF(KP,I)=-QT(KP)*DRR(I,K) 374 : 80064 CONTINUE 375 : ELSE 376 : DO 81064 I=1,NP 377 : DO 81064 K=KDEB,KFIN 378 : KP=K-KDEB+1 379 : SBF(KP,I)=-QT(KP)*WT(KP,I) 380 : 81064 CONTINUE 381 : ENDIF 382 : 383 : ENDIF 384 : 385 : C Le coeur du calcul ... 386 : 387 : IF(IKOMP.EQ.0)THEN 388 : 389 : DO 8014 I=1,NP 390 : DO 8014 J= 1,NP 391 : DO 8014 K=KDEB,KFIN 392 : KP=K-KDEB+1 393 : 394 : ZVGG=AIRE(KP)*CHGLP(KP)*QGGT(J,I) 395 : 396 : ZVGT=AIRE(KP)*( 397 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 398 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 399 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 400 : &+ CHGLD(KP)*QGGT(J,I) 401 : 402 : V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)+UMI(KP,3)*HR(K,J,3) 403 : 404 : SBF(KP,I)=SBF(KP,I) 405 : & +TETAC(KP,J)*(ZVGG+V2*WT(KP,I))+TETAD(KP,J)*ZVGT 406 : 407 : 8014 CONTINUE 408 : 409 : ELSEIF(IKOMP.EQ.1)THEN 410 : 411 : DO 8015 I=1,NP 412 : DO 8015 J= 1,NP 413 : DO 8015 K=KDEB,KFIN 414 : 415 : KP=K-KDEB+1 416 : 417 : ZVGG=AIRE(KP)*CHGLP(KP)*QGGT(J,I) 418 : 419 : ZVGT=AIRE(KP)*( 420 : & HR(K,I,1)*HR(K,J,1)*ALFT(KP) 421 : &+ HR(K,I,2)*HR(K,J,2)*ALFT(KP) 422 : &+ HR(K,I,3)*HR(K,J,3)*ALFT(KP) ) 423 : &+ CHGLD(KP)*QGGT(J,I) 424 : 425 : V2=UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+UIZ(KP,J)*HR(K,J,3) 426 : 427 : SBF(KP,I)=SBF(KP,I) 428 : & +TETAC(KP,J)*(ZVGG+V2*WT(KP,I))+TETAD(KP,J)*ZVGT 429 : 430 : 8015 CONTINUE 431 : 432 : ENDIF 433 : 434 : C 435 : C Fin de l'accumulation dans SBF. 436 : C On ajoute ces incréments G. 437 : C 438 : DO 8017 I=1,NP 439 : DO 8017 K=KDEB,KFIN 440 : KP=K-KDEB+1 441 : NF=IPADS(LE(I,K)) 442 : G(NF) = G(NF)-SBF(KP,I) 443 : 8017 CONTINUE 444 : 445 : 8001 CONTINUE 446 : 447 : C WRITE(6,*)' G DANS YCTSCL ' 448 : C WRITE(6,1984)(M,G(M),M=1,NPTS) 449 : 449 : C CALL ARRET(0) 450 : IPAS=1 451 : RETURN 452 : 1002 FORMAT(10(1X,1PE11.4)) 453 : END 454 : 455 : 456 : 457 : 458 : 459 : 460 : 461 :
© Cast3M 2003 - All rights reserved.
Disclaimer