1 : C KSUPG1 SOURCE CHAT 05/01/13 01:09:14 5004 2 : SUBROUTINE KSUPG1(WTILDE,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV, 3 : & TN,IPADI,UN,COEFT,NPTD,NEL,NP,DRR,HR,FN, 4 : & AIRE,AL,AH,AP,IDCEN,IPADU,LE,CHGL,IKOMP, 5 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8 (A-H,O-Z) 8 : -INC CCOPTIO 9 : -INC CCREEL 10 : DIMENSION IPADU(*),LE(NP,1),IPADI(*) 11 : DIMENSION UN(NPTD,IDIM),TN(*) 12 : DIMENSION AIRE(LRV),AL(LRV),AH(LRV),AP(LRV),WTILDE(LRV,NP) 13 : DIMENSION COEFT(LRV),CHGLD(LRV),CHGLP(LRV) 14 : DIMENSION HR(NEL,NP,IDIM),DRR(NP,NEL),FN(NP,1) 15 : PARAMETER (LRV1=64) 16 : DIMENSION UM(LRV1),UMI(LRV,3),UMJ(LRV1,3) 17 : DIMENSION GRAD(LRV1,3) 18 : DIMENSION UP(LRV1),UPI(LRV1,3) 19 : DIMENSION BM(LRV1),BP(LRV1) 20 : DIMENSION TO1(LRV1),TO2(LRV1) 21 : DIMENSION CXT(LRV1),CYT(LRV1),CZT(LRV1) 22 : DIMENSION CXY(LRV1),CXZ(LRV1),CYZ(LRV1),CFM(LRV1) 23 : DIMENSION AL2(LRV1),AH2(LRV1),AP2(LRV1),XMB(LRV1) 24 : DATA IKV/1/,NPG/1/ 25 : CC0=1.D0 26 : xpeti2 = sqrt(XPETIT) 27 : IF(IKOMP.EQ.1)CC0=0.2D0 28 : IF(LRV.NE.LRV1)STOP 29 : DO 7005 N=1,IDIM 30 : DO 70051 K=KDEB,KFIN 31 : KP=K-KDEB+1 32 : UMI(KP,N)=xpeti2 33 : UMJ(KP,N)=xpeti2 34 : UPI(KP,N)=xpeti2 35 : GRAD(KP,N)=xpeti2 36 : 70051 CONTINUE 37 : 7005 CONTINUE 38 : IF(IKV.EQ.1)THEN 39 : DO 70060 N=1,IDIM 40 : DO 70061 I=1,NP 41 : C*IBMDIR* PREFER VECTOR 42 : DO 70062 K=KDEB,KFIN 43 : KP=K-KDEB+1 44 : NF=IPADU(LE(I,K)) 45 : NG=IPADI(LE(I,K)) 46 : UMI(KP,N)=UMI(KP,N)+UN(NF,N)*DRR(I,K) 47 : UMJ(KP,N)=UMJ(KP,N)+UN(NF,N)*DRR(I,K) 48 : GRAD(KP,N)=GRAD(KP,N)+HR(K,I,N)*TN(NG) 49 : 70062 CONTINUE 50 : 70061 CONTINUE 51 : 70060 CONTINUE 52 : DO 7017 K=KDEB,KFIN 53 : KP=K-KDEB+1 54 : UMJ(KP,1)=UMJ(KP,1)/AIRE(KP) 55 : UMJ(KP,2)=UMJ(KP,2)/AIRE(KP) 56 : UM(KP)=UMJ(KP,1)*UMJ(KP,1)+UMJ(KP,2)*UMJ(KP,2) 57 : UM(KP)=SQRT(UM(KP)) 58 : 7017 CONTINUE 59 : ELSE 60 : DO 70160 N=1,IDIM 61 : DO 70161 I=1,NP 62 : C*IBMDIR* PREFER VECTOR 63 : DO 70162 K=KDEB,KFIN 64 : KP=K-KDEB+1 65 : NF=IPADU(LE(I,K)) 66 : NG=IPADI(LE(I,K)) 67 : UMI(KP,N)=UMI(KP,N)+UN(NF,N)*DRR(I,K) 68 : GRAD(KP,N)=GRAD(KP,N)+HR(K,I,N)*TN(NG) 69 : U=0.D0 70 : DO 7026 L=1,NPG 71 : U=U+UN(NF,N)*FN(I,L) 72 : 7026 CONTINUE 73 : UMJ(KP,N)=UMJ(KP,N)+U 74 : 70162 CONTINUE 75 : 70161 CONTINUE 76 : 70160 CONTINUE 77 : DO 7027 K=KDEB,KFIN 78 : KP=K-KDEB+1 79 : UM(KP)=UMJ(KP,1)*UMJ(KP,1)+UMJ(KP,2)*UMJ(KP,2) 80 : UM(KP)=SQRT(UM(KP)) 81 : 7027 CONTINUE 82 : ENDIF 83 : C *********** 84 : C * 2 D * 85 : C *********** 86 : IF(IDIM.EQ.2)THEN 87 : QUA4=CHGL 88 : DO 7004 K=KDEB,KFIN 89 : KP=K-KDEB+1 90 : AL2(KP)=1.D0/AL(KP)/AL(KP) 91 : AH2(KP)=1.D0/AH(KP)/AH(KP) 92 : XMB(KP)=(AL(KP)+AH(KP))/2.D0 93 : 7004 CONTINUE 94 : DO 7007 K=KDEB,KFIN 95 : KP=K-KDEB+1 96 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 97 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 98 : AX=GRAD(KP,1)*GRAD(KP,1)+GRAD(KP,2)*GRAD(KP,2) 99 : AX=SQRT(AX)+xpeti2 100 : GRAD(KP,1)=GRAD(KP,1)/AX 101 : GRAD(KP,2)=GRAD(KP,2)/AX 102 : UP(KP)=(GRAD(KP,1)*UMJ(KP,1)+GRAD(KP,2)*UMJ(KP,2)) 103 : UPI(KP,1)=UP(KP)*GRAD(KP,1) 104 : UPI(KP,2)=UP(KP)*GRAD(KP,2) 105 : 7007 CONTINUE 106 : DO 70074 K=KDEB,KFIN 107 : KP=K-KDEB+1 108 : BMX=UMJ(KP,1)/AL(KP) 109 : BMY=UMJ(KP,2)/AH(KP) 110 : BM(KP)=BMX*BMX+BMY*BMY 111 : BM(KP)=SQRT(BM(KP))+xpeti2 112 : BPX=UPI(KP,1)/AL(KP) 113 : BPY=UPI(KP,2)/AH(KP) 114 : BP(KP)=BPX*BPX+BPY*BPY 115 : BP(KP)=SQRT(BP(KP))+xpeti2 116 : 70074 CONTINUE 117 : C 118 : C DECENTREMENT 119 : C 120 : IF(IDCEN.EQ.1)THEN 121 : C CENTREE 122 : DO 70081 K=KDEB,KFIN 123 : KP=K-KDEB+1 124 : TO1(KP)=0.D0 125 : TO2(KP)=0.D0 126 : CXT(KP)=0.D0 127 : CYT(KP)=0.D0 128 : CXY(KP)=0.D0 129 : 70081 CONTINUE 130 : ELSEIF(IDCEN.EQ.2)THEN 131 : C SUPGDC 132 : DO 7008 K=KDEB,KFIN 133 : KP=K-KDEB+1 134 : HMK=2.D0*UM(KP)/BM(KP) 135 : XMB(KP)=HMK 136 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 137 : AKSI=ALFA/3.D0 138 : IF(ALFA.GT.3.D0)AKSI=1.D0 139 : CCT=AKSI/BM(KP) 140 : HPK=2.D0*UP(KP)/BP(KP) 141 : ALFA=UP(KP)*HPK/COEFT(KP)/2.D0 142 : AKSI=ALFA/3.D0 143 : IF(ALFA.GT.3.D0)AKSI=1.D0 144 : CCP=AKSI/BP(KP) 145 : CPT=CCP-CCT 146 : CC2=0.D0 147 : IF(CPT.GE.0.D0)CC2=CPT 148 : CC2=CC2*CC0 149 : TO1(KP)=CCT 150 : TO2(KP)=CC2 151 : CXT(KP)=(UMJ(KP,1)*UMJ(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2) 152 : CYT(KP)=(UMJ(KP,2)*UMJ(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2) 153 : CXY(KP)=(UMJ(KP,1)*UMJ(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2) 154 : 7008 CONTINUE 155 : ELSEIF(IDCEN.EQ.3)THEN 156 : C SUPG 157 : DO 7018 K=KDEB,KFIN 158 : KP=K-KDEB+1 159 : HMK=2.D0*UM(KP)/BM(KP) 160 : XMB(KP)=HMK 161 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 162 : AKSI=ALFA/3.D0 163 : IF(ALFA.GT.3.D0)AKSI=1.D0 164 : CCT=AKSI/BM(KP) 165 : TO1(KP)=CCT 166 : TO2(KP)=0.D0 167 : CXT(KP)=UMJ(KP,1)*UMJ(KP,1)*CCT 168 : CYT(KP)=UMJ(KP,2)*UMJ(KP,2)*CCT 169 : CXY(KP)=UMJ(KP,1)*UMJ(KP,2)*CCT 170 : 7018 CONTINUE 171 : ELSEIF(IDCEN.EQ.4)THEN 172 : C Tenseur Visqueux 173 : DT19=DTM1*0.5D0 174 : DO 7009 K=KDEB,KFIN 175 : KP=K-KDEB+1 176 : TO1(KP)=DT19 177 : TO2(KP)=0.D0 178 : CXT(KP)=UMJ(KP,1)*UMJ(KP,1)*DT19 179 : CYT(KP)=UMJ(KP,2)*UMJ(KP,2)*DT19 180 : CXY(KP)=UMJ(KP,1)*UMJ(KP,2)*DT19 181 : 7009 CONTINUE 182 : ENDIF 183 : C 184 : C CALCUL DU PAS DE TEMPS 185 : C 186 : C*IBMDIR* PREFER SCALAR 187 : DO 7010 K=KDEB,KFIN 188 : KP=K-KDEB+1 189 : DT0=DT 190 : DT1=2.D0/ 191 : & (UMJ(KP,1)*UMJ(KP,1)/(CXT(KP)+COEFT(KP)) 192 : & +UMJ(KP,2)*UMJ(KP,2)/(CYT(KP)+COEFT(KP))) 193 : DT2=0.5D0/ 194 : & ((CXT(KP)+COEFT(KP))*AL2(KP) 195 : & +(CYT(KP)+COEFT(KP))*AH2(KP)) 196 : IF(DT1.LT.DT)DT=DT1 197 : IF(DT2.LT.DT)DT=DT2 198 : IF(DT.NE.DT0) THEN 199 : DTT1=DT1 200 : DTT2=DT2 201 : DIAEL=XMB(KP) 202 : NUEL=K 203 : END IF 204 : 7010 CONTINUE 205 : DO 7011 I=1,NP 206 : DO 70111 K=KDEB,KFIN 207 : KP=K-KDEB+1 208 : WTILDE(KP,I)=DRR(I,K) 209 : & +(TO1(KP)*(UMJ(KP,1)*HR(K,I,1)+UMJ(KP,2)*HR(K,I,2)) 210 : & + TO2(KP)*(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2)) 211 : & )*AIRE(KP) 212 : 70111 CONTINUE 213 : 7011 CONTINUE 214 : DO 7012 K=KDEB,KFIN 215 : KP=K-KDEB+1 216 : CHGLD(KP)=COEFT(KP)/12.D0*QUA4*(AL2(KP)+AH2(KP)) 217 : CHGLP(KP)=(CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0*QUA4 218 : 7012 CONTINUE 219 : C *********** 220 : C * 3 D * 221 : C *********** 222 : ELSEIF(IDIM.EQ.3)THEN 223 : CUB8=CHGL 224 : DO 8004 K=KDEB,KFIN 225 : KP=K-KDEB+1 226 : AL2(KP)=1.D0/AL(KP)/AL(KP) 227 : AH2(KP)=1.D0/AH(KP)/AH(KP) 228 : AP2(KP)=1.D0/AP(KP)/AP(KP) 229 : XMB(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0 230 : CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ 231 : & AP(KP)*AH(KP)/AL(KP) 232 : 8004 CONTINUE 233 : DO 8007 K=KDEB,KFIN 234 : KP=K-KDEB+1 235 : UMI(KP,1)=UMI(KP,1)/AIRE(KP) 236 : UMI(KP,2)=UMI(KP,2)/AIRE(KP) 237 : UMI(KP,3)=UMI(KP,3)/AIRE(KP) 238 : UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2) 239 : & +UMI(KP,3)*UMI(KP,3) 240 : UM(KP)=SQRT(UM(KP)) 241 : AX=GRAD(KP,1)*GRAD(KP,1)+GRAD(KP,2)*GRAD(KP,2) 242 : & +GRAD(KP,3)*GRAD(KP,3) 243 : AX=SQRT(AX)+xpeti2 244 : GRAD(KP,1)=GRAD(KP,1)/AX 245 : GRAD(KP,2)=GRAD(KP,2)/AX 246 : GRAD(KP,3)=GRAD(KP,3)/AX 247 : UP(KP)=GRAD(KP,1)*UMI(KP,1)+GRAD(KP,2)*UMI(KP,2) 248 : & +GRAD(KP,3)*UMI(KP,3) 249 : UPI(KP,1)=UP(KP)*GRAD(KP,1) 250 : UPI(KP,2)=UP(KP)*GRAD(KP,2) 251 : UPI(KP,3)=UP(KP)*GRAD(KP,3) 252 : 8007 CONTINUE 253 : DO 80074 K=KDEB,KFIN 254 : KP=K-KDEB+1 255 : BMX=UMI(KP,1)/AL(KP) 256 : BMY=UMI(KP,2)/AH(KP) 257 : BMZ=UMI(KP,3)/AP(KP) 258 : BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ 259 : BM(KP)=SQRT(BM(KP))+xpeti2 260 : BPX=UPI(KP,1)/AL(KP) 261 : BPY=UPI(KP,2)/AH(KP) 262 : BPZ=UPI(KP,3)/AP(KP) 263 : BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ 264 : BP(KP)=SQRT(BP(KP))+xpeti2 265 : 80074 CONTINUE 266 : C 267 : C DECENTREMENT 268 : C 269 : IF(IDCEN.EQ.1)THEN 270 : C CENTREE 271 : DO 80081 K=KDEB,KFIN 272 : KP=K-KDEB+1 273 : TO1(KP)=0.D0 274 : TO2(KP)=0.D0 275 : CXT(KP)=0.D0 276 : CYT(KP)=0.D0 277 : CZT(KP)=0.D0 278 : CXY(KP)=0.D0 279 : CXZ(KP)=0.D0 280 : CYZ(KP)=0.D0 281 : 80081 CONTINUE 282 : ELSEIF(IDCEN.EQ.2)THEN 283 : C SUPGDC 284 : DO 8008 K=KDEB,KFIN 285 : KP=K-KDEB+1 286 : HMK=2.D0*UM(KP)/BM(KP) 287 : XMB(KP)=HMK 288 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 289 : AKSI=ALFA/3.D0 290 : IF(ALFA.GT.3.D0)AKSI=1.D0 291 : CCT=AKSI/BM(KP) 292 : HPK=2.D0*UP(KP)/BP(KP) 293 : ALFA=UP(KP)*HPK/COEFT(KP)/2.D0 294 : AKSI=ALFA/3.D0 295 : IF(ALFA.GT.3.D0)AKSI=1.D0 296 : CCP=AKSI/BP(KP) 297 : CPT=CCP-CCT 298 : CC2=0.D0 299 : IF(CPT.GE.0.D0)CC2=CPT 300 : CC2=CC2*CC0 301 : TO1(KP)=CCT 302 : TO2(KP)=CC2 303 : CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2) 304 : CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2) 305 : CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPI(KP,3)*UPI(KP,3)*CC2) 306 : CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2) 307 : CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPI(KP,1)*UPI(KP,3)*CC2) 308 : CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPI(KP,2)*UPI(KP,3)*CC2) 309 : 8008 CONTINUE 310 : ELSEIF(IDCEN.EQ.3)THEN 311 : C SUPG 312 : DO 8018 K=KDEB,KFIN 313 : KP=K-KDEB+1 314 : HMK=2.D0*UM(KP)/BM(KP) 315 : XMB(KP)=HMK 316 : ALFA=UM(KP)*HMK/COEFT(KP)/2.D0 317 : AKSI=ALFA/3.D0 318 : IF(ALFA.GT.3.D0)AKSI=1.D0 319 : CCT=AKSI/BM(KP) 320 : TO1(KP)=CCT 321 : TO2(KP)=0.D0 322 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT 323 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT 324 : CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT 325 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT 326 : CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT 327 : CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT 328 : 8018 CONTINUE 329 : ELSEIF(IDCEN.EQ.4)THEN 330 : C Tenseur Visqueux 331 : DT19=DTM1*0.5D0 332 : DO 8009 K=KDEB,KFIN 333 : KP=K-KDEB+1 334 : TO1(KP)=DT19 335 : TO2(KP)=0.D0 336 : CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19 337 : CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19 338 : CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19 339 : CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19 340 : CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19 341 : CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19 342 : 8009 CONTINUE 343 : ENDIF 344 : C 345 : C CALCUL DU PAS DE TEMPS 346 : C 347 : C*IBMDIR* PREFER SCALAR 348 : DO 8010 K=KDEB,KFIN 349 : KP=K-KDEB+1 350 : DT0=DT 351 : DT1=2.D0/ 352 : & (UMI(KP,1)*UMI(KP,1)/(CXT(KP)+COEFT(KP)) 353 : & +UMI(KP,2)*UMI(KP,2)/(CYT(KP)+COEFT(KP)) 354 : & +UMI(KP,3)*UMI(KP,3)/(CZT(KP)+COEFT(KP))) 355 : DT2=0.5D0/ 356 : & ((CXT(KP)+COEFT(KP))*AL2(KP) 357 : & +(CYT(KP)+COEFT(KP))*AH2(KP) 358 : & +(CZT(KP)+COEFT(KP))*AP2(KP)) 359 : IF(DT1.LT.DT)DT=DT1 360 : IF(DT2.LT.DT)DT=DT2 361 : IF(DT.NE.DT0) THEN 362 : DTT1=DT1 363 : DTT2=DT2 364 : DIAEL=XMB(KP) 365 : NUEL=K 366 : END IF 367 : 8010 CONTINUE 368 : DO 8011 I=1,NP 369 : DO 80111 K=KDEB,KFIN 370 : KP=K-KDEB+1 371 : WTILDE(KP,I)=DRR(I,K) 372 : & +(TO1(KP)*(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2) 373 : & +UMI(KP,3)*HR(K,I,3)) 374 : & + TO2(KP)*(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2) 375 : & +UPI(KP,3)*HR(K,I,3)) 376 : & )*AIRE(KP) 377 : 80111 CONTINUE 378 : 8011 CONTINUE 379 : DO 8012 K=KDEB,KFIN 380 : KP=K-KDEB+1 381 : GEO1=CFM(KP)*CUB8 382 : CHGLD(KP)=COEFT(KP)*GEO1 383 : CHGLP(KP)=(CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1 384 : 8012 CONTINUE 385 : ENDIF 386 : RETURN 387 : 1002 FORMAT(10(1X,1PE11.4)) 388 : END 389 : 390 : 391 : 392 : 393 : 394 : 395 : 396 : 397 : 398 :
© Cast3M 2003 - All rights reserved.
Disclaimer