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