1 : C YCVIT SOURCE CHAT 05/01/13 04:16:48 5004 2 : SUBROUTINE YCVIT 3 : & (HR,RPG,DRR,LE,NEL,K0,NPTI,NPTS,NPTU,IES,NP,IAXI, 4 : & IPADI,IPADS,IPADU,IKOMP,IKAS, 5 : & RO,IKR, 6 : & COEFF,IK1,RGE,IKG,NELG,TN,IKT,TREF,IKREF, 7 : & AMUT,GN,UN,TK,TE, 8 : & F,GK,GE,GE1, 9 : & VOLU,COTE,NELZ,IDCEN,IMODEL, 10 : & DT,DTT1,DTT2,NUEL,DIAEL,DIAM) 11 : 12 : IMPLICIT INTEGER(I-N) 13 : IMPLICIT REAL*8 (A-H,O-Z) 14 : 15 : C*********************************************************************** 16 : C APPELE PAR NSKE 17 : C CE SP DISCRETISE LES EQUATIONS DE NAVIER STOKES COUPLEES 18 : C AUX DEUX EQUATIONS DU MODELE K - EPSILON 19 : C EN 2D SUR LES ELEMENTS QUA4 ET TRI3 PLAN OU AXI 20 : C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 21 : C LES OPERATEURS SONT "SOUS-INTEGRES" 22 : C 23 : C*** 24 : C 25 : C MODELE K - EPSILON 26 : C 27 : C 28 : C CONSTANTES : CNU=0.09D0 C1=1.41D0 C2=1.92D0 29 : C SGT=0.9D0 SGK=1.D0 SGE=1.3D0 30 : C 31 : C P --> TABLEAU P 32 : C 33 : C LE NU TURBULENT EST CALCULE PAR ELEMENT 34 : C IL N EST PAS TRANSMIS EN ARGUMENT 35 : C 36 : C 37 : C*********************************************************************** 38 : 39 : DIMENSION UN(NPTU,IES),GN(NPTI,IES),TK(NPTI),TE(NPTI) 40 : DIMENSION AMUT(*),RO(*) 41 : DIMENSION COTE(NELZ,IES),VOLU(NELZ) 42 : 43 : C NOMBRE DE NOEUDS MAX PAR ELEMENT NPX 44 : PARAMETER (NPX=9) 45 : 46 : C LONGUEUR DES REGISTRES VECTORIELS DE LA MACHINE CIBLE 47 : 48 : PARAMETER (LRV=64) 49 : C 50 : DIMENSION IPADI(*),IPADS(*),IPADU(*),LE(NP,*) 51 : DIMENSION HR(NEL,NP,IES),RPG(*),DRR(NP,NEL) 52 : DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) 53 : SAVE QGGT,Q1,Q2,Q3 54 : 55 : DIMENSION F(NPTS,IES),GK(NPTS),GE(NPTS),GE1(NPTS) 56 : 57 : -INC CCVQUA4 58 : -INC CCREEL 59 : *- 60 : -INC CCOPTIO 61 : -INC SMCOORD 62 : C* -TABLEAUX ADDITIONNELS POUR LA VECTORISATION ********************** 63 : C 64 : C ---- ILS CONDUISENT A UN SURCOUT DE LRV*(59+25NPX) MOTS. AVEC NPX=9 65 : C ET LRV=128 ON OBTIENT 37KMOTS (0.3 MEGAOCTETS) ET CECI QUELLE 66 : C SOIT LA NOMBRE TOTAL D'ELEMENTS DU MAILLAGE 67 : C ------------------------------------------------------------------ 68 : DIMENSION AIRE(LRV),COEF(LRV),AL(LRV),AH(LRV),AP(LRV) 69 : DIMENSION XMB(LRV),XMH (LRV) 70 : DIMENSION CFM(LRV) 71 : C DIMENSION CF1(LRV),CF2(LRV),CF3(LRV) 72 : C 73 : DIMENSION DR(LRV,NPX) 74 : DIMENSION UIX(LRV,NPX),UIY(LRV,NPX),UIZ(LRV,NPX) 75 : DIMENSION GIX(LRV,NPX),GIY(LRV,NPX),GIZ(LRV,NPX) 76 : DIMENSION AK(LRV,NPX),AE(LRV,NPX) 77 : C 78 : DIMENSION UM(LRV),UMI(LRV,3),AKM(LRV),AEM(LRV),ARO(LRV) 79 : DIMENSION P(LRV),PG(LRV),ANUT(LRV), 80 : ; DUDX(LRV),DUDY(LRV),DUDZ(LRV),DVDX(LRV),DVDY(LRV), 81 : ; DVDZ(LRV),DWDX(LRV),DWDY(LRV),DWDZ(LRV), 82 : ; DTDX(LRV),DTDY(LRV),DTDZ(LRV) 83 : C 84 : DIMENSION C1P(LRV),C1G(LRV),COEFT(LRV),COEFK(LRV),COEFE(LRV) 85 : C 86 : DIMENSION BE1(LRV,NPX) 87 : C 88 : C 89 : C ---- TABLEAUX POUR L'OPTION RAPIDE ------------------------------- 90 : C DANS LA TRIPLE BOUCLE SUR K,I,J ON PEUT SORTIR UN MAXIMUM 91 : C D'OPERATIONS DE LA BOUCLES LA PLUS INTERNE. LES TABLEAUX CI- 92 : C DESSOUS PERMETTENT DE STOCKER LES VALEURS CALCULEES 93 : C ------------------------------------------------------------------ 94 : C*** 95 : DIMENSION COEFF(*),RGE(NELG,IES) 96 : DIMENSION TN(*),TREF(*) 97 : 98 : DIMENSION RGX(LRV),RGY(LRV),RGZ(LRV) 99 : DIMENSION SAF1(LRV,9),SAF2(LRV,9),SAF3(LRV,9) 100 : DIMENSION SAFK(LRV,9),SAFE(LRV,9) 101 : DIMENSION AL2(LRV),AH2(LRV),AP2(LRV) 102 : DIMENSION GRADGX(LRV,3),GRADGY(LRV,3),GRADGZ(LRV,3) 103 : DIMENSION GRADK (LRV,3),GRADE (LRV,3) 104 : DIMENSION UPGX(LRV),UPIGX(LRV,3),UPGY(LRV),UPIGY(LRV,3) 105 : DIMENSION UPGZ(LRV),UPIGZ(LRV,3) 106 : DIMENSION UPK (LRV),UPIK (LRV,3),UPE (LRV),UPIE (LRV,3) 107 : DIMENSION BM(LRV),BPGX(LRV),BPGY(LRV),BPGZ(LRV) 108 : DIMENSION BPK (LRV),BPE (LRV) 109 : DIMENSION GXCXT(LRV),GXCYT(LRV),GXCZT(LRV) 110 : DIMENSION GXCXY(LRV),GXCXZ(LRV),GXCYZ(LRV) 111 : DIMENSION GYCXT(LRV),GYCYT(LRV),GYCZT(LRV) 112 : DIMENSION GYCXY(LRV),GYCXZ(LRV),GYCYZ(LRV) 113 : DIMENSION GZCXT(LRV),GZCYT(LRV),GZCZT(LRV) 114 : DIMENSION GZCXY(LRV),GZCXZ(LRV),GZCYZ(LRV) 115 : DIMENSION GKCXT(LRV),GKCYT(LRV),GKCZT(LRV) 116 : DIMENSION GKCXY(LRV),GKCXZ(LRV),GKCYZ(LRV) 117 : DIMENSION GECXT(LRV),GECYT(LRV),GECZT(LRV) 118 : DIMENSION GECXY(LRV),GECXZ(LRV),GECYZ(LRV) 119 : 120 : DIMENSION CONST(2,7) 121 : 122 : SAVE IPAS 123 : DATA IPAS/0/ 124 : 125 : C Constantes du modèle RNG 126 : C DATA CNU,C1,C2,C3/0.0845D0,1.42D0,1.68D0,0.8D0/ 127 : C DATA SGT,SGK,SGE /0.9D0,0.7179D0,0.7179D0/ 128 : C 129 : C Modèle standard 130 : C DATA CNU,C1,C2,C3,SGT,SGE/0.09D0,1.41D0,1.92D0,0.8D0,0.9D0,1.3D0/ 131 : C 132 : C IMODEL=1 modele Standard IMODEL=2 modele RNG 133 : C 134 : C 135 : C*** Variation des constantes en fonction du modèle 136 : DATA (CONST(1,j),j=1,4) /0.09D0,1.41D0,1.92D0,0.8D0/ 137 : DATA (CONST(1,j),j=5,7) /0.9D0,1.D0,1.3D0/ 138 : DATA (CONST(2,j),j=1,4) /0.0845D0,1.42D0,1.68D0,0.8D0/ 139 : DATA (CONST(2,j),j=5,7) /0.9D0,0.7179D0,0.7179D0/ 140 : C*** 141 : C 142 : CDIR$ VECTOR 143 : CDIR$ IVDEP 144 : 145 : C WRITE(IOIMP,*)' Debut YCVIT' 146 : C WRITE(IOIMP,*)' NELZ,ies=',nelz,ies 147 : C WRITE(IOIMP,1002)cote 148 : C WRITE(IOIMP,*)' Avt precaution ies=',ies 149 : C 150 : C INITIALISATIONS DIVERSES 151 : C 152 : zxzini = sqrt(xpetit) 153 : DTM1=0.D0 154 : 154 : C IF(IPAS.EQ.0)CALL CALHRG(QGGT,IES) 154 : IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) 155 : NK=K0 156 : C ******** 157 : C * 2D * 158 : C ******** 159 : C 160 : C IMODEL=3 161 : C Suite de la definition des constantes 162 : C IF(IMODEL.EQ.3) THEN 163 : C IMO=2 164 : C ELSE 165 : C IMO=1 166 : C ENDIF 167 : C 168 : C (IMODEL=1 modèle Standard ; IMODEL=2 modèle RNG) 169 : CNU=CONST(IMODEL,1) 170 : C1=CONST(IMODEL,2) 171 : C2=CONST(IMODEL,3) 172 : C3=CONST(IMODEL,4) 173 : SGT=CONST(IMODEL,5) 174 : SGK=CONST(IMODEL,6) 175 : SGE=CONST(IMODEL,7) 176 : C 177 : IF(IES.EQ.3)GO TO 10 178 : C 179 : 180 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 181 : 182 : C DIFFERENCES TRIANGLE / QUADRANGLE 183 : QUA4=0.D0 184 : IF(NP.EQ.4)QUA4=1.D0 185 : 186 : C 187 : C K EST LE NUMERO DE L'ELEMENT 188 : C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS 189 : C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS --- 190 : NNN=MOD(NEL,LRV) 191 : IF(NNN.EQ.0) NPACK=NEL/LRV 192 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 193 : KPACKD=1 194 : KPACKF=NPACK 195 : C 196 : C*** BOUCLE SUR LES PAQUETS D'ELEMENTS *** 197 : C 198 : DO 70001 KPACK=KPACKD,KPACKF 199 : C 200 : C --- POUR CHAQUE PAQUET DE LRV ELEMENTS 201 : KDEB=1+(KPACK-1)*LRV 202 : KFIN=MIN(NEL,KDEB+LRV-1) 203 : C WRITE(IOIMP,*)' KPACK=',KPACK 204 : C 205 : DO 70002 K=KDEB,KFIN 206 : KP=K-KDEB+1 207 : NK=K+K0 208 : K1=1+(1-IK1)*(NK-1) 209 : AIRE(KP)=VOLU(NK) 210 : COEF(KP)=COEFF(K1) 211 : AL(KP)=COTE(NK,1) 212 : AH(KP)=COTE(NK,2) 213 : AL2(KP)=1.D0/AL(KP)/AL(KP) 214 : AH2(KP)=1.D0/AH(KP)/AH(KP) 215 : XMB (KP)=(AL(KP)+AH(KP))*0.5D0 216 : 70002 CONTINUE 217 : C WRITE(IOIMP,*)' Apr bcl 70002,,,,,,,,,,,,,,,,,,,,,,,, ' 218 : 219 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 220 : C Initialisation des UMI,UMIT avant accumulation 221 : DO 7005 K=KDEB,KFIN 222 : KP=K-KDEB+1 223 : UMI(KP,1)=zxzini 224 : UMI(KP,2)=zxzini 225 : GRADGX(KP,1)=zxzini 226 : GRADGX(KP,2)=zxzini 227 : GRADGY(KP,1)=zxzini 228 : GRADGY(KP,2)=zxzini 229 : AKM(KP)=zxzini 230 : AEM(KP)=zxzini 231 : UPIK(KP,1)=zxzini 232 : UPIK(KP,2)=zxzini 233 : UPIE(KP,1)=zxzini 234 : UPIE(KP,2)=zxzini 235 : GRADK(KP,1)=zxzini 236 : GRADK(KP,2)=zxzini 237 : GRADE(KP,1)=zxzini 238 : GRADE(KP,2)=zxzini 239 : DUDX(KP)=0.D0 240 : DUDY(KP)=0.D0 241 : DVDX(KP)=0.D0 242 : DVDY(KP)=0.D0 243 : DTDX(KP)=0.D0 244 : DTDY(KP)=0.D0 245 : C1P(KP)=C1 246 : C1G(KP)=0.D0 247 : PG(KP)=0.D0 248 : 7005 CONTINUE 249 : 250 : DO 7006 I=1,NP 251 : DO 6006 K=KDEB,KFIN 252 : KP=K-KDEB+1 253 : NF=IPADI(LE(I,K)) 254 : NFU=IPADU(LE(I,K)) 255 : DR(KP,I)=DRR(I,K) 256 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 257 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 258 : UIX(KP,I)=UN(NFU,1) 259 : UIY(KP,I)=UN(NFU,2) 260 : GIX(KP,I)=GN(NF,1) 261 : GIY(KP,I)=GN(NF,2) 262 : AK(KP,I)=TK(NF) 263 : AE(KP,I)=TE(NF) 264 : GRADGX(KP,1)=GRADGX(KP,1)+HR(K,I,1)*GIX(KP,I) 265 : GRADGX(KP,2)=GRADGX(KP,2)+HR(K,I,2)*GIX(KP,I) 266 : GRADGY(KP,1)=GRADGY(KP,1)+HR(K,I,1)*GIY(KP,I) 267 : GRADGY(KP,2)=GRADGY(KP,2)+HR(K,I,2)*GIY(KP,I) 268 : AKM(KP)=AKM(KP)+TK(NF)*DR(KP,I) 269 : AEM(KP)=AEM(KP)+TE(NF)*DR(KP,I) 270 : GRADK(KP,1)=GRADK(KP,1)+HR(K,I,1)*AK(KP,I) 271 : GRADK(KP,2)=GRADK(KP,2)+HR(K,I,2)*AK(KP,I) 272 : GRADE(KP,1)=GRADE(KP,1)+HR(K,I,1)*AE(KP,I) 273 : GRADE(KP,2)=GRADE(KP,2)+HR(K,I,2)*AE(KP,I) 274 : C *** CALCUL DU TERME DE PRODUCTION **************************** 275 : DUDX(KP) = DUDX(KP)+HR(K,I,1)*UN(NFU,1) 276 : DUDY(KP) = DUDY(KP)+HR(K,I,2)*UN(NFU,1) 277 : DVDX(KP) = DVDX(KP)+HR(K,I,1)*UN(NFU,2) 278 : DVDY(KP) = DVDY(KP)+HR(K,I,2)*UN(NFU,2) 279 : 6006 CONTINUE 280 : 7006 CONTINUE 281 : 282 : C WRITE(IOIMP,*)'Initialisation SAF1,SAF2,SBF ' 283 : C 284 : C Initialisation des variables d'accumulation SAF1,SAF2,SBF 285 : C 286 : IF(IKOMP.EQ.0)THEN 287 : 288 : DO 70060 K=KDEB,KFIN 289 : KP=K-KDEB+1 290 : ARO(KP)=1.D0 291 : 70060 CONTINUE 292 : 293 : IF(IKAS.EQ.2)THEN 294 : DO 70061 I=1,NP 295 : DO 60061 K=KDEB,KFIN 296 : KP=K-KDEB+1 297 : SAF1(KP,I)=0.D0 298 : SAF2(KP,I)=0.D0 299 : SAFK(KP,I)=0.D0 300 : SAFE(KP,I)=0.D0 301 : 60061 CONTINUE 302 : 70061 CONTINUE 303 : ELSEIF(IKAS.EQ.3)THEN 304 : DO 70021 K=KDEB,KFIN 305 : KP=K-KDEB+1 306 : NK=K+K0 307 : NKG=1+(1-IKG)*(NK-1) 308 : RGX(KP)=RGE(NKG,1) 309 : RGY(KP)=RGE(NKG,2) 310 : 70021 CONTINUE 311 : DO 70062 I=1,NP 312 : DO 60062 K=KDEB,KFIN 313 : KP=K-KDEB+1 314 : SAF1(KP,I)=(-RGX(KP))*DR(KP,I) 315 : SAF2(KP,I)=(-RGY(KP))*DR(KP,I) 316 : SAFK(KP,I)=0.D0 317 : SAFE(KP,I)=0.D0 318 : NKG=1+(1-IKG)*(NK-1) 319 : DTDX(KP)=DTDX(KP)+HR(K,I,1)*RGE(NKG,1) 320 : DTDY(KP)=DTDY(KP)+HR(K,I,2)*RGE(NKG,2) 321 : 60062 CONTINUE 322 : 70062 CONTINUE 323 : ELSEIF(IKAS.EQ.5)THEN 324 : DO 70022 K=KDEB,KFIN 325 : KP=K-KDEB+1 326 : NK=K+K0 327 : NKG=1+(1-IKG)*(NK-1) 328 : RGX(KP)=RGE(NKG,1) 329 : RGY(KP)=RGE(NKG,2) 330 : 70022 CONTINUE 331 : DO 70063 I=1,NP 332 : DO 60063 K=KDEB,KFIN 333 : KP=K-KDEB+1 334 : NF=1+(1-IKT)*(IPADS(LE(I,K))-1) 335 : NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1) 336 : SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I) 337 : SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I) 338 : SAFK(KP,I)=0.D0 339 : SAFE(KP,I)=0.D0 340 : DTDX(KP)=DTDX(KP)+HR(K,I,1)*TN(NF) 341 : DTDY(KP)=DTDY(KP)+HR(K,I,2)*TN(NF) 342 : 60063 CONTINUE 343 : 70063 CONTINUE 344 : ENDIF 345 : 346 : ELSEIF(IKOMP.EQ.1)THEN 347 : 348 : DO 70020 K=KDEB,KFIN 349 : KP=K-KDEB+1 350 : NK=K+K0 351 : NKR=1+(1-IKR)*(NK-1) 352 : ARO(KP)=RO(NKR) 353 : 70020 CONTINUE 354 : 355 : IF(IKAS.EQ.4)THEN 356 : DO 70064 I=1,NP 357 : DO 60064 K=KDEB,KFIN 358 : KP=K-KDEB+1 359 : SAF1(KP,I)=0.D0 360 : SAF2(KP,I)=0.D0 361 : SAFK(KP,I)=0.D0 362 : SAFE(KP,I)=0.D0 363 : 60064 CONTINUE 364 : 70064 CONTINUE 365 : ELSEIF(IKAS.EQ.5)THEN 366 : DO 70024 K=KDEB,KFIN 367 : KP=K-KDEB+1 368 : NK=K+K0 369 : NKG=1+(1-IKG)*(NK-1) 370 : RGX(KP)=RGE(NKG,1) 371 : RGY(KP)=RGE(NKG,2) 372 : 70024 CONTINUE 373 : DO 70065 I=1,NP 374 : DO 60065 K=KDEB,KFIN 375 : KP=K-KDEB+1 376 : SAF1(KP,I)=-RGX(KP)*DR(KP,I) 377 : SAF2(KP,I)=-RGY(KP)*DR(KP,I) 378 : SAFK(KP,I)=0.D0 379 : SAFE(KP,I)=0.D0 380 : 60065 CONTINUE 381 : 70065 CONTINUE 382 : ENDIF 383 : ENDIF 384 : 385 : DO 7007 K=KDEB,KFIN 386 : KP=K-KDEB+1 387 : NK=K+K0 388 : AKM(KP)=ABS(AKM(KP))/AIRE(KP)+zxzini 389 : AEM(KP)=ABS(AEM(KP))/AIRE(KP)+zxzini 390 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 391 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 392 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 393 : UM(KP)=SQRT(UM(KP)) 394 : AX=GRADGX(KP,1)*GRADGX(KP,1)+GRADGX(KP,2)*GRADGX(KP,2) 395 : AX=SQRT(AX)+zxzini 396 : GRADGX(KP,1)=GRADGX(KP,1)/AX 397 : GRADGX(KP,2)=GRADGX(KP,2)/AX 398 : AY=GRADGY(KP,1)*GRADGY(KP,1)+GRADGY(KP,2)*GRADGY(KP,2) 399 : AY=SQRT(AY)+zxzini 400 : GRADGY(KP,1)=GRADGY(KP,1)/AY 401 : GRADGY(KP,2)=GRADGY(KP,2)/AY 402 : UPGX(KP)=GRADGX(KP,1)*UMI(KP,1)+GRADGX(KP,2)*UMI(KP,2) 403 : UPIGX(KP,1)=UPGX(KP)*GRADGX(KP,1) 404 : UPIGX(KP,2)=UPGX(KP)*GRADGX(KP,2) 405 : UPGY(KP)=GRADGY(KP,1)*UMI(KP,1)+GRADGY(KP,2)*UMI(KP,2) 406 : UPIGY(KP,1)=UPGY(KP)*GRADGY(KP,1) 407 : UPIGY(KP,2)=UPGY(KP)*GRADGY(KP,2) 408 : 409 : AX=GRADK(KP,1)*GRADK(KP,1)+GRADK(KP,2)*GRADK(KP,2) 410 : AX=SQRT(AX)+zxzini 411 : GRADK(KP,1)=GRADK(KP,1)/AX 412 : GRADK(KP,2)=GRADK(KP,2)/AX 413 : UPK(KP)=GRADK(KP,1)*UMI(KP,1)+GRADK(KP,2)*UMI(KP,2) 414 : UPIK(KP,1)=UPK(KP)*GRADK(KP,1) 415 : UPIK(KP,2)=UPK(KP)*GRADK(KP,2) 416 : 417 : AX=GRADE(KP,1)*GRADE(KP,1)+GRADE(KP,2)*GRADE(KP,2) 418 : AX=SQRT(AX)+zxzini 419 : GRADE(KP,1)=GRADE(KP,1)/AX 420 : GRADE(KP,2)=GRADE(KP,2)/AX 421 : UPE(KP)=GRADE(KP,1)*UMI(KP,1)+GRADE(KP,2)*UMI(KP,2) 422 : UPIE(KP,1)=UPE(KP)*GRADE(KP,1) 423 : UPIE(KP,2)=UPE(KP)*GRADE(KP,2) 424 : 425 : ANUT(KP)=CNU*AKM(KP)*AKM(KP)/AEM(KP) 426 : AMUT(NK)=ARO(KP)*ANUT(KP) 427 : COEFT(KP)=COEF(KP)+AMUT(NK) 428 : COEFK(KP)=COEF(KP)+ANUT(KP)/SGK 429 : COEFE(KP)=COEF(KP)+ANUT(KP)/SGE 430 : C --- 2. CALCUL DU TERME DE PRODUCTION PROPREMENT DIT ---------- 431 : PTEMP=( (DUDY(KP)+DVDX(KP))**2 432 : $ + 2.D0*(DUDX(KP)**2+DVDY(KP)**2) ) 433 : P(KP)=PTEMP 434 : 7007 CONTINUE 435 : C------------ANCIENNE VERSION-------------------- 436 : C IF(IAXI.NE.0) THEN 437 : C DO 70012 K=KDEB,KFIN 438 : C KP=K+1-KDEB 439 : C P(KP)=P(KP)+2.D0*(UMI(KP,1)*UMI(KP,1))/(RPG(K)*RPG(K)) 440 : C70012 CONTINUE 441 : C END IF 442 : C 443 : C IF(IKAS.EQ.5.AND.IKOMP.EQ.0)THEN 444 : C IF(IMODEL.EQ.1)THEN 445 : C DO 70113 K=KDEB,KFIN 446 : C KP=K+1-KDEB 447 : C PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP))/SGT 448 : C IF(PG(KP).GT.0.D0)C1G(KP)=C1 449 : C70113 CONTINUE 450 : C ELSEIF(IMODEL.EQ.2)THEN 451 : C DO 70213 K=KDEB,KFIN 452 : C KP=K+1-KDEB 453 : C PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP))/SGT 454 : C RF =-PG(KP)/(P(KP)+PG(KP)+zxzini) 455 : C IF(PG(KP).GT.0.D0)RF=0.D0 456 : C C1G(KP)=C1*(1.D0+C3*RF) 457 : C C1P(KP)=C1*(1.D0+C3*RF) 458 : C70213 CONTINUE 459 : C ENDIF 460 : C ELSEIF(IKAS.EQ.2.AND.IKOMP.EQ.0)THEN 461 : C IF(IMODEL.EQ.1)THEN 462 : C DO 70114 K=KDEB,KFIN 463 : C KP=K+1-KDEB 464 : C PG(KP)=(DTDX(KP)+DTDY(KP))/SGT 465 : C IF(PG(KP).GT.0.D0)C1G(KP)=C1 466 : C70114 CONTINUE 467 : C ELSEIF(IMODEL.EQ.2)THEN 468 : C DO 70214 K=KDEB,KFIN 469 : C KP=K+1-KDEB 470 : C PG(KP)=(DTDX(KP)+DTDY(KP))/SGT 471 : C RF =-PG(KP)/(P(KP)+PG(KP)+zxzini) 472 : C IF(PG(KP).GT.0.D0)RF=0.D0 473 : C C1G(KP)=C1*(1.D0+C3*RF ) 474 : C C1P(KP)=C1*(1.D0+C3*RF ) 475 : C70214 CONTINUE 476 : C END IF 477 : C END IF 478 : C-------------VERSION POUR RNG--------------- 479 : IF(IKAS.EQ.5.AND.IKOMP.EQ.0)THEN 480 : DO 70213 K=KDEB,KFIN 481 : KP=K+1-KDEB 482 : PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP))/SGT 483 : RF =-PG(KP)/(P(KP)+PG(KP)+zxzini) 484 : IF(PG(KP).GT.0.D0)RF=0.D0 485 : C1G(KP)=C1*(1.D0+C3*RF) 486 : C1P(KP)=C1*(1.D0+C3*RF) 487 : 70213 CONTINUE 488 : ELSEIF(IKAS.EQ.2.AND.IKOMP.EQ.0)THEN 489 : DO 70214 K=KDEB,KFIN 490 : KP=K+1-KDEB 491 : PG(KP)=(DTDX(KP)+DTDY(KP))/SGT 492 : RF =-PG(KP)/(P(KP)+PG(KP)+zxzini) 493 : IF(PG(KP).GT.0.D0)RF=0.D0 494 : C1G(KP)=C1*(1.D0+C3*RF ) 495 : C1P(KP)=C1*(1.D0+C3*RF ) 496 : 70214 CONTINUE 497 : END IF 498 : 499 : IF(IMODEL.EQ.2) THEN 500 : ETA0=4.377D0 501 : BETA0=0.012D0 502 : ETA=(P(KP)**0.5D0)*AKM(KP)/AEM(KP) 503 : RETA= ETA*(1.D0-(ETA/ETA0))/(1.D0+BETA0*(ETA**3.D0)) 504 : C1P(KP)=C1P(KP)-RETA 505 : END IF 506 : C -------------------------------------------------------------- 507 : C *** FIN DU CALCUL DU TERME DE PRODUCTION ********************* 508 : 509 : DO 70074 K=KDEB,KFIN 510 : KP=K-KDEB+1 511 : BMX=UMI(KP,1)/AL(KP) 512 : BMY=UMI(KP,2)/AH(KP) 513 : BM(KP)=BMX*BMX+BMY*BMY 514 : BM(KP)=SQRT(BM(KP))+zxzini 515 : BPGXX=UPIGX(KP,1)/AL(KP) 516 : BPGXY=UPIGX(KP,2)/AH(KP) 517 : BPGX(KP)=BPGXX*BPGXX+BPGXY*BPGXY 518 : BPGX(KP)=SQRT(BPGX(KP))+zxzini 519 : BPGYX=UPIGY(KP,1)/AL(KP) 520 : BPGYY=UPIGY(KP,2)/AH(KP) 521 : BPGY(KP)=BPGYX*BPGYX+BPGYY*BPGYY 522 : BPGY(KP)=SQRT(BPGY(KP))+zxzini 523 : BPKX=UPIK(KP,1)/AL(KP) 524 : BPKY=UPIK(KP,2)/AH(KP) 525 : BPK(KP)=BPKX*BPKX+BPKY*BPKY 526 : BPK(KP)=SQRT(BPK(KP))+zxzini 527 : BPEX=UPIE(KP,1)/AL(KP) 528 : BPEY=UPIE(KP,2)/AH(KP) 529 : BPE(KP)=BPEX*BPEX+BPEY*BPEY 530 : BPE(KP)=SQRT(BPE(KP))+zxzini 531 : 70074 CONTINUE 532 : C 533 : C AVANT DECENTREMENT 534 : C 535 : C ALFA,AKSI,CCU utilisées seulement dans la boucle 7008 536 : IF(IDCEN.EQ.1)THEN 537 : DO 70081 K=KDEB,KFIN 538 : KP=K-KDEB+1 539 : GXCXT(KP)=0.D0 540 : GXCYT(KP)=0.D0 541 : GXCXY(KP)=0.D0 542 : GYCXT(KP)=0.D0 543 : GYCYT(KP)=0.D0 544 : GYCXY(KP)=0.D0 545 : GKCXT(KP)=0.D0 546 : GKCYT(KP)=0.D0 547 : GKCXY(KP)=0.D0 548 : GECXT(KP)=0.D0 549 : GECYT(KP)=0.D0 550 : GECXY(KP)=0.D0 551 : 70081 CONTINUE 552 : 553 : ELSEIF(IDCEN.EQ.2)THEN 554 : DO 7008 K=KDEB,KFIN 555 : KP=K-KDEB+1 556 : HMK=2.D0*UM(KP)/BM(KP) 557 : XMB(KP)=HMK 558 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 559 : AKSI=ALFA/3.D0 560 : IF(ALFA.GT.3.D0)AKSI=1.D0 561 : CCT=AKSI/BM(KP) 562 : 563 : HPK=2.D0*UPGX(KP)/BPGX(KP) 564 : ALFA=UPGX(KP)*HPK/COEFT(KP)/2.D0 565 : AKSI=ALFA/3.D0 566 : IF(ALFA.GT.3.D0)AKSI=1.D0 567 : CCP=AKSI/BPGX(KP) 568 : CPT=CCP-CCT 569 : CC2X=0.D0 570 : IF(CPT.GE.0.D0)CC2X=CPT 571 : 572 : HPK=2.D0*UPGY(KP)/BPGY(KP) 573 : ALFA=UPGY(KP)*HPK/COEFT(KP)/2.D0 574 : AKSI=ALFA/3.D0 575 : IF(ALFA.GT.3.D0)AKSI=1.D0 576 : CCP=AKSI/BPGY(KP) 577 : CPT=CCP-CCT 578 : CC2Y=0.D0 579 : IF(CPT.GE.0.D0)CC2Y=CPT 580 : 581 : HPK=2.D0*UPK(KP)/BPK(KP) 582 : ALFA=UPK(KP)*HPK/COEFK(KP)/2.D0 583 : AKSI=ALFA/3.D0 584 : IF(ALFA.GT.3.D0)AKSI=1.D0 585 : CCP=AKSI/BPK(KP) 586 : CPT=CCP-CCT 587 : CC2K=0.D0 588 : IF(CPT.GE.0.D0)CC2K=CPT 589 : 590 : HPK=2.D0*UPE(KP)/BPE(KP) 591 : ALFA=UPE(KP)*HPK/COEFE(KP)/2.D0 592 : AKSI=ALFA/3.D0 593 : IF(ALFA.GT.3.D0)AKSI=1.D0 594 : CCP=AKSI/BPE(KP) 595 : CPT=CCP-CCT 596 : CC2E=0.D0 597 : IF(CPT.GE.0.D0)CC2E=CPT 598 : 599 : GXCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT 600 : $ +UPIGX(KP,1)*UPIGX(KP,1)*CC2X) 601 : GXCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT 602 : $ +UPIGX(KP,2)*UPIGX(KP,2)*CC2X) 603 : GXCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT 604 : $ +UPIGX(KP,1)*UPIGX(KP,2)*CC2X) 605 : 606 : GYCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT 607 : $ +UPIGY(KP,1)*UPIGY(KP,1)*CC2Y) 608 : GYCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT 609 : $ +UPIGY(KP,2)*UPIGY(KP,2)*CC2Y) 610 : GYCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT 611 : $ +UPIGY(KP,1)*UPIGY(KP,2)*CC2Y) 612 : 613 : GKCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT 614 : $ +UPIK(KP,1)*UPIK(KP,1)*CC2K) 615 : GKCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT 616 : $ +UPIK(KP,2)*UPIK(KP,2)*CC2K) 617 : GKCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT 618 : $ +UPIK(KP,1)*UPIK(KP,2)*CC2K) 619 : 620 : GECXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT 621 : $ +UPIE(KP,1)*UPIE(KP,1)*CC2E) 622 : GECYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT 623 : $ +UPIE(KP,2)*UPIE(KP,2)*CC2E) 624 : GECXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT 625 : $ +UPIE(KP,1)*UPIE(KP,2)*CC2E) 626 : 627 : 7008 CONTINUE 628 : 629 : ELSEIF(IDCEN.EQ.3)THEN 630 : DO 7018 K=KDEB,KFIN 631 : KP=K-KDEB+1 632 : HMK=2.D0*UM(KP)/BM(KP) 633 : XMB(KP)=HMK 634 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 635 : AKSI=ALFA/3.D0 636 : IF(ALFA.GT.3.D0)AKSI=1.D0 637 : CCT=AKSI/BM(KP) 638 : 639 : GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 640 : GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 641 : GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 642 : 643 : GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 644 : GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 645 : GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 646 : 647 : GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 648 : GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 649 : GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 650 : 651 : GECXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 652 : GECYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 653 : GECXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 654 : 655 : 7018 CONTINUE 656 : 657 : ELSEIF(IDCEN.EQ.4)THEN 658 : DT19=DTM1*0.5D0 659 : DO 7009 K=KDEB,KFIN 660 : KP=K-KDEB+1 661 : XMB(KP)=(AL(KP)+AH(KP))/2.D0 662 : GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 663 : GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 664 : GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 665 : GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 666 : GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 667 : GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 668 : GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 669 : GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 670 : GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 671 : GECXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 672 : GECYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 673 : GECXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 674 : 7009 CONTINUE 675 : 676 : ENDIF 677 : 678 : C WRITE(IOIMP,*)' Avt Calcul DT' 679 : C 680 : C AVANT CALCUL DT 681 : C 682 : DO 7010 K=KDEB,KFIN 683 : KP=K-KDEB+1 684 : DT0=DT 685 : DT1X=2.D0/ 686 : & (UMI(KP,1)*UMI(KP,1)/(GXCXT(KP)+COEFT(KP)) 687 : & +UMI(KP,2)*UMI(KP,2)/(GXCYT(KP)+COEFT(KP))) 688 : 689 : DT1Y=2.D0/ 690 : & (UMI(KP,1)*UMI(KP,1)/(GYCXT(KP)+COEFT(KP)) 691 : & +UMI(KP,2)*UMI(KP,2)/(GYCYT(KP)+COEFT(KP))) 692 : DT2X=0.5D0/ 693 : & ((GXCXT(KP)+COEFT(KP))*AL2(KP) 694 : & +(GXCYT(KP)+COEFT(KP))*AH2(KP)) 695 : DT2Y=0.5D0/ 696 : & ((GYCXT(KP)+COEFT(KP))*AL2(KP) 697 : & +(GYCYT(KP)+COEFT(KP))*AH2(KP)) 698 : IF(DT1X.LT.DT)DT=DT1X 699 : IF(DT1Y.LT.DT)DT=DT1Y 700 : IF(DT2X.LT.DT)DT=DT2X 701 : IF(DT2Y.LT.DT)DT=DT2Y 702 : C IF(DT3.LT.DT)DT=DT3 703 : IF(DT.NE.DT0) THEN 704 : DTT1=MIN(DT1X,DT1Y) 705 : DTT2=MIN(DT2X,DT2Y) 706 : C DTT3=DT3 707 : DIAEL=XMB(KP) 708 : NUEL=K 709 : END IF 710 : 7010 CONTINUE 711 : 712 : C WRITE(IOIMP,*)' le coeur du calcul ' 713 : C Le coeur du calcul ... 714 : 715 : IF(IKOMP.EQ.0)THEN 716 : 717 : DO 7014 I=1,NP 718 : DO 70141 J=1,NP 719 : DO 70142 K=KDEB,KFIN 720 : KP=K-KDEB+1 721 : ZVGGX=AIRE(KP)*( 722 : & HR(K,I,1)*HR(K,J,1)*GXCXT(KP) 723 : & + HR(K,I,1)*HR(K,J,2)*GXCXY(KP) 724 : & + HR(K,I,2)*HR(K,J,1)*GXCXY(KP) 725 : & + HR(K,I,2)*HR(K,J,2)*GXCYT(KP) 726 : & + (GXCXT(KP)*AL2(KP)+GXCYT(KP)*AH2(KP)) 727 : $ /12.D0*VGGT(J,I)*QUA4 ) 728 : 729 : ZVGGY=AIRE(KP)*( 730 : & HR(K,I,1)*HR(K,J,1)*GYCXT(KP) 731 : & + HR(K,I,1)*HR(K,J,2)*GYCXY(KP) 732 : & + HR(K,I,2)*HR(K,J,1)*GYCXY(KP) 733 : & + HR(K,I,2)*HR(K,J,2)*GYCYT(KP) 734 : & + (GYCXT(KP)*AL2(KP)+GYCYT(KP)*AH2(KP))/12 735 : $ .D0*VGGT(J,I)*QUA4 ) 736 : 737 : ZVGK =AIRE(KP)*( 738 : & HR(K,I,1)*HR(K,J,1)*GKCXT(KP) 739 : & + HR(K,I,1)*HR(K,J,2)*GKCXY(KP) 740 : & + HR(K,I,2)*HR(K,J,1)*GKCXY(KP) 741 : & + HR(K,I,2)*HR(K,J,2)*GKCYT(KP) 742 : & + (GKCXT(KP)*AL2(KP)+GKCYT(KP)*AH2(KP))/12 743 : $ .D0*VGGT(J,I)*QUA4 ) 744 : 745 : ZVGE =AIRE(KP)*( 746 : & HR(K,I,1)*HR(K,J,1)*GECXT(KP) 747 : & + HR(K,I,1)*HR(K,J,2)*GECXY(KP) 748 : & + HR(K,I,2)*HR(K,J,1)*GECXY(KP) 749 : & + HR(K,I,2)*HR(K,J,2)*GECYT(KP) 750 : & + (GECXT(KP)*AL2(KP)+GECYT(KP)*AH2(KP))/12 751 : $ .D0*VGGT(J,I)*QUA4 ) 752 : 753 : ZVGD=AIRE(KP)*( 754 : & HR(K,I,1)*HR(K,J,1)*COEFT(KP) 755 : & + HR(K,I,2)*HR(K,J,2)*COEFT(KP) 756 : & + (COEFT(KP)*AL2(KP)+COEFT(KP)*AH2(KP))/12 757 : $ .D0*VGGT(J,I)*QUA4 ) 758 : 759 : ZVGDK=AIRE(KP)*( 760 : & HR(K,I,1)*HR(K,J,1)*COEFK(KP) 761 : & + HR(K,I,2)*HR(K,J,2)*COEFK(KP) 762 : & + (COEFK(KP)*AL2(KP)+COEFK(KP)*AH2(KP))/12 763 : $ .D0*VGGT(J,I)*QUA4 ) 764 : 765 : ZVGDE=AIRE(KP)*( 766 : & HR(K,I,1)*HR(K,J,1)*COEFE(KP) 767 : & + HR(K,I,2)*HR(K,J,2)*COEFE(KP) 768 : & + (COEFE(KP)*AL2(KP)+COEFE(KP)*AH2(KP))/12 769 : $ .D0*VGGT(J,I)*QUA4 ) 770 : 771 : Cnon V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I) 772 : 773 : V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I) 774 : $ *HR(K,J,2) 775 : 776 : SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD 777 : $ *UIX(KP,J) 778 : SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD 779 : $ *UIY(KP,J) 780 : 781 : SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK 782 : $ *AK(KP,J) 783 : SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE 784 : $ *AE(KP,J) 785 : 786 : 70142 CONTINUE 787 : 70141 CONTINUE 788 : 7014 CONTINUE 789 : 790 : ELSEIF(IKOMP.EQ.1)THEN 791 : 792 : DO 7015 I=1,NP 793 : DO 70151 J=1,NP 794 : DO 70152 K=KDEB,KFIN 795 : KP=K-KDEB+1 796 : ZVGGX=AIRE(KP)*( 797 : & HR(K,I,1)*HR(K,J,1)*GXCXT(KP) 798 : & + HR(K,I,1)*HR(K,J,2)*GXCXY(KP) 799 : & + HR(K,I,2)*HR(K,J,1)*GXCXY(KP) 800 : & + HR(K,I,2)*HR(K,J,2)*GXCYT(KP) 801 : & + (GXCXT(KP)*AL2(KP)+GXCYT(KP)*AH2(KP))/12 802 : $ .D0*VGGT(J,I)*QUA4 ) 803 : 804 : ZVGGY=AIRE(KP)*( 805 : & HR(K,I,1)*HR(K,J,1)*GYCXT(KP) 806 : & + HR(K,I,1)*HR(K,J,2)*GYCXY(KP) 807 : & + HR(K,I,2)*HR(K,J,1)*GYCXY(KP) 808 : & + HR(K,I,2)*HR(K,J,2)*GYCYT(KP) 809 : & + (GYCXT(KP)*AL2(KP)+GYCYT(KP)*AH2(KP))/12 810 : $ .D0*VGGT(J,I)*QUA4 ) 811 : 812 : ZVGK =AIRE(KP)*( 813 : & HR(K,I,1)*HR(K,J,1)*GKCXT(KP) 814 : & + HR(K,I,1)*HR(K,J,2)*GKCXY(KP) 815 : & + HR(K,I,2)*HR(K,J,1)*GKCXY(KP) 816 : & + HR(K,I,2)*HR(K,J,2)*GKCYT(KP) 817 : & + (GKCXT(KP)*AL2(KP)+GKCYT(KP)*AH2(KP))/12 818 : $ .D0*VGGT(J,I)*QUA4 ) 819 : 820 : ZVGE =AIRE(KP)*( 821 : & HR(K,I,1)*HR(K,J,1)*GECXT(KP) 822 : & + HR(K,I,1)*HR(K,J,2)*GECXY(KP) 823 : & + HR(K,I,2)*HR(K,J,1)*GECXY(KP) 824 : & + HR(K,I,2)*HR(K,J,2)*GECYT(KP) 825 : & + (GECXT(KP)*AL2(KP)+GECYT(KP)*AH2(KP))/12 826 : $ .D0*VGGT(J,I)*QUA4 ) 827 : 828 : ZVGD=AIRE(KP)*( 829 : & HR(K,I,1)*HR(K,J,1)*COEFT(KP) 830 : & + HR(K,I,2)*HR(K,J,2)*COEFT(KP) 831 : & + (COEFT(KP)*AL2(KP)+COEFT(KP)*AH2(KP))/12 832 : $ .D0*VGGT(J,I)*QUA4 ) 833 : 834 : ZVGDK=AIRE(KP)*( 835 : & HR(K,I,1)*HR(K,J,1)*COEFK(KP) 836 : & + HR(K,I,2)*HR(K,J,2)*COEFK(KP) 837 : & + (COEFK(KP)*AL2(KP)+COEFK(KP)*AH2(KP))/12 838 : $ .D0*VGGT(J,I)*QUA4 ) 839 : 840 : ZVGDE=AIRE(KP)*( 841 : & HR(K,I,1)*HR(K,J,1)*COEFE(KP) 842 : & + HR(K,I,2)*HR(K,J,2)*COEFE(KP) 843 : & + (COEFE(KP)*AL2(KP)+COEFE(KP)*AH2(KP))/12 844 : $ .D0*VGGT(J,I)*QUA4 ) 845 : 846 : COEFT3=COEFT(KP)/3.D0 847 : ZVGUU=AIRE(KP)*( HR(K,I,1)*HR(K,J,1) 848 : & + AL2(KP)/12.D0*VGGT(J,I)*QUA4 )*COEFT3 849 : $ *UIX(KP,J) 850 : ZVGUV=AIRE(KP)*( HR(K,I,1)*HR(K,J,2))*COEFT3*UIY(KP 851 : $ ,J) 852 : ZVGVU=AIRE(KP)*( HR(K,I,2)*HR(K,J,1))*COEFT3*UIX(KP 853 : $ ,J) 854 : ZVGVV=AIRE(KP)*( HR(K,I,2)*HR(K,J,2) 855 : & + AH2(KP)/12.D0*VGGT(J,I)*QUA4 )*COEFT3 856 : $ *UIY(KP,J) 857 : 858 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2))*DR(KP 859 : $ ,I) 860 : 861 : SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD 862 : $ *UIX(KP,J)+ ZVGUU + ZVGUV 863 : SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD 864 : $ *UIY(KP,J)+ ZVGVU + ZVGVV 865 : 866 : SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK 867 : $ *AK(KP,J) 868 : SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE 869 : $ *AE(KP,J) 870 : 871 : 70152 CONTINUE 872 : 70151 CONTINUE 873 : 7015 CONTINUE 874 : 875 : ENDIF 876 : 877 : IF(IAXI.NE.0) THEN 878 : DO 7016 I=1,NP 879 : DO 70161 J=1,NP 880 : DO 70162 K=KDEB,KFIN 881 : KP=K-KDEB+1 882 : C IP=LE(I,K) 883 : C XG=XCOOR((IP-1)*(IDIM+1) +1)+zxzini 884 : C R2=1.D0/XG/XG*DR(KP,I) 885 : 886 : R2=1.D0/RPG(K)/RPG(K)*DR(KP,I) 887 : AX=(GXCXT(KP)+GXCYT(KP))*0.5D0 888 : SAF1(KP,I)=SAF1(KP,I)+R2*(COEFT(KP)*UIX(KP,J)+AX 889 : $ *GIX(KP,J)) 890 : 891 : 70162 CONTINUE 892 : 70161 CONTINUE 893 : 7016 CONTINUE 894 : 895 : IF(IKOMP.EQ.1)THEN 896 : 897 : DO 7118 I=1,NP 898 : DO 71181 K=KDEB,KFIN 899 : KP=K-KDEB+1 900 : 901 : R1=1.D0/RPG(K)*DR(KP,I) 902 : SAF1(KP,I)=SAF1(KP,I)+R1*UMI(KP,1)*GIX(KP,I) 903 : SAF2(KP,I)=SAF2(KP,I)+R1*UMI(KP,1)*GIY(KP,I) 904 : 905 : 71181 CONTINUE 906 : 7118 CONTINUE 907 : 908 : ENDIF 909 : 910 : ENDIF 911 : C 912 : C -------------------------------------------------------------- 913 : C 914 : C 915 : DO 7017 I=1,NP 916 : DO 70171 K=KDEB,KFIN 917 : KP=K-KDEB+1 918 : TSK =(ANUT(KP)*(P(KP)+PG(KP))-AEM(KP))*DR(KP,I) 919 : SAFK(KP,I)=SAFK(KP,I)-TSK 920 : 921 : TSE =(CNU*(C1P(KP)*P(KP)+C1G(KP)*PG(KP))*AKM(KP)) 922 : $ *DR(KP,I) 923 : SAFE(KP,I)=SAFE(KP,I)-TSE 924 : 925 : TSE1 =C2*AEM(KP)/AKM(KP)*DR(KP,I) 926 : BE1(KP,I) =TSE1 927 : 70171 CONTINUE 928 : 7017 CONTINUE 929 : 930 : C 931 : C 932 : C --- ECRITURE DANS LES TABLEAUX F(1,*),F(2,*),G,GK,GE,GE1 (SCALAIRE) 933 : DO 70019 I=1,NP 934 : DO 60019 K=KDEB,KFIN 935 : KP=K+1-KDEB 936 : NF=IPADS(LE(I,K)) 937 : F(NF,1)=F(NF,1)-SAF1(KP,I) 938 : F(NF,2)=F(NF,2)-SAF2(KP,I) 939 : GK(NF) =GK(NF) -SAFK(KP,I) 940 : GE1(NF)=GE1(NF)+BE1(KP,I) 941 : GE(NF) =GE (NF)-SAFE(KP,I) 942 : 60019 CONTINUE 943 : 70019 CONTINUE 944 : 70001 CONTINUE 945 : C 946 : C *** FIN DE LA BOUCLE SUR LES PAQUETS D'ELEMENTS *** 947 : C 948 : C WRITE(IOIMP,*)' FIN YCVIT GK=' 949 : C WRITE(IOIMP,1002)(GK(i),i=1,nptd) 950 : 951 : IPAS=1 952 : RETURN 953 : 954 : C ******** 955 : C * 3D * 956 : C ******** 957 : 958 : 10 CONTINUE 959 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 960 : 961 : C DIFFERENCES HEXAHEDRES / AUTRES ELEMENTS 962 : CUB8=0.D0 963 : IF(NP.EQ.8)CUB8=1.D0 964 : 965 : C 966 : C K EST LE NUMERO DE L'ELEMENT 967 : C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS 968 : C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS --- 969 : NNN=MOD(NEL,LRV) 970 : IF(NNN.EQ.0) NPACK=NEL/LRV 971 : IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV 972 : KPACKD=1 973 : KPACKF=NPACK 974 : C 975 : C*** BOUCLE SUR LES PAQUETS D'ELEMENTS *** 976 : C 977 : DO 80001 KPACK=KPACKD,KPACKF 978 : C 979 : C --- POUR CHAQUE PAQUET DE LRV ELEMENTS 980 : KDEB=1+(KPACK-1)*LRV 981 : KFIN=MIN(NEL,KDEB+LRV-1) 982 : C WRITE(IOIMP,*)' KPACK=',KPACK 983 : C 984 : DO 80002 K=KDEB,KFIN 985 : KP=K-KDEB+1 986 : NK=K+K0 987 : K1=1+(1-IK1)*(NK-1) 988 : AIRE(KP)=VOLU(NK) 989 : COEF(KP)=COEFF(K1) 990 : AL(KP)=COTE(NK,1) 991 : AH(KP)=COTE(NK,2) 992 : AP(KP)=COTE(NK,3) 993 : AL2(KP)=1.D0/AL(KP)/AL(KP) 994 : AH2(KP)=1.D0/AH(KP)/AH(KP) 995 : AP2(KP)=1.D0/AP(KP)/AP(KP) 996 : CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ 997 : & AP(KP)*AH(KP)/AL(KP) 998 : C CF1(KP)=AL(KP)*AH(KP)/AP(KP) 999 : C CF2(KP)=AL(KP)*AP(KP)/AH(KP) 1000 : C CF3(KP)=AP(KP)*AH(KP)/AL(KP) 1001 : XMH (KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0 1002 : XMB (KP)=XMH (KP) 1003 : 80002 CONTINUE 1004 : C WRITE(IOIMP,*)' Apr bcl 80002,,,,,,,,,,,,,,,,,,,,,,,, ' 1005 : 1006 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1007 : C Initialisation des UMI,UMIT avant accumulation 1008 : DO 8005 K=KDEB,KFIN 1009 : KP=K-KDEB+1 1010 : UMI(KP,1)=zxzini 1011 : UMI(KP,2)=zxzini 1012 : UMI(KP,3)=zxzini 1013 : GRADGX(KP,1)=zxzini 1014 : GRADGX(KP,2)=zxzini 1015 : GRADGX(KP,3)=zxzini 1016 : GRADGY(KP,1)=zxzini 1017 : GRADGY(KP,2)=zxzini 1018 : GRADGY(KP,3)=zxzini 1019 : GRADGZ(KP,1)=zxzini 1020 : GRADGZ(KP,2)=zxzini 1021 : GRADGZ(KP,3)=zxzini 1022 : AKM(KP)=zxzini 1023 : AEM(KP)=zxzini 1024 : UPIK(KP,1)=zxzini 1025 : UPIK(KP,2)=zxzini 1026 : UPIK(KP,3)=zxzini 1027 : UPIE(KP,1)=zxzini 1028 : UPIE(KP,2)=zxzini 1029 : UPIE(KP,3)=zxzini 1030 : GRADK(KP,1)=zxzini 1031 : GRADK(KP,2)=zxzini 1032 : GRADK(KP,3)=zxzini 1033 : GRADE(KP,1)=zxzini 1034 : GRADE(KP,2)=zxzini 1035 : GRADE(KP,3)=zxzini 1036 : DUDX(KP)=0.D0 1037 : DUDY(KP)=0.D0 1038 : DUDZ(KP)=0.D0 1039 : DVDX(KP)=0.D0 1040 : DVDY(KP)=0.D0 1041 : DVDZ(KP)=0.D0 1042 : DWDX(KP)=0.D0 1043 : DWDY(KP)=0.D0 1044 : DWDZ(KP)=0.D0 1045 : DTDX(KP)=0.D0 1046 : DTDY(KP)=0.D0 1047 : DTDZ(KP)=0.D0 1048 : C1P(KP)=C1 1049 : C1G(KP)=0.D0 1050 : PG(KP)=0.D0 1051 : 8005 CONTINUE 1052 : 1053 : DO 8006 I=1,NP 1054 : DO 5006 K=KDEB,KFIN 1055 : KP=K-KDEB+1 1056 : NF=IPADI(LE(I,K)) 1057 : NFU=IPADU(LE(I,K)) 1058 : DR(KP,I)=DRR(I,K) 1059 : UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I) 1060 : UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I) 1061 : UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I) 1062 : UIX(KP,I)=UN(NFU,1) 1063 : UIY(KP,I)=UN(NFU,2) 1064 : UIZ(KP,I)=UN(NFU,3) 1065 : GIX(KP,I)=GN(NF,1) 1066 : GIY(KP,I)=GN(NF,2) 1067 : GIZ(KP,I)=GN(NF,3) 1068 : AK(KP,I)=TK(NF) 1069 : AE(KP,I)=TE(NF) 1070 : GRADGX(KP,1)=GRADGX(KP,1)+HR(K,I,1)*GIX(KP,I) 1071 : GRADGX(KP,2)=GRADGX(KP,2)+HR(K,I,2)*GIX(KP,I) 1072 : GRADGX(KP,3)=GRADGX(KP,3)+HR(K,I,3)*GIX(KP,I) 1073 : GRADGY(KP,1)=GRADGY(KP,1)+HR(K,I,1)*GIY(KP,I) 1074 : GRADGY(KP,2)=GRADGY(KP,2)+HR(K,I,2)*GIY(KP,I) 1075 : GRADGY(KP,3)=GRADGY(KP,3)+HR(K,I,3)*GIY(KP,I) 1076 : GRADGZ(KP,1)=GRADGZ(KP,1)+HR(K,I,1)*GIZ(KP,I) 1077 : GRADGZ(KP,2)=GRADGZ(KP,2)+HR(K,I,2)*GIZ(KP,I) 1078 : GRADGZ(KP,3)=GRADGZ(KP,3)+HR(K,I,3)*GIZ(KP,I) 1079 : AKM(KP)=AKM(KP)+TK(NF)*DR(KP,I) 1080 : AEM(KP)=AEM(KP)+TE(NF)*DR(KP,I) 1081 : GRADK(KP,1)=GRADK(KP,1)+HR(K,I,1)*AK(KP,I) 1082 : GRADK(KP,2)=GRADK(KP,2)+HR(K,I,2)*AK(KP,I) 1083 : GRADK(KP,3)=GRADK(KP,3)+HR(K,I,3)*AK(KP,I) 1084 : GRADE(KP,1)=GRADE(KP,1)+HR(K,I,1)*AE(KP,I) 1085 : GRADE(KP,2)=GRADE(KP,2)+HR(K,I,2)*AE(KP,I) 1086 : GRADE(KP,3)=GRADE(KP,3)+HR(K,I,3)*AE(KP,I) 1087 : C *** CALCUL DU TERME DE PRODUCTION **************************** 1088 : DUDX(KP) = DUDX(KP)+HR(K,I,1)*UN(NFU,1) 1089 : DUDY(KP) = DUDY(KP)+HR(K,I,2)*UN(NFU,1) 1090 : DUDZ(KP) = DUDZ(KP)+HR(K,I,3)*UN(NFU,1) 1091 : DVDX(KP) = DVDX(KP)+HR(K,I,1)*UN(NFU,2) 1092 : DVDY(KP) = DVDY(KP)+HR(K,I,2)*UN(NFU,2) 1093 : DVDZ(KP) = DVDZ(KP)+HR(K,I,3)*UN(NFU,2) 1094 : DWDX(KP) = DWDX(KP)+HR(K,I,1)*UN(NFU,3) 1095 : DWDY(KP) = DWDY(KP)+HR(K,I,2)*UN(NFU,3) 1096 : DWDZ(KP) = DWDZ(KP)+HR(K,I,3)*UN(NFU,3) 1097 : 5006 CONTINUE 1098 : 8006 CONTINUE 1099 : 1100 : C WRITE(IOIMP,*)'Initialisation SAF1,SAF2,SBF ' 1101 : C 1102 : C Initialisation des variables d'accumulation SAF1,SAF2,SBF 1103 : C 1104 : IF(IKOMP.EQ.0)THEN 1105 : 1106 : DO 80060 K=KDEB,KFIN 1107 : KP=K-KDEB+1 1108 : ARO(KP)=1.D0 1109 : 80060 CONTINUE 1110 : 1111 : IF(IKAS.EQ.2)THEN 1112 : DO 80061 I=1,NP 1113 : DO 50061 K=KDEB,KFIN 1114 : KP=K-KDEB+1 1115 : SAF1(KP,I)=0.D0 1116 : SAF2(KP,I)=0.D0 1117 : SAF3(KP,I)=0.D0 1118 : SAFK(KP,I)=0.D0 1119 : SAFE(KP,I)=0.D0 1120 : 50061 CONTINUE 1121 : 80061 CONTINUE 1122 : ELSEIF(IKAS.EQ.3)THEN 1123 : DO 80021 K=KDEB,KFIN 1124 : KP=K-KDEB+1 1125 : NK=K+K0 1126 : NKG=1+(1-IKG)*(NK-1) 1127 : RGX(KP)=RGE(NKG,1) 1128 : RGY(KP)=RGE(NKG,2) 1129 : RGZ(KP)=RGE(NKG,3) 1130 : 80021 CONTINUE 1131 : DO 80062 I=1,NP 1132 : DO 50062 K=KDEB,KFIN 1133 : KP=K-KDEB+1 1134 : SAF1(KP,I)=(-RGX(KP))*DR(KP,I) 1135 : SAF2(KP,I)=(-RGY(KP))*DR(KP,I) 1136 : SAF3(KP,I)=(-RGZ(KP))*DR(KP,I) 1137 : SAFK(KP,I)=0.D0 1138 : SAFE(KP,I)=0.D0 1139 : NKG=1+(1-IKG)*(NK-1) 1140 : DTDX(KP)=DTDX(KP)+HR(K,I,1)*RGE(NKG,1) 1141 : DTDY(KP)=DTDY(KP)+HR(K,I,2)*RGE(NKG,2) 1142 : DTDZ(KP)=DTDZ(KP)+HR(K,I,3)*RGE(NKG,3) 1143 : 50062 CONTINUE 1144 : 80062 CONTINUE 1145 : ELSEIF(IKAS.EQ.5)THEN 1146 : DO 80022 K=KDEB,KFIN 1147 : KP=K-KDEB+1 1148 : NK=K+K0 1149 : NKG=1+(1-IKG)*(NK-1) 1150 : RGX(KP)=RGE(NKG,1) 1151 : RGY(KP)=RGE(NKG,2) 1152 : RGZ(KP)=RGE(NKG,3) 1153 : 80022 CONTINUE 1154 : DO 80063 I=1,NP 1155 : DO 50063 K=KDEB,KFIN 1156 : KP=K-KDEB+1 1157 : NF=1+(1-IKT)*(IPADS(LE(I,K))-1) 1158 : NFR=1+(1-IKREF)*(IPADS(LE(I,K))-1) 1159 : SAF1(KP,I)=(-RGX(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I) 1160 : SAF2(KP,I)=(-RGY(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I) 1161 : SAF3(KP,I)=(-RGZ(KP)*(TREF(NFR)-TN(NF)))*DR(KP,I) 1162 : SAFK(KP,I)=0.D0 1163 : SAFE(KP,I)=0.D0 1164 : DTDX(KP)=DTDX(KP)+HR(K,I,1)*TN(NF) 1165 : DTDY(KP)=DTDY(KP)+HR(K,I,2)*TN(NF) 1166 : DTDZ(KP)=DTDZ(KP)+HR(K,I,3)*TN(NF) 1167 : 50063 CONTINUE 1168 : 80063 CONTINUE 1169 : ENDIF 1170 : 1171 : ELSEIF(IKOMP.EQ.1)THEN 1172 : 1173 : DO 80020 K=KDEB,KFIN 1174 : KP=K-KDEB+1 1175 : NK=K+K0 1176 : NKR=1+(1-IKR)*(NK-1) 1177 : ARO(KP)=RO(NKR) 1178 : 80020 CONTINUE 1179 : 1180 : IF(IKAS.EQ.4)THEN 1181 : DO 80064 I=1,NP 1182 : DO 50064 K=KDEB,KFIN 1183 : KP=K-KDEB+1 1184 : SAF1(KP,I)=0.D0 1185 : SAF2(KP,I)=0.D0 1186 : SAF3(KP,I)=0.D0 1187 : SAFK(KP,I)=0.D0 1188 : SAFE(KP,I)=0.D0 1189 : 50064 CONTINUE 1190 : 80064 CONTINUE 1191 : ELSEIF(IKAS.EQ.5)THEN 1192 : DO 80024 K=KDEB,KFIN 1193 : KP=K-KDEB+1 1194 : NK=K+K0 1195 : NKG=1+(1-IKG)*(NK-1) 1196 : RGX(KP)=RGE(NKG,1) 1197 : RGY(KP)=RGE(NKG,2) 1198 : RGZ(KP)=RGE(NKG,3) 1199 : 80024 CONTINUE 1200 : DO 80065 I=1,NP 1201 : DO 50065 K=KDEB,KFIN 1202 : KP=K-KDEB+1 1203 : SAF1(KP,I)=-RGX(KP)*DR(KP,I) 1204 : SAF2(KP,I)=-RGY(KP)*DR(KP,I) 1205 : SAF3(KP,I)=-RGZ(KP)*DR(KP,I) 1206 : SAFK(KP,I)=0.D0 1207 : SAFE(KP,I)=0.D0 1208 : 50065 CONTINUE 1209 : 80065 CONTINUE 1210 : ENDIF 1211 : ENDIF 1212 : 1213 : DO 8007 K=KDEB,KFIN 1214 : KP=K-KDEB+1 1215 : NK=K+K0 1216 : AKM(KP)=ABS(AKM(KP))/AIRE(KP)+zxzini 1217 : AEM(KP)=ABS(AEM(KP))/AIRE(KP)+zxzini 1218 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 1219 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 1220 : UMI(KP,3)=UMI(KP,3)/AIRE(KP) 1221 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 1222 : & +UMI(KP,3)*UMI(KP,3) 1223 : UM(KP)=SQRT(UM(KP)) 1224 : AX=GRADGX(KP,1)*GRADGX(KP,1)+GRADGX(KP,2)*GRADGX(KP,2) 1225 : & +GRADGX(KP,3)*GRADGX(KP,3) 1226 : AX=SQRT(AX)+zxzini 1227 : GRADGX(KP,1)=GRADGX(KP,1)/AX 1228 : GRADGX(KP,2)=GRADGX(KP,2)/AX 1229 : GRADGX(KP,3)=GRADGX(KP,3)/AX 1230 : AY=GRADGY(KP,1)*GRADGY(KP,1)+GRADGY(KP,2)*GRADGY(KP,2) 1231 : & +GRADGY(KP,3)*GRADGY(KP,3) 1232 : AY=SQRT(AY)+zxzini 1233 : GRADGY(KP,1)=GRADGY(KP,1)/AY 1234 : GRADGY(KP,2)=GRADGY(KP,2)/AY 1235 : GRADGY(KP,3)=GRADGY(KP,3)/AY 1236 : AZ=GRADGZ(KP,1)*GRADGZ(KP,1)+GRADGZ(KP,2)*GRADGZ(KP,2) 1237 : & +GRADGZ(KP,3)*GRADGZ(KP,3) 1238 : AZ=SQRT(AZ)+zxzini 1239 : GRADGZ(KP,1)=GRADGZ(KP,1)/AZ 1240 : GRADGZ(KP,2)=GRADGZ(KP,2)/AZ 1241 : GRADGZ(KP,3)=GRADGZ(KP,3)/AZ 1242 : 1243 : UPGX(KP)=GRADGX(KP,1)*UMI(KP,1)+GRADGX(KP,2)*UMI(KP,2) 1244 : & +GRADGX(KP,3)*UMI(KP,3) 1245 : UPIGX(KP,1)=UPGX(KP)*GRADGX(KP,1) 1246 : UPIGX(KP,2)=UPGX(KP)*GRADGX(KP,2) 1247 : UPIGX(KP,3)=UPGX(KP)*GRADGX(KP,3) 1248 : UPGY(KP)=GRADGY(KP,1)*UMI(KP,1)+GRADGY(KP,2)*UMI(KP,2) 1249 : & +GRADGY(KP,3)*UMI(KP,3) 1250 : UPIGY(KP,1)=UPGY(KP)*GRADGY(KP,1) 1251 : UPIGY(KP,2)=UPGY(KP)*GRADGY(KP,2) 1252 : UPIGY(KP,3)=UPGY(KP)*GRADGY(KP,3) 1253 : UPGZ(KP)=GRADGZ(KP,1)*UMI(KP,1)+GRADGZ(KP,2)*UMI(KP,2) 1254 : & +GRADGZ(KP,3)*UMI(KP,3) 1255 : UPIGZ(KP,1)=UPGZ(KP)*GRADGZ(KP,1) 1256 : UPIGZ(KP,2)=UPGZ(KP)*GRADGZ(KP,2) 1257 : UPIGZ(KP,3)=UPGZ(KP)*GRADGZ(KP,3) 1258 : 1259 : AX=GRADK(KP,1)*GRADK(KP,1)+GRADK(KP,2)*GRADK(KP,2) 1260 : & +GRADK(KP,3)*GRADK(KP,3) 1261 : AX=SQRT(AX)+zxzini 1262 : GRADK(KP,1)=GRADK(KP,1)/AX 1263 : GRADK(KP,2)=GRADK(KP,2)/AX 1264 : GRADK(KP,3)=GRADK(KP,3)/AX 1265 : UPK(KP)=GRADK(KP,1)*UMI(KP,1)+GRADK(KP,2)*UMI(KP,2) 1266 : & +GRADK(KP,3)*UMI(KP,3) 1267 : UPIK(KP,1)=UPK(KP)*GRADK(KP,1) 1268 : UPIK(KP,2)=UPK(KP)*GRADK(KP,2) 1269 : UPIK(KP,3)=UPK(KP)*GRADK(KP,3) 1270 : 1271 : AX=GRADE(KP,1)*GRADE(KP,1)+GRADE(KP,2)*GRADE(KP,2) 1272 : & +GRADE(KP,3)*GRADE(KP,3) 1273 : AX=SQRT(AX)+zxzini 1274 : GRADE(KP,1)=GRADE(KP,1)/AX 1275 : GRADE(KP,2)=GRADE(KP,2)/AX 1276 : GRADE(KP,3)=GRADE(KP,3)/AX 1277 : UPE(KP)=GRADE(KP,1)*UMI(KP,1)+GRADE(KP,2)*UMI(KP,2) 1278 : & +GRADE(KP,3)*UMI(KP,3) 1279 : UPIE(KP,1)=UPE(KP)*GRADE(KP,1) 1280 : UPIE(KP,2)=UPE(KP)*GRADE(KP,2) 1281 : UPIE(KP,3)=UPE(KP)*GRADE(KP,3) 1282 : 1283 : ANUT(KP)=CNU*AKM(KP)*AKM(KP)/AEM(KP) 1284 : AMUT(NK)=ARO(KP)*ANUT(KP) 1285 : COEFT(KP)=COEF(KP)+AMUT(NK) 1286 : COEFK(KP)=COEF(KP)+ANUT(KP) 1287 : COEFE(KP)=COEF(KP)+ANUT(KP)/SGE 1288 : C --- 2.D0 CALCUL DU TERME DE PRODUCTION PROPREMENT DIT ---------- 1289 : PTEMP=( (DUDY(KP)+DVDX(KP))**2 + (DUDZ(KP)+DWDX(KP))**2 1290 : & + (DVDZ(KP)+DWDY(KP))**2 1291 : & + 2.D0*(DUDX(KP)**2+DVDY(KP)**2+DWDZ(KP)**2) ) 1292 : P(KP)=PTEMP 1293 : 8007 CONTINUE 1294 : 1295 : IF(IKAS.EQ.5.AND.IKOMP.EQ.0)THEN 1296 : IF(IMODEL.EQ.1)THEN 1297 : DO 80113 K=KDEB,KFIN 1298 : KP=K+1-KDEB 1299 : PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP) 1300 : & +RGZ(KP)*DTDZ(KP) )/SGT 1301 : IF(PG(KP).GT.0.D0)C1G(KP)=C1 1302 : 80113 CONTINUE 1303 : ELSEIF(IMODEL.EQ.2)THEN 1304 : DO 80213 K=KDEB,KFIN 1305 : KP=K+1-KDEB 1306 : PG(KP)=(RGX(KP)*DTDX(KP)+RGY(KP)*DTDY(KP) 1307 : & +RGZ(KP)*DTDZ(KP) )/SGT 1308 : RF =-PG(KP)/(P(KP)+PG(KP)+zxzini) 1309 : IF(PG(KP).GT.0.D0)RF=0.D0 1310 : C1G(KP)=C1*(1.D0+C3*RF) 1311 : C1P(KP)=C1*(1.D0+C3*RF) 1312 : 80213 CONTINUE 1313 : ENDIF 1314 : ELSEIF(IKAS.EQ.2.AND.IKOMP.EQ.0)THEN 1315 : IF(IMODEL.EQ.1)THEN 1316 : DO 80114 K=KDEB,KFIN 1317 : KP=K+1-KDEB 1318 : PG(KP)=(DTDX(KP)+DTDY(KP)+DTDZ(KP))/SGT 1319 : IF(PG(KP).GT.0.D0)C1G(KP)=C1 1320 : 80114 CONTINUE 1321 : ELSEIF(IMODEL.EQ.2)THEN 1322 : DO 80214 K=KDEB,KFIN 1323 : KP=K+1-KDEB 1324 : PG(KP)=(DTDX(KP)+DTDY(KP)+DTDZ(KP))/SGT 1325 : RF =-PG(KP)/(P(KP)+PG(KP)+zxzini) 1326 : IF(PG(KP).GT.0.D0)RF=0.D0 1327 : C1G(KP)=C1*(1.D0+C3*RF ) 1328 : C1P(KP)=C1*(1.D0+C3*RF ) 1329 : 80214 CONTINUE 1330 : END IF 1331 : END IF 1332 : 1333 : C -------------------------------------------------------------- 1334 : C *** FIN DU CALCUL DU TERME DE PRODUCTION ********************* 1335 : 1336 : DO 80074 K=KDEB,KFIN 1337 : KP=K-KDEB+1 1338 : BMX=UMI(KP,1)/AL(KP) 1339 : BMY=UMI(KP,2)/AH(KP) 1340 : BMZ=UMI(KP,3)/AP(KP) 1341 : BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ 1342 : BM(KP)=SQRT(BM(KP))+zxzini 1343 : 1344 : BPGXX=UPIGX(KP,1)/AL(KP) 1345 : BPGXY=UPIGX(KP,2)/AH(KP) 1346 : BPGXZ=UPIGX(KP,3)/AP(KP) 1347 : BPGX(KP)=BPGXX*BPGXX+BPGXY*BPGXY+BPGXZ*BPGXZ 1348 : BPGX(KP)=SQRT(BPGX(KP))+zxzini 1349 : BPGYX=UPIGY(KP,1)/AL(KP) 1350 : BPGYY=UPIGY(KP,2)/AH(KP) 1351 : BPGYZ=UPIGY(KP,3)/AP(KP) 1352 : BPGY(KP)=BPGYX*BPGYX+BPGYY*BPGYY+BPGYZ*BPGYZ 1353 : BPGY(KP)=SQRT(BPGY(KP))+zxzini 1354 : BPGZX=UPIGZ(KP,1)/AL(KP) 1355 : BPGZY=UPIGZ(KP,2)/AH(KP) 1356 : BPGZZ=UPIGZ(KP,3)/AP(KP) 1357 : BPGZ(KP)=BPGZX*BPGZX+BPGZY*BPGZY+BPGZZ*BPGZZ 1358 : BPGZ(KP)=SQRT(BPGZ(KP))+zxzini 1359 : 1360 : BPKX=UPIK(KP,1)/AL(KP) 1361 : BPKY=UPIK(KP,2)/AH(KP) 1362 : BPKZ=UPIK(KP,3)/AP(KP) 1363 : BPK(KP)=BPKX*BPKX+BPKY*BPKY+BPKZ*BPKZ 1364 : BPK(KP)=SQRT(BPK(KP))+zxzini 1365 : BPEX=UPIE(KP,1)/AL(KP) 1366 : BPEY=UPIE(KP,2)/AH(KP) 1367 : BPEZ=UPIE(KP,3)/AP(KP) 1368 : BPE(KP)=BPEX*BPEX+BPEY*BPEY+BPEZ*BPEZ 1369 : BPE(KP)=SQRT(BPE(KP))+zxzini 1370 : 80074 CONTINUE 1371 : C 1372 : C AVANT DECENTREMENT 1373 : C 1374 : C ALFA,AKSI,CCU utilisées seulement dans la boucle 7008 1375 : IF(IDCEN.EQ.1)THEN 1376 : DO 80081 K=KDEB,KFIN 1377 : KP=K-KDEB+1 1378 : GXCXT(KP)=0.D0 1379 : GXCYT(KP)=0.D0 1380 : GXCZT(KP)=0.D0 1381 : GXCXY(KP)=0.D0 1382 : GXCXZ(KP)=0.D0 1383 : GXCYZ(KP)=0.D0 1384 : GYCXT(KP)=0.D0 1385 : GYCYT(KP)=0.D0 1386 : GYCZT(KP)=0.D0 1387 : GYCXY(KP)=0.D0 1388 : GYCXZ(KP)=0.D0 1389 : GYCYZ(KP)=0.D0 1390 : GZCXT(KP)=0.D0 1391 : GZCYT(KP)=0.D0 1392 : GZCZT(KP)=0.D0 1393 : GZCXY(KP)=0.D0 1394 : GZCXZ(KP)=0.D0 1395 : GZCYZ(KP)=0.D0 1396 : GKCXT(KP)=0.D0 1397 : GKCYT(KP)=0.D0 1398 : GKCZT(KP)=0.D0 1399 : GKCXY(KP)=0.D0 1400 : GKCXZ(KP)=0.D0 1401 : GKCYZ(KP)=0.D0 1402 : GECXT(KP)=0.D0 1403 : GECYT(KP)=0.D0 1404 : GECZT(KP)=0.D0 1405 : GECXY(KP)=0.D0 1406 : GECXZ(KP)=0.D0 1407 : GECYZ(KP)=0.D0 1408 : 80081 CONTINUE 1409 : 1410 : ELSEIF(IDCEN.EQ.2)THEN 1411 : DO 8008 K=KDEB,KFIN 1412 : KP=K-KDEB+1 1413 : HMK=2.D0*UM(KP)/BM(KP) 1414 : XMB(KP)=HMK 1415 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 1416 : AKSI=ALFA/3.D0 1417 : IF(ALFA.GT.3.D0)AKSI=1.D0 1418 : CCT=AKSI/BM(KP) 1419 : 1420 : HPK=2.D0*UPGX(KP)/BPGX(KP) 1421 : ALFA=UPGX(KP)*HPK/COEFT(KP)/2.D0 1422 : AKSI=ALFA/3.D0 1423 : IF(ALFA.GT.3.D0)AKSI=1.D0 1424 : CCP=AKSI/BPGX(KP) 1425 : CPT=CCP-CCT 1426 : CC2X=0.D0 1427 : IF(CPT.GE.0.D0)CC2X=CPT 1428 : 1429 : HPK=2.D0*UPGY(KP)/BPGY(KP) 1430 : ALFA=UPGY(KP)*HPK/COEFT(KP)/2.D0 1431 : AKSI=ALFA/3.D0 1432 : IF(ALFA.GT.3.D0)AKSI=1.D0 1433 : CCP=AKSI/BPGY(KP) 1434 : CPT=CCP-CCT 1435 : CC2Y=0.D0 1436 : IF(CPT.GE.0.D0)CC2Y=CPT 1437 : 1438 : HPK=2.D0*UPGZ(KP)/BPGZ(KP) 1439 : ALFA=UPGZ(KP)*HPK/COEFT(KP)/2.D0 1440 : AKSI=ALFA/3.D0 1441 : IF(ALFA.GT.3.D0)AKSI=1.D0 1442 : CCP=AKSI/BPGZ(KP) 1443 : CPT=CCP-CCT 1444 : CC2Z=0.D0 1445 : IF(CPT.GE.0.D0)CC2Z=CPT 1446 : 1447 : HPK=2.D0*UPK(KP)/BPK(KP) 1448 : ALFA=UPK(KP)*HPK/COEFK(KP)/2.D0 1449 : AKSI=ALFA/3.D0 1450 : IF(ALFA.GT.3.D0)AKSI=1.D0 1451 : CCP=AKSI/BPK(KP) 1452 : CPT=CCP-CCT 1453 : CC2K=0.D0 1454 : IF(CPT.GE.0.D0)CC2K=CPT 1455 : 1456 : HPK=2.D0*UPE(KP)/BPE(KP) 1457 : ALFA=UPE(KP)*HPK/COEFE(KP)/2.D0 1458 : AKSI=ALFA/3.D0 1459 : IF(ALFA.GT.3.D0)AKSI=1.D0 1460 : CCP=AKSI/BPE(KP) 1461 : CPT=CCP-CCT 1462 : CC2E=0.D0 1463 : IF(CPT.GE.0.D0)CC2E=CPT 1464 : 1465 : GXCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIGX(KP,1)*UPIGX(KP,1 1466 : $ )*CC2X) 1467 : GXCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIGX(KP,2)*UPIGX(KP,2 1468 : $ )*CC2X) 1469 : GXCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIGX(KP,3)*UPIGX(KP,3 1470 : $ )*CC2X) 1471 : GXCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIGX(KP,1)*UPIGX(KP,2 1472 : $ )*CC2X) 1473 : GXCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIGX(KP,1)*UPIGX(KP,3 1474 : $ )*CC2X) 1475 : GXCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIGX(KP,2)*UPIGX(KP,3 1476 : $ )*CC2X) 1477 : 1478 : GYCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIGY(KP,1)*UPIGY(KP,1 1479 : $ )*CC2Y) 1480 : GYCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIGY(KP,2)*UPIGY(KP,2 1481 : $ )*CC2Y) 1482 : GYCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIGY(KP,3)*UPIGY(KP,3 1483 : $ )*CC2Y) 1484 : GYCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIGY(KP,1)*UPIGY(KP,2 1485 : $ )*CC2Y) 1486 : GYCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIGY(KP,1)*UPIGY(KP,3 1487 : $ )*CC2Y) 1488 : GYCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIGY(KP,2)*UPIGY(KP,3 1489 : $ )*CC2Y) 1490 : 1491 : GZCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIGZ(KP,1)*UPIGZ(KP,1 1492 : $ )*CC2Z) 1493 : GZCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIGZ(KP,2)*UPIGZ(KP,2 1494 : $ )*CC2Z) 1495 : GZCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIGZ(KP,3)*UPIGZ(KP,3 1496 : $ )*CC2Z) 1497 : GZCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIGZ(KP,1)*UPIGZ(KP,2 1498 : $ )*CC2Z) 1499 : GZCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIGZ(KP,1)*UPIGZ(KP,3 1500 : $ )*CC2Z) 1501 : GZCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIGZ(KP,2)*UPIGZ(KP,3 1502 : $ )*CC2Z) 1503 : 1504 : GKCXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIK(KP,1)*UPIK(KP,1) 1505 : $ *CC2K) 1506 : GKCYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIK(KP,2)*UPIK(KP,2) 1507 : $ *CC2K) 1508 : GKCZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIK(KP,3)*UPIK(KP,3) 1509 : $ *CC2K) 1510 : GKCXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIK(KP,1)*UPIK(KP,2) 1511 : $ *CC2K) 1512 : GKCXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIK(KP,1)*UPIK(KP,3) 1513 : $ *CC2K) 1514 : GKCYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIK(KP,2)*UPIK(KP,3) 1515 : $ *CC2K) 1516 : 1517 : GECXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPIE(KP,1)*UPIE(KP,1) 1518 : $ *CC2E) 1519 : GECYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPIE(KP,2)*UPIE(KP,2) 1520 : $ *CC2E) 1521 : GECZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPIE(KP,3)*UPIE(KP,3) 1522 : $ *CC2E) 1523 : GECXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPIE(KP,1)*UPIE(KP,2) 1524 : $ *CC2E) 1525 : GECXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPIE(KP,1)*UPIE(KP,3) 1526 : $ *CC2E) 1527 : GECYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPIE(KP,2)*UPIE(KP,3) 1528 : $ *CC2E) 1529 : 1530 : 8008 CONTINUE 1531 : 1532 : ELSEIF(IDCEN.EQ.3)THEN 1533 : DO 8018 K=KDEB,KFIN 1534 : KP=K-KDEB+1 1535 : HMK=2.D0*UM(KP)/BM(KP) 1536 : XMB(KP)=HMK 1537 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 1538 : AKSI=ALFA/3.D0 1539 : IF(ALFA.GT.3.D0)AKSI=1.D0 1540 : CCT=AKSI/BM(KP) 1541 : 1542 : GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 1543 : GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 1544 : GXCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 1545 : GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 1546 : GXCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 1547 : GXCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 1548 : 1549 : GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 1550 : GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 1551 : GYCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 1552 : GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 1553 : GYCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 1554 : GYCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 1555 : 1556 : GZCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 1557 : GZCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 1558 : GZCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 1559 : GZCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 1560 : GZCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 1561 : GZCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 1562 : 1563 : GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 1564 : GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 1565 : GKCZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 1566 : GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 1567 : GKCXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 1568 : GKCYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 1569 : 1570 : GECXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 1571 : GECYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 1572 : GECZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 1573 : GECXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 1574 : GECXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 1575 : GECYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 1576 : 1577 : 8018 CONTINUE 1578 : 1579 : ELSEIF(IDCEN.EQ.4)THEN 1580 : DT19=DTM1*0.5D0 1581 : DO 8009 K=KDEB,KFIN 1582 : KP=K-KDEB+1 1583 : XMB(KP)=(AL(KP)+AH(KP))/2.D0 1584 : 1585 : GXCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 1586 : GXCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 1587 : GXCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 1588 : GXCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 1589 : GXCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 1590 : GXCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 1591 : 1592 : GYCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 1593 : GYCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 1594 : GYCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 1595 : GYCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 1596 : GYCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 1597 : GYCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 1598 : 1599 : GZCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 1600 : GZCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 1601 : GZCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 1602 : GZCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 1603 : GZCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 1604 : GZCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 1605 : 1606 : GKCXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 1607 : GKCYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 1608 : GKCZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 1609 : GKCXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 1610 : GKCXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 1611 : GKCYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 1612 : 1613 : GECXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 1614 : GECYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 1615 : GECZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 1616 : GECXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 1617 : GECXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 1618 : GECYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 1619 : 1620 : 8009 CONTINUE 1621 : 1622 : ENDIF 1623 : 1624 : C WRITE(IOIMP,*)' Avt Calcul DT' 1625 : C 1626 : C AVANT CALCUL DT 1627 : C 1628 : DO 8010 K=KDEB,KFIN 1629 : KP=K-KDEB+1 1630 : DT0=DT 1631 : DT1X=2.D0/ 1632 : & (UMI(KP,1)*UMI(KP,1)/(GXCXT(KP)+COEFT(KP)) 1633 : & +UMI(KP,2)*UMI(KP,2)/(GXCYT(KP)+COEFT(KP)) 1634 : & +UMI(KP,3)*UMI(KP,3)/(GXCZT(KP)+COEFT(KP))) 1635 : DT1Y=2.D0/ 1636 : & (UMI(KP,1)*UMI(KP,1)/(GYCXT(KP)+COEFT(KP)) 1637 : & +UMI(KP,2)*UMI(KP,2)/(GYCYT(KP)+COEFT(KP)) 1638 : & +UMI(KP,3)*UMI(KP,3)/(GYCZT(KP)+COEFT(KP))) 1639 : DT1Z=2.D0/ 1640 : & (UMI(KP,1)*UMI(KP,1)/(GZCXT(KP)+COEFT(KP)) 1641 : & +UMI(KP,2)*UMI(KP,2)/(GZCYT(KP)+COEFT(KP)) 1642 : & +UMI(KP,3)*UMI(KP,3)/(GZCZT(KP)+COEFT(KP))) 1643 : DT2X=0.5D0/ 1644 : & ((GXCXT(KP)+COEFT(KP))*AL2(KP) 1645 : & +(GXCYT(KP)+COEFT(KP))*AH2(KP) 1646 : & +(GXCZT(KP)+COEFT(KP))*AP2(KP)) 1647 : DT2Y=0.5D0/ 1648 : & ((GYCXT(KP)+COEFT(KP))*AL2(KP) 1649 : & +(GYCYT(KP)+COEFT(KP))*AH2(KP) 1650 : & +(GYCZT(KP)+COEFT(KP))*AP2(KP)) 1651 : DT2Z=0.5D0/ 1652 : & ((GZCXT(KP)+COEFT(KP))*AL2(KP) 1653 : & +(GZCYT(KP)+COEFT(KP))*AH2(KP) 1654 : & +(GZCZT(KP)+COEFT(KP))*AP2(KP)) 1655 : IF(DT1X.LT.DT)DT=DT1X 1656 : IF(DT1Y.LT.DT)DT=DT1Y 1657 : IF(DT1Z.LT.DT)DT=DT1Z 1658 : IF(DT2X.LT.DT)DT=DT2X 1659 : IF(DT2Y.LT.DT)DT=DT2Y 1660 : IF(DT2Z.LT.DT)DT=DT2Z 1661 : C IF(DT3.LT.DT)DT=DT3 1662 : IF(DT.NE.DT0) THEN 1663 : DTT1=MIN(DT1X,DT1Y,DT1Z) 1664 : DTT2=MIN(DT2X,DT2Y,DT2Z) 1665 : C DTT3=DT3 1666 : DIAEL=XMB(KP) 1667 : NUEL=K 1668 : END IF 1669 : 8010 CONTINUE 1670 : 1671 : C WRITE(IOIMP,*)' le coeur du calcul ' 1672 : C Le coeur du calcul ... 1673 : 1674 : IF(IKOMP.EQ.0)THEN 1675 : 1676 : DO 8014 I=1,NP 1677 : DO 80141 J=1,NP 1678 : DO 80142 K=KDEB,KFIN 1679 : KP=K-KDEB+1 1680 : C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 1681 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 1682 : ZVGGX=AIRE(KP)*( 1683 : & HR(K,I,1)*HR(K,J,1)*GXCXT(KP) 1684 : & + HR(K,I,2)*HR(K,J,2)*GXCYT(KP) 1685 : & + HR(K,I,3)*HR(K,J,3)*GXCZT(KP) 1686 : & + HR(K,I,1)*HR(K,J,2)*GXCXY(KP) 1687 : & + HR(K,I,2)*HR(K,J,1)*GXCXY(KP) 1688 : & + HR(K,I,1)*HR(K,J,3)*GXCXZ(KP) 1689 : & + HR(K,I,3)*HR(K,J,1)*GXCXZ(KP) 1690 : & + HR(K,I,2)*HR(K,J,3)*GXCYZ(KP) 1691 : & + HR(K,I,3)*HR(K,J,2)*GXCYZ(KP) ) 1692 : & + (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*GEO1 1693 : C &+ (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1694 : 1695 : ZVGGY=AIRE(KP)*( 1696 : & HR(K,I,1)*HR(K,J,1)*GYCXT(KP) 1697 : & + HR(K,I,2)*HR(K,J,2)*GYCYT(KP) 1698 : & + HR(K,I,3)*HR(K,J,3)*GYCZT(KP) 1699 : & + HR(K,I,1)*HR(K,J,2)*GYCXY(KP) 1700 : & + HR(K,I,2)*HR(K,J,1)*GYCXY(KP) 1701 : & + HR(K,I,1)*HR(K,J,3)*GYCXZ(KP) 1702 : & + HR(K,I,3)*HR(K,J,1)*GYCXZ(KP) 1703 : & + HR(K,I,2)*HR(K,J,3)*GYCYZ(KP) 1704 : & + HR(K,I,3)*HR(K,J,2)*GYCYZ(KP) ) 1705 : & + (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*GEO1 1706 : C &+ (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1707 : 1708 : ZVGGZ=AIRE(KP)*( 1709 : & HR(K,I,1)*HR(K,J,1)*GZCXT(KP) 1710 : & + HR(K,I,2)*HR(K,J,2)*GZCYT(KP) 1711 : & + HR(K,I,3)*HR(K,J,3)*GZCZT(KP) 1712 : & + HR(K,I,1)*HR(K,J,2)*GZCXY(KP) 1713 : & + HR(K,I,2)*HR(K,J,1)*GZCXY(KP) 1714 : & + HR(K,I,1)*HR(K,J,3)*GZCXZ(KP) 1715 : & + HR(K,I,3)*HR(K,J,1)*GZCXZ(KP) 1716 : & + HR(K,I,2)*HR(K,J,3)*GZCYZ(KP) 1717 : & + HR(K,I,3)*HR(K,J,2)*GZCYZ(KP) ) 1718 : & + (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*GEO1 1719 : C &+ (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1720 : 1721 : ZVGD=AIRE(KP)*( 1722 : & HR(K,I,1)*HR(K,J,1)*COEFT(KP) 1723 : & + HR(K,I,2)*HR(K,J,2)*COEFT(KP) 1724 : & + HR(K,I,3)*HR(K,J,3)*COEFT(KP) ) 1725 : & + COEFT(KP)*XMH(KP)*GEO1 1726 : C &+ COEFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 1727 : 1728 : ZVGK =AIRE(KP)*( 1729 : & HR(K,I,1)*HR(K,J,1)*GKCXT(KP) 1730 : & + HR(K,I,2)*HR(K,J,2)*GKCYT(KP) 1731 : & + HR(K,I,3)*HR(K,J,3)*GKCZT(KP) 1732 : & + HR(K,I,1)*HR(K,J,2)*GKCXY(KP) 1733 : & + HR(K,I,2)*HR(K,J,1)*GKCXY(KP) 1734 : & + HR(K,I,1)*HR(K,J,3)*GKCXZ(KP) 1735 : & + HR(K,I,3)*HR(K,J,1)*GKCXZ(KP) 1736 : & + HR(K,I,2)*HR(K,J,3)*GKCYZ(KP) 1737 : & + HR(K,I,3)*HR(K,J,2)*GKCYZ(KP) ) 1738 : & + (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*GEO1 1739 : C &+ (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1740 : 1741 : ZVGE =AIRE(KP)*( 1742 : & HR(K,I,1)*HR(K,J,1)*GECXT(KP) 1743 : & + HR(K,I,2)*HR(K,J,2)*GECYT(KP) 1744 : & + HR(K,I,3)*HR(K,J,3)*GECZT(KP) 1745 : & + HR(K,I,1)*HR(K,J,2)*GECXY(KP) 1746 : & + HR(K,I,2)*HR(K,J,1)*GECXY(KP) 1747 : & + HR(K,I,1)*HR(K,J,3)*GECXZ(KP) 1748 : & + HR(K,I,3)*HR(K,J,1)*GECXZ(KP) 1749 : & + HR(K,I,2)*HR(K,J,3)*GECYZ(KP) 1750 : & + HR(K,I,3)*HR(K,J,2)*GECYZ(KP) ) 1751 : & + (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*GEO1 1752 : C &+ (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1753 : 1754 : ZVGDK=AIRE(KP)*( 1755 : & HR(K,I,1)*HR(K,J,1)*COEFK(KP) 1756 : & + HR(K,I,2)*HR(K,J,2)*COEFK(KP) 1757 : & + HR(K,I,3)*HR(K,J,3)*COEFK(KP) ) 1758 : & + COEFK(KP)*XMH(KP)*GEO1 1759 : C &+ COEFK(KP)*XMH(KP)*QGGT(J,I)*CUB8 1760 : 1761 : ZVGDE=AIRE(KP)*( 1762 : & HR(K,I,1)*HR(K,J,1)*COEFE(KP) 1763 : & + HR(K,I,2)*HR(K,J,2)*COEFE(KP) 1764 : & + HR(K,I,3)*HR(K,J,3)*COEFE(KP) ) 1765 : & + COEFE(KP)*XMH(KP)*GEO1 1766 : C &+ COEFE(KP)*XMH(KP)*QGGT(J,I)*CUB8 1767 : 1768 : Cnon V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2)+UIZ(KP,I)*HR(K,J,3)) 1769 : C & *DR(KP,I) 1770 : 1771 : V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I) 1772 : $ *HR(K,J,2)+UMI(KP,3)*DR(KP,I)*HR(K,J,3) 1773 : 1774 : SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD 1775 : $ *UIX(KP,J) 1776 : SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD 1777 : $ *UIY(KP,J) 1778 : SAF3(KP,I)=SAF3(KP,I)+(V2+ZVGGZ)*GIZ(KP,J)+ZVGD 1779 : $ *UIZ(KP,J) 1780 : 1781 : SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK 1782 : $ *AK(KP,J) 1783 : SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE 1784 : $ *AE(KP,J) 1785 : 1786 : 80142 CONTINUE 1787 : 80141 CONTINUE 1788 : 8014 CONTINUE 1789 : 1790 : ELSEIF(IKOMP.EQ.1)THEN 1791 : C-> 1792 : 1793 : DO 8015 I=1,NP 1794 : DO 80151 J=1,NP 1795 : DO 80152 K=KDEB,KFIN 1796 : KP=K-KDEB+1 1797 : C GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8 1798 : GEO1=CFM(KP)*QGGT(J,I)*CUB8 1799 : ZVGGX=AIRE(KP)*( 1800 : & HR(K,I,1)*HR(K,J,1)*GXCXT(KP) 1801 : & + HR(K,I,2)*HR(K,J,2)*GXCYT(KP) 1802 : & + HR(K,I,3)*HR(K,J,3)*GXCZT(KP) 1803 : & + HR(K,I,1)*HR(K,J,2)*GXCXY(KP) 1804 : & + HR(K,I,2)*HR(K,J,1)*GXCXY(KP) 1805 : & + HR(K,I,1)*HR(K,J,3)*GXCXZ(KP) 1806 : & + HR(K,I,3)*HR(K,J,1)*GXCXZ(KP) 1807 : & + HR(K,I,2)*HR(K,J,3)*GXCYZ(KP) 1808 : & + HR(K,I,3)*HR(K,J,2)*GXCYZ(KP) ) 1809 : & + (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*GEO1 1810 : C &+ (GXCXT(KP)+GXCYT(KP)+GXCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1811 : 1812 : ZVGGY=AIRE(KP)*( 1813 : & HR(K,I,1)*HR(K,J,1)*GYCXT(KP) 1814 : & + HR(K,I,2)*HR(K,J,2)*GYCYT(KP) 1815 : & + HR(K,I,3)*HR(K,J,3)*GYCZT(KP) 1816 : & + HR(K,I,1)*HR(K,J,2)*GYCXY(KP) 1817 : & + HR(K,I,2)*HR(K,J,1)*GYCXY(KP) 1818 : & + HR(K,I,1)*HR(K,J,3)*GYCXZ(KP) 1819 : & + HR(K,I,3)*HR(K,J,1)*GYCXZ(KP) 1820 : & + HR(K,I,2)*HR(K,J,3)*GYCYZ(KP) 1821 : & + HR(K,I,3)*HR(K,J,2)*GYCYZ(KP) ) 1822 : & + (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*GEO1 1823 : C &+ (GYCXT(KP)+GYCYT(KP)+GYCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1824 : 1825 : ZVGGZ=AIRE(KP)*( 1826 : & HR(K,I,1)*HR(K,J,1)*GZCXT(KP) 1827 : & + HR(K,I,2)*HR(K,J,2)*GZCYT(KP) 1828 : & + HR(K,I,3)*HR(K,J,3)*GZCZT(KP) 1829 : & + HR(K,I,1)*HR(K,J,2)*GZCXY(KP) 1830 : & + HR(K,I,2)*HR(K,J,1)*GZCXY(KP) 1831 : & + HR(K,I,1)*HR(K,J,3)*GZCXZ(KP) 1832 : & + HR(K,I,3)*HR(K,J,1)*GZCXZ(KP) 1833 : & + HR(K,I,2)*HR(K,J,3)*GZCYZ(KP) 1834 : & + HR(K,I,3)*HR(K,J,2)*GZCYZ(KP) ) 1835 : & + (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*GEO1 1836 : C &+ (GZCXT(KP)+GZCYT(KP)+GZCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1837 : 1838 : ZVGD=AIRE(KP)*( 1839 : & HR(K,I,1)*HR(K,J,1)*COEFT(KP) 1840 : & + HR(K,I,2)*HR(K,J,2)*COEFT(KP) 1841 : & + HR(K,I,3)*HR(K,J,3)*COEFT(KP) ) 1842 : & + COEFT(KP)*GEO1 1843 : C &+ COEFT(KP)*XMH(KP)*GEO1 1844 : C &+ COEFT(KP)*XMH(KP)*QGGT(J,I)*CUB8 1845 : 1846 : ZVGK =AIRE(KP)*( 1847 : & HR(K,I,1)*HR(K,J,1)*GKCXT(KP) 1848 : & + HR(K,I,2)*HR(K,J,2)*GKCYT(KP) 1849 : & + HR(K,I,3)*HR(K,J,3)*GKCZT(KP) 1850 : & + HR(K,I,1)*HR(K,J,2)*GKCXY(KP) 1851 : & + HR(K,I,2)*HR(K,J,1)*GKCXY(KP) 1852 : & + HR(K,I,1)*HR(K,J,3)*GKCXZ(KP) 1853 : & + HR(K,I,3)*HR(K,J,1)*GKCXZ(KP) 1854 : & + HR(K,I,2)*HR(K,J,3)*GKCYZ(KP) 1855 : & + HR(K,I,3)*HR(K,J,2)*GKCYZ(KP) ) 1856 : & + (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*GEO1 1857 : C &+ (GKCXT(KP)+GKCYT(KP)+GKCZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1858 : 1859 : ZVGE =AIRE(KP)*( 1860 : & HR(K,I,1)*HR(K,J,1)*GECXT(KP) 1861 : & + HR(K,I,2)*HR(K,J,2)*GECYT(KP) 1862 : & + HR(K,I,3)*HR(K,J,3)*GECZT(KP) 1863 : & + HR(K,I,1)*HR(K,J,2)*GECXY(KP) 1864 : & + HR(K,I,2)*HR(K,J,1)*GECXY(KP) 1865 : & + HR(K,I,1)*HR(K,J,3)*GECXZ(KP) 1866 : & + HR(K,I,3)*HR(K,J,1)*GECXZ(KP) 1867 : & + HR(K,I,2)*HR(K,J,3)*GECYZ(KP) 1868 : & + HR(K,I,3)*HR(K,J,2)*GECYZ(KP) ) 1869 : & + (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*GEO1 1870 : C &+ (GECXT(KP)+GECYT(KP)+GECZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8 1871 : 1872 : ZVGDK=AIRE(KP)*( 1873 : & HR(K,I,1)*HR(K,J,1)*COEFK(KP) 1874 : & + HR(K,I,2)*HR(K,J,2)*COEFK(KP) 1875 : & + HR(K,I,3)*HR(K,J,3)*COEFK(KP) ) 1876 : & + COEFK(KP)*GEO1 1877 : 1878 : ZVGDE=AIRE(KP)*( 1879 : & HR(K,I,1)*HR(K,J,1)*COEFE(KP) 1880 : & + HR(K,I,2)*HR(K,J,2)*COEFE(KP) 1881 : & + HR(K,I,3)*HR(K,J,3)*COEFE(KP) ) 1882 : & + COEFE(KP)*GEO1 1883 : C &+ COEFE(KP)*XMH(KP)*QGGT(J,I)*CUB8 1884 : 1885 : COEFT3=COEFT(KP)/3.D0 1886 : ZVGUU=(AIRE(KP)*(HR(K,I,1)*HR(K,J,1))+GEO1)*COEFT3 1887 : $ *UIX(KP,J) 1888 : ZVGUV=AIRE(KP)*( HR(K,I,1)*HR(K,J,2))*COEFT3*UIY(KP 1889 : $ ,J) 1890 : ZVGUW=AIRE(KP)*( HR(K,I,1)*HR(K,J,3))*COEFT3*UIZ(KP 1891 : $ ,J) 1892 : ZVGVU=AIRE(KP)*( HR(K,I,2)*HR(K,J,1))*COEFT3*UIX(KP 1893 : $ ,J) 1894 : ZVGVV=(AIRE(KP)*(HR(K,I,2)*HR(K,J,2))+GEO1)*COEFT3 1895 : $ *UIY(KP,J) 1896 : ZVGVW=AIRE(KP)*( HR(K,I,2)*HR(K,J,3))*COEFT3*UIZ(KP 1897 : $ ,J) 1898 : ZVGWU=AIRE(KP)*( HR(K,I,3)*HR(K,J,1))*COEFT3*UIX(KP 1899 : $ ,J) 1900 : ZVGWV=AIRE(KP)*( HR(K,I,3)*HR(K,J,2))*COEFT3*UIY(KP 1901 : $ ,J) 1902 : ZVGWW=(AIRE(KP)*(HR(K,I,3)*HR(K,J,3))+GEO1)*COEFT3 1903 : $ *UIZ(KP,J) 1904 : 1905 : V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+UIZ(KP 1906 : $ ,J)*HR(K,J,3))*DR(KP,I) 1907 : 1908 : C V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)*HR(K,J,2) 1909 : 1910 : SAF1(KP,I)=SAF1(KP,I)+(V2+ZVGGX)*GIX(KP,J)+ZVGD 1911 : $ *UIX(KP,J)+ ZVGUU + ZVGUV + ZVGUW 1912 : SAF2(KP,I)=SAF2(KP,I)+(V2+ZVGGY)*GIY(KP,J)+ZVGD 1913 : $ *UIY(KP,J)+ ZVGVU + ZVGUV + ZVGVW 1914 : SAF3(KP,I)=SAF3(KP,I)+(V2+ZVGGZ)*GIZ(KP,J)+ZVGD 1915 : $ *UIZ(KP,J)+ ZVGWU + ZVGWV + ZVGWW 1916 : 1917 : SAFK(KP,I)=SAFK(KP,I)+(V2+ZVGK)*AK(KP,J)+ZVGDK 1918 : $ *AK(KP,J) 1919 : SAFE(KP,I)=SAFE(KP,I)+(V2+ZVGE)*AE(KP,J)+ZVGDE 1920 : $ *AE(KP,J) 1921 : 1922 : 80152 CONTINUE 1923 : 80151 CONTINUE 1924 : 8015 CONTINUE 1925 : 1926 : ENDIF 1927 : 1928 : C 1929 : C -------------------------------------------------------------- 1930 : C 1931 : C 1932 : DO 8017 I=1,NP 1933 : DO 80171 K=KDEB,KFIN 1934 : KP=K-KDEB+1 1935 : TSK =(ANUT(KP)*(P(KP)+PG(KP))-AEM(KP))*DR(KP,I) 1936 : SAFK(KP,I)=SAFK(KP,I)-TSK 1937 : 1938 : TSE =(CNU*(C1P(KP)*P(KP)+C1G(KP)*PG(KP))*AKM(KP)) 1939 : $ *DR(KP,I) 1940 : SAFE(KP,I)=SAFE(KP,I)-TSE 1941 : 1942 : TSE1 =C2*AEM(KP)/AKM(KP)*DR(KP,I) 1943 : BE1(KP,I) =TSE1 1944 : 80171 CONTINUE 1945 : 8017 CONTINUE 1946 : 1947 : C 1948 : C 1949 : C --- ECRITURE DANS LES TABLEAUX F(1,*),F(2,*),G,GK,GE,GE1 (SCALAIRE) 1950 : DO 80019 I=1,NP 1951 : DO 50019 K=KDEB,KFIN 1952 : KP=K+1-KDEB 1953 : NF=IPADS(LE(I,K)) 1954 : F(NF,1)=F(NF,1)-SAF1(KP,I) 1955 : F(NF,2)=F(NF,2)-SAF2(KP,I) 1956 : F(NF,3)=F(NF,3)-SAF3(KP,I) 1957 : GK(NF) =GK(NF) -SAFK(KP,I) 1958 : GE1(NF)=GE1(NF)+BE1(KP,I) 1959 : GE(NF) =GE (NF)-SAFE(KP,I) 1960 : 50019 CONTINUE 1961 : 80019 CONTINUE 1962 : 80001 CONTINUE 1963 : C 1964 : C *** FIN DE LA BOUCLE SUR LES PAQUETS D'ELEMENTS *** 1965 : C 1966 : C WRITE(IOIMP,*)' FIN YCVIT GK=' 1967 : C WRITE(IOIMP,1002)(GK(i),i=1,nptd) 1968 : 1969 : IPAS=1 1970 : RETURN 1971 : 1972 : 1002 FORMAT(10(1X,1PE11.4)) 1973 : END 1974 : 1975 : 1976 : 1977 : 1978 : 1979 : 1980 :
© Cast3M 2003 - All rights reserved.
Disclaimer