1 : C CLI281 SOURCE CHAT 05/01/12 22:07:46 5004 2 : SUBROUTINE CLI281(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO, 3 : & ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYN,ICHLIM,ICHRES,ICHRLI) 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : CLI281 9 : C 10 : C DESCRIPTION : Subroutine appellée par CLIM22 11 : C calcul de RESIDU et CLIM at the board 12 : C OPTION: 'INOU' 2D 13 : C 14 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 15 : C 16 : C AUTEUR : S.Kudriakov, DEN/DM2S/SFME/LTMF 17 : C 18 : C************************************************************************ 19 : C 20 : C APPELES (Calcul) : 21 : C 22 : C************************************************************************ 23 : C 24 : C HISTORIQUE (Anomalies et modifications éventuelles) 25 : C 26 : C HISTORIQUE : 27 : C 28 : C************************************************************************ 29 : C 30 : IMPLICIT INTEGER(I-N) 31 : -INC CCOPTIO 32 : -INC SMLMOTS 33 : -INC SMELEME 34 : POINTEUR MELEFC.MELEME 35 : -INC SMLENTI 36 : POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI 37 : -INC SMCHPOI 38 : POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL, 39 : & MPVC.MPOVAL, MPPC.MPOVAL, MPYN.MPOVAL, MPLIM.MPOVAL, 40 : & MPRES.MPOVAL, MPRLI.MPOVAL 41 : C---------------------------------------- 42 : C**** Variables de COOPTIO 43 : C---------------------------------------- 44 : INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 45 : & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 46 : & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 47 : & ,IECHO, IIMPI, IOSPI 48 : & ,IDIM, IFICLE, IPREFI 49 : & ,MCOORD 50 : & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 51 : & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 52 : & ,NORINC,NORVAL,NORIND,NORVAD 53 : & ,NUCROU, IPSAUV 54 : C 55 : INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC 56 : & ,IYN,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC 57 : & ,NGF,NGC,NLF,NLC,NLCB,LRECP,LRECV,I,NSP,NESP 58 : REAL*8 VOLU,SURF,RC,PC,UXC,UYC,UZC,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ 59 : & ,CT2X,CT2Y,CT2Z,RF,PF,UXF,UYF,UZF,TOP,BOT 60 : & ,UNC,UNF,UTF,UT2F,SF,GAMF,ECIN,PSRF,HTF,GM1 61 : & ,CELLT,UT2C,UTC,TEMP 62 : CHARACTER*(8) TYPE 63 : C------------------------------------------------------------ 64 : -INC SMLREEL 65 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL 66 : C------------------------------------------------------- 67 : C********** Les CP's and CV's *********************** 68 : C------------------------------------------------------- 69 : SEGMENT GCONST 70 : REAL*8 GC(NSP) 71 : ENDSEGMENT 72 : POINTEUR CP.GCONST, CV.GCONST 73 : C------------------------------------------------------------- 74 : C******* Les fractionines massiques ************************** 75 : C------------------------------------------------------------- 76 : SEGMENT FRAMAS 77 : REAL*8 YET(NSP) 78 : ENDSEGMENT 79 : POINTEUR YC.FRAMAS, YF.FRAMAS 80 : C------------------------------------------------------------- 81 : C********** Segments for the flux-vector ******************* 82 : C------------------------------------------------------------- 83 : SEGMENT FUNEL 84 : REAL*8 FU(4+NSP) 85 : ENDSEGMENT 86 : POINTEUR flux2D.funel, flux3D.funel 87 : SEGINI FLUX2D 88 : SEGINI FLUX3D 89 : C------------------------------------------------------ 90 : C**** KRIPAD pour la correspondance global/local 91 : C------------------------------------------------------ 91 : CALL KRIPAD(MELEMC,MLEMC) 91 : CALL KRIPAD(MELECB,MLEMCB) 91 : CALL KRIPAD(MELEMF,MLEMF) 92 : C------------------------------------------------------ 93 : C**** CHPOINTs de la table DOMAINE 94 : C------------------------------------------------------ 94 : CALL LICHT(INORM,MPNORM,TYPE,ICEL) 94 : CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL) 94 : CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL) 95 : C------------------------------------------------------ 96 : C**** CHPOINTs des variables 97 : C------------------------------------------------------ 97 : CALL LICHT(IROC,MPRC,TYPE,ICEL) 97 : CALL LICHT(IVITC,MPVC,TYPE,ICEL) 97 : CALL LICHT(IPC,MPPC,TYPE,ICEL) 97 : CALL LICHT(IYN,MPYN,TYPE,ICEL) 97 : CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL) 97 : CALL LICHT(ICHRES,MPRES,TYPE,ICEL) 97 : CALL LICHT(ICHRLI,MPRLI,TYPE,ICEL) 98 : C--------------------------------------------------------- 99 : C**** Boucle sur le face pour le calcul des invariants de 100 : C Riemann et du flux 101 : C--------------------------------------------------------- 102 : SEGACT MELEFC 103 : NFAC=MELEFC.NUM(/2) 104 : c write(*,*) 'hell', NFAC 105 : UZC=0.0D0 106 : CNZ=0.0D0 107 : CTZ=0.0D0 108 : CT2X=0.0D0 109 : CT2Y=0.0D0 110 : CT2Z=0.0D0 111 : DO 1 IFAC=1,NFAC,1 112 : NGF=MELEFC.NUM(1,IFAC) 113 : NGC=MELEFC.NUM(2,IFAC) 114 : NLF=MLEMF.LECT(NGF) 115 : NLC=MLEMC.LECT(NGC) 116 : NLCB=MLEMCB.LECT(NGF) 117 : VOLU=MPVOL.VPOCHA(NLC,1) 118 : SURF=MPSURF.VPOCHA(NLF,1) 119 : C---------------------------------------------- 120 : C In CASTEM les normales sont sortantes 121 : C---------------------------------------------- 122 : CNX=MPNORM.VPOCHA(NLF,1) 123 : CNY=MPNORM.VPOCHA(NLF,2) 124 : IF(IDIM.EQ.2)THEN 125 : CTX=-1.0D0*CNY 126 : CTY=CNX 127 : ELSE 128 : CNZ=MPNORM.VPOCHA(NLF,3) 129 : CTX=MPNORM.VPOCHA(NLF,4) 130 : CTY=MPNORM.VPOCHA(NLF,5) 131 : CTZ=MPNORM.VPOCHA(NLF,6) 132 : CT2X=MPNORM.VPOCHA(NLF,7) 133 : CT2Y=MPNORM.VPOCHA(NLF,8) 134 : CT2Z=MPNORM.VPOCHA(NLF,9) 135 : ENDIF 136 : C---------------------------------------- 137 : SEGINI CP, CV 138 : MLRECP = LRECP 139 : MLRECV = LRECV 140 : SEGACT MLRECP, MLRECV 141 : DO 10 I=1,(NSP-1) 142 : CP.GC(I)=MLRECP.PROG(I) 143 : CV.GC(I)=MLRECV.PROG(I) 144 : 10 CONTINUE 145 : CP.GC(NSP)=MLRECP.PROG(NSP) 146 : CV.GC(NSP)=MLRECV.PROG(NSP) 147 : C---------------------------- 148 : C Variables au centre 149 : C---------------------------- 150 : RC=MPRC.VPOCHA(NLC,1) 151 : PC=MPPC.VPOCHA(NLC,1) 152 : UXC=MPVC.VPOCHA(NLC,1) 153 : UYC=MPVC.VPOCHA(NLC,2) 154 : IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3) 155 : SEGINI YC 156 : SEGACT MPYN 157 : DO 100 I=1,(NSP-1) 158 : YC.YET(I)=MPYN.VPOCHA(NLC,I) 159 : 100 CONTINUE 160 : c------------------------------------------------------------- 161 : c Computing GAMMA at the cell-center 162 : c------------------------------------------------------------- 163 : top=0.0D0 164 : bot=0.0D0 165 : do 102 i=1,(nsp-1) 166 : top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp)) 167 : bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp)) 168 : 102 continue 169 : top=cp.gc(nsp)+top 170 : bot=cv.gc(nsp)+bot 171 : GAMC=top/bot 172 : C------------------------------------------------------------- 173 : C We decide whether the boundary conditions are 'inlet' or 174 : C 'outlet' by looking at the UNC (UNC > 0.0 -- 'outlet') 175 : C------------------------------------------------------------- 176 : UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ) 177 : c write(*,*) 'before', UNC, NLC,NGC 178 : IF(UNC .GE. 0.0D0) THEN 179 : C-------------------------------------- 180 : C 'OUTP' -- outlet with given pressure 181 : C-------------------------------------- 182 : C----------------------------------------- 183 : C Variables à la face 184 : C----------------------------------------- 185 : PF=MPLIM.VPOCHA(NLCB,NSP+2) 186 : C--------------------------------------- 187 : C******* On calcule UN, UT, UT2 188 : C--------------------------------------- 189 : UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ) 190 : UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z) 191 : C----------------------------------------------- 192 : C******* Densite, vitesse, pression sur le bord 193 : C----------------------------------------------- 194 : MPRLI.VPOCHA(NLCB,1)=RC 195 : MPRLI.VPOCHA(NLCB,2)=UXC 196 : MPRLI.VPOCHA(NLCB,3)=UYC 197 : IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZC 198 : MPRLI.VPOCHA(NLCB,IDIM+2)=PF 199 : do 104 i=1,(nsp-1) 200 : MPRLI.VPOCHA(NLCB,IDIM+2+I)=YC.YET(I) 201 : 104 continue 202 : C--------------------------------------------------- 203 : C******* Probleme de Riemann entre l'etat gauche 204 : C RC,UNC,UTC,UT2C,PC et l'etat droite 205 : C RC,UNC,UTC,UT2C,PF 206 : C On utilise AUSM+ 207 : C Flux dans le repaire normale 208 : C--------------------------------------------------- 209 : NESP=NSP-1 210 : IF(IDIM.EQ.2)THEN 210 : CALL FAUSMP(NESP, 211 : & GAMC,RC,PC,UNC,UTC, 212 : & GAMC,RC,PF,UNC,UTC, 213 : & YC.YET,YC.YET, 214 : & FLUX2D.FU, 215 : & CELLT) 216 : C------------------------------------------------------- 217 : C******* Residuum (son SPG a le meme ordre que MELEFC) 218 : C------------------------------------------------------- 219 : MPRES.VPOCHA(IFAC,1)=-1*FLUX2D.FU(1)*SURF/VOLU 220 : MPRES.VPOCHA(IFAC,2)=-1*((FLUX2D.FU(2)*CNX)+ 221 : & (FLUX2D.FU(3)*CTX))*SURF/VOLU 222 : MPRES.VPOCHA(IFAC,3)=-1*((FLUX2D.FU(2)*CNY)+ 223 : & (FLUX2D.FU(3)*CTY))*SURF/VOLU 224 : MPRES.VPOCHA(IFAC,4)=-1*FLUX2D.FU(4)*SURF/VOLU 225 : do 105 i=1,(nsp-1) 226 : MPRES.VPOCHA(IFAC,4+I)=-1*FLUX2D.FU(4+I)*SURF/VOLU 227 : 105 continue 228 : ELSE 228 : CALL FAUSM3(NESP, 229 : & GAMC,RC,PC,UNC,UTC,UT2C, 230 : & GAMC,RC,PF,UNC,UTC,UT2C, 231 : & YC.YET,YC.YET, 232 : & FLUX3D.FU, 233 : & CELLT) 234 : C------------------------------------------------------ 235 : C******* Residuum (son SPG a le meme ordre que MELEFC) 236 : C------------------------------------------------------ 237 : MPRES.VPOCHA(IFAC,1)=-1*FLUX3D.FU(1)*SURF/VOLU 238 : MPRES.VPOCHA(IFAC,2)=-1*((FLUX3D.FU(2)*CNX)+ 239 : & (FLUX3D.FU(3)*CTX)+(FLUX3D.FU(4)*CT2X))*SURF/VOLU 240 : MPRES.VPOCHA(IFAC,3)=-1*((FLUX3D.FU(2)*CNY)+ 241 : & (FLUX3D.FU(3)*CTY)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU 242 : MPRES.VPOCHA(IFAC,4)=-1*((FLUX3D.FU(2)*CNZ)+ 243 : & (FLUX3D.FU(3)*CTZ)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU 244 : MPRES.VPOCHA(IFAC,5)=-1*FLUX3D.FU(5)*SURF/VOLU 245 : do 106 i=1,(nsp-1) 246 : MPRES.VPOCHA(IFAC,5+I)=-1*FLUX3D.FU(5+I)*SURF/VOLU 247 : 106 continue 248 : ENDIF 249 : TEMP=-1000.0 250 : ELSE 251 : c write(*,*) 'insuuuuuuuu', NLCB 252 : C------------------------------------------------------- 253 : C 'INSU' -- 'inlet' b.c. with given Ht and S (entropy) 254 : C------------------------------------------------------- 255 : UZF=0.0D0 256 : UT2F=0.0D0 257 : C---------------------------- 258 : C Variables à la face 259 : C---------------------------- 260 : HTF=MPLIM.VPOCHA(NLCB,1) 261 : SF=MPLIM.VPOCHA(NLCB,2) 262 : SEGINI YF 263 : TEMP=1000.0 264 : DO 101 I=1,(NSP-1) 265 : YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I) 266 : 101 CONTINUE 267 : UTF=0.0D0 268 : c------------------------------------------------------------- 269 : c Computing GAMMA at the face-center 270 : c------------------------------------------------------------- 271 : top=0.0D0 272 : bot=0.0D0 273 : do 103 i=1,(nsp-1) 274 : top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp)) 275 : bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp)) 276 : 103 continue 277 : top=cp.gc(nsp)+top 278 : bot=cv.gc(nsp)+bot 279 : GAMF=top/bot 280 : GM1=GAMF-1.0D0 281 : C----------------------------------------- 282 : C changing direction of the normal vector 283 : C----------------------------------------- 284 : CNX=-1.0D0*CNX 285 : CNY=-1.0D0*CNY 286 : IF(IDIM.EQ.2)THEN 287 : CTX=-1.0D0*CNY 288 : CTY=CNX 289 : ELSE 290 : CNZ=-1.0D0*CNZ 291 : CTX=-1.0D0*CTX 292 : CTY=-1.0D0*CTY 293 : CTZ=-1.0D0*CTZ 294 : CT2X=-1.0D0*CT2X 295 : CT2Y=-1.0D0*CT2Y 296 : CT2Z=-1.0D0*CT2Z 297 : ENDIF 298 : C--------------------------------------- 299 : C******* On calcule UN, UT, UT2, ASON, S 300 : C--------------------------------------- 301 : UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ) 302 : UNF=UNC 303 : UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ) 304 : C---------------------------------- 305 : UXF=UNF*CNX+UTF*CTX+UT2F*CT2X 306 : UYF=UNF*CNY+UTF*CTY+UT2F*CT2Y 307 : UZF=UNF*CNZ+UTF*CTZ+UT2F*CT2Z 308 : C---------------------------------- 309 : ECIN=0.5D0*((UXF*UXF)+(UYF*UYF)+(UZF*UZF)) 310 : PSRF=(GM1/GAMF)*(HTF-ECIN) 311 : RF=PSRF/SF 312 : RF=RF**(1.0D0/GM1) 313 : PF=SF*(RF**GAMF) 314 : C----------------------------------------------- 315 : C******* Densite, vitesse, pression sur le bord 316 : C----------------------------------------------- 317 : MPRLI.VPOCHA(NLCB,1)=RF 318 : MPRLI.VPOCHA(NLCB,2)=UXF 319 : MPRLI.VPOCHA(NLCB,3)=UYF 320 : IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZF 321 : MPRLI.VPOCHA(NLCB,IDIM+2)=PF 322 : do 1040 i=1,(nsp-1) 323 : MPRLI.VPOCHA(NLCB,IDIM+2+I)=YF.YET(I) 324 : 1040 continue 325 : C--------------------------------------------------- 326 : C******* Probleme de Riemann entre l'etat gauche 327 : C RF,UNC,UTF,UT2F,PF et l'etat droite 328 : C RC,UNC,UTC,UT2C,PC 329 : C On utilise AUSM+ 330 : C Flux dans le repaire normale 331 : C--------------------------------------------------- 332 : NESP=NSP-1 333 : IF(IDIM.EQ.2)THEN 333 : CALL FAUSMP(NESP, 334 : & GAMF,RF,PF,UNC,UTF, 335 : & GAMC,RC,PC,UNC,UTC, 336 : & YF.YET,YC.YET, 337 : & FLUX2D.FU, 338 : & CELLT) 339 : C------------------------------------------------------- 340 : C******* Residuum (son SPG a le meme ordre que MELEFC) 341 : C------------------------------------------------------- 342 : MPRES.VPOCHA(IFAC,1)=FLUX2D.FU(1)*SURF/VOLU 343 : MPRES.VPOCHA(IFAC,2)=((FLUX2D.FU(2)*CNX)+(FLUX2D.FU(3)*CTX)) 344 : & *SURF/VOLU 345 : MPRES.VPOCHA(IFAC,3)=((FLUX2D.FU(2)*CNY)+(FLUX2D.FU(3)*CTY)) 346 : & *SURF/VOLU 347 : MPRES.VPOCHA(IFAC,4)=FLUX2D.FU(4)*SURF/VOLU 348 : do 1050 i=1,(nsp-1) 349 : MPRES.VPOCHA(IFAC,4+I)=FLUX2D.FU(4+I)*SURF/VOLU 350 : 1050 continue 351 : ELSE 351 : CALL FAUSM3(NESP, 352 : & GAMF,RF,PF,UNC,UTF,UT2F, 353 : & GAMC,RC,PC,UNC,UTC,UT2C, 354 : & YF.YET,YC.YET, 355 : & FLUX3D.FU, 356 : & CELLT) 357 : C------------------------------------------------------ 358 : C******* Residuum (son SPG a le meme ordre que MELEFC) 359 : C------------------------------------------------------ 360 : MPRES.VPOCHA(IFAC,1)=FLUX3D.FU(1)*SURF/VOLU 361 : MPRES.VPOCHA(IFAC,2)=((FLUX3D.FU(2)*CNX)+(FLUX3D.FU(3)*CTX)+ 362 : & (FLUX3D.FU(4)*CT2X))*SURF/VOLU 363 : MPRES.VPOCHA(IFAC,3)=((FLUX3D.FU(2)*CNY)+(FLUX3D.FU(3)*CTY)+ 364 : & (FLUX3D.FU(4)*CT2Z))*SURF/VOLU 365 : MPRES.VPOCHA(IFAC,4)=((FLUX3D.FU(2)*CNZ)+(FLUX3D.FU(3)*CTZ)+ 366 : & (FLUX3D.FU(4)*CT2Z))*SURF/VOLU 367 : MPRES.VPOCHA(IFAC,5)=FLUX3D.FU(5)*SURF/VOLU 368 : do 1060 i=1,(nsp-1) 369 : MPRES.VPOCHA(IFAC,5+I)=FLUX3D.FU(5+I)*SURF/VOLU 370 : 1060 continue 371 : ENDIF 372 : ENDIF 373 : 1 CONTINUE 374 : C 375 : SEGDES MELEFC 376 : C 377 : SEGDES MLEMC 378 : SEGDES MLEMCB 379 : SEGDES MLEMF 380 : C 381 : SEGDES MPNORM 382 : SEGDES MPVOL 383 : SEGDES MPSURF 384 : SEGDES MPRC 385 : SEGDES MPPC 386 : SEGDES MPVC 387 : SEGDES MPYN 388 : SEGDES MPLIM 389 : SEGDES MPRES 390 : SEGDES MPRLI 391 : SEGDES MLRECP 392 : SEGDES MLRECV 393 : SEGDES YC 394 : IF(TEMP .GT. 0.0D0) SEGDES YF 395 : SEGDES FLUX2D 396 : SEGDES FLUX3D 397 : C 398 : 9999 CONTINUE 399 : RETURN 400 : END 401 : 402 : 403 : 404 : 405 : 406 : 407 : 408 : 409 :
© Cast3M 2003 - All rights reserved.
Disclaimer