1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23!> \file cfxtcl.f90 24!> \brief Handle boundary condition type code (\ref itypfb) when the 25!> compressible model is enabled. 26!> 27!> Please refer to the 28!> <a href="../../theory.pdf#cfxtcl"><b>cfxtcl</b></a> 29!> section of the theory guide for more informations. 30!------------------------------------------------------------------------------- 31 32!------------------------------------------------------------------------------ 33! Arguments 34!------------------------------------------------------------------------------ 35! mode name role 36!------------------------------------------------------------------------------ 37!> \param[in] nvar total number of variables 38!> \param[in,out] icodcl face boundary condition code: 39!> - 1 Dirichlet 40!> - 2 Radiative outlet 41!> - 3 Neumann 42!> - 4 sliding and 43!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 44!> - 5 smooth wall and 45!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 46!> - 6 rough wall and 47!> \f$ \vect{u} \cdot \vect{n} = 0 \f$ 48!> - 9 free inlet/outlet 49!> (input mass flux blocked to 0) 50!> - 13 Dirichlet for the advection operator and 51!> Neumann for the diffusion operator 52!> \param[in] itypfb boundary face types 53!> \param[in] dt time step (per cell) 54!> \param[in,out] rcodcl boundary condition values: 55!> - rcodcl(1) value of the dirichlet 56!> - rcodcl(2) value of the exterior exchange 57!> coefficient (infinite if no exchange) 58!> - rcodcl(3) value flux density 59!> (negative if gain) in w/m2 or roughness 60!> in m if icodcl=6 61!> -# for the velocity \f$ (\mu+\mu_T) 62!> \gradv \vect{u} \cdot \vect{n} \f$ 63!> -# for the pressure \f$ \Delta t 64!> \grad P \cdot \vect{n} \f$ 65!> -# for a scalar \f$ cp \left( K + 66!> \dfrac{K_T}{\sigma_T} \right) 67!> \grad T \cdot \vect{n} \f$ 68!______________________________________________________________________________ 69 70subroutine cfxtcl & 71 ( nvar , & 72 icodcl , itypfb , dt , rcodcl ) 73 74!=============================================================================== 75! Module files 76!=============================================================================== 77 78use paramx 79use numvar 80use optcal 81use cstphy 82use cstnum 83use entsor 84use parall 85use ppppar 86use ppthch 87use ppincl 88use cfpoin 89use mesh 90use field 91use cs_cf_bindings 92 93!=============================================================================== 94 95implicit none 96 97! Arguments 98 99integer nvar 100 101integer icodcl(nfabor,nvar) 102integer itypfb(nfabor) 103 104double precision dt(ncelet) 105double precision rcodcl(nfabor,nvar,3) 106 107! Local variables 108 109integer ifac , iel 110integer ii , iccfth 111integer icalep 112integer ien , itk, niv 113integer nvarcf 114 115integer nvcfmx 116parameter (nvcfmx=6) 117integer ivarcf(nvcfmx) 118 119double precision hint, bmasfl, drom 120 121double precision, allocatable, dimension(:) :: w5 122double precision, allocatable, dimension(:) :: w7 123double precision, allocatable, dimension(:) :: wbfb, wbfa 124double precision, allocatable, dimension(:) :: bc_en, bc_pr, bc_tk 125double precision, allocatable, dimension(:) :: bc_fracv, bc_fracm, bc_frace 126double precision, allocatable, dimension(:,:) :: bc_vel 127 128double precision, dimension(:), pointer :: coefbp 129double precision, dimension(:), pointer :: crom, brom, cpro_cv, cvar_en 130double precision, dimension(:,:), pointer :: vel 131double precision, dimension(:), pointer :: cvar_fracv, cvar_fracm, cvar_frace 132 133!=============================================================================== 134 135! Map field arrays 136 call field_get_val_v(ivarfl(iu), vel) 137 138!=============================================================================== 139! 1. Initializations 140!=============================================================================== 141 142! Allocate temporary arrays 143allocate(w5(ncelet)) 144allocate(w7(nfabor), wbfa(nfabor), wbfb(nfabor)) 145allocate(bc_en(nfabor)) 146allocate(bc_pr(nfabor)) 147allocate(bc_tk(nfabor)) 148allocate(bc_fracv(nfabor)) 149allocate(bc_fracm(nfabor)) 150allocate(bc_frace(nfabor)) 151allocate(bc_vel(3,nfabor)) 152 153ien = isca(ienerg) 154itk = isca(itempk) 155 156call field_get_val_s(icrom, crom) 157call field_get_val_s(ibrom, brom) 158 159call field_get_val_s(ivarfl(ien), cvar_en) 160 161if (icv.ge.0) call field_get_val_s(icv, cpro_cv) 162 163! mixture fractions for the homogeneous two-phase flows 164if (ippmod(icompf).eq.2) then 165 call field_get_val_s(ivarfl(isca(ifracv)), cvar_fracv) 166 call field_get_val_s(ivarfl(isca(ifracm)), cvar_fracm) 167 call field_get_val_s(ivarfl(isca(ifrace)), cvar_frace) 168endif 169 170! list of the variables of the compressible model 171ivarcf(1) = ipr 172ivarcf(2) = iu 173ivarcf(3) = iv 174ivarcf(4) = iw 175ivarcf(5) = ien 176ivarcf(6) = itk 177nvarcf = 6 178 179call field_get_coefb_s(ivarfl(ipr), coefbp) 180do ifac = 1, nfabor 181 wbfb(ifac) = coefbp(ifac) 182enddo 183 184! Computation of epsilon_sup = e - CvT 185! Needed if walls with imposed temperature are set. 186 187icalep = 0 188do ifac = 1, nfabor 189 if(icodcl(ifac,itk).eq.5) then 190 icalep = 1 191 endif 192enddo 193if(icalep.ne.0) then 194 ! At cell centers 195 call cs_cf_thermo_eps_sup(crom, w5, ncel) 196 197 ! At boundary faces centers 198 call cs_cf_thermo_eps_sup(brom, w7, nfabor) 199endif 200 201! Loop on all boundary faces and treatment of types of BCs given by itypfb 202 203do ifac = 1, nfabor 204 iel = ifabor(ifac) 205 206!=============================================================================== 207! 2. Treatment of all wall boundary faces 208!=============================================================================== 209 210 if (itypfb(ifac).eq.iparoi .or. itypfb(ifac).eq.iparug) then 211 212 ! pressure : 213 ! if the gravity is prevailing: hydrostatic pressure 214 ! (warning: the density is here explicit and the term is an approximation) 215 216 if(icfgrp.eq.1) then 217 218 icodcl(ifac,ipr) = 15 219 hint = dt(iel)/distb(ifac) 220 rcodcl(ifac,ipr,3) = -hint & 221 * ( gx*(cdgfbo(1,ifac)-xyzcen(1,iel)) & 222 + gy*(cdgfbo(2,ifac)-xyzcen(2,iel)) & 223 + gz*(cdgfbo(3,ifac)-xyzcen(3,iel)) ) & 224 * crom(iel) 225 226 else 227 228 ! generally proportional to the bulk value 229 ! (Pboundary = COEFB*Pi) 230 ! The part deriving from pinf in stiffened gas is set in explicit for now 231 call cs_cf_thermo_wall_bc(wbfa, wbfb, ifac-1) 232 233 if (wbfb(ifac).lt.rinfin*0.5d0.and.wbfb(ifac).gt.0.d0) then 234 icodcl(ifac,ipr) = 12 235 rcodcl(ifac,ipr,1) = wbfa(ifac) 236 rcodcl(ifac,ipr,2) = wbfb(ifac) 237 else 238 ! If rarefaction is too strong : homogeneous Dirichlet 239 icodcl(ifac,ipr) = 13 240 rcodcl(ifac,ipr,1) = 0.d0 241 endif 242 243 endif 244 245 ! Velocity and turbulence are treated in a standard manner in condli. 246 247 ! For thermal B.C., a pre-treatment has be done here since the solved 248 ! variable is the total energy 249 ! (internal energy + epsilon_sup + cinetic energy). 250 ! Especially, when a temperature is imposed on a wall, clptur treatment 251 ! has to be prepared. Except for the solved energy all the variables rho 252 ! and s will take arbitrarily a zero flux B.C. (their B.C. are only used 253 ! for the gradient reconstruction and imposing something else than zero 254 ! flux could bring out spurious values near the boundary layer). 255 256 ! adiabatic by default 257 if( icodcl(ifac,itk).eq.0.and. & 258 icodcl(ifac,ien).eq.0) then 259 icodcl(ifac,itk) = 3 260 rcodcl(ifac,itk,3) = 0.d0 261 endif 262 263 ! imposed temperature 264 if(icodcl(ifac,itk).eq.5) then 265 266 ! The value of the energy that leads to the right flux is imposed. 267 ! However it should be noted that it is the B.C. for the diffusion 268 ! flux. For the gradient reconstruction, something else will be 269 ! needed. For example, a zero flux or an other B.C. respecting a 270 ! profile: it may be possible to treat the total energy as the 271 ! temperature, keeping in mind that the total energy contains 272 ! the cinetic energy, which could make the choice of the profile more 273 ! difficult. 274 275 icodcl(ifac,ien) = 5 276 if(icv.eq.-1) then 277 rcodcl(ifac,ien,1) = cv0*rcodcl(ifac,itk,1) 278 else 279 rcodcl(ifac,ien,1) = cpro_cv(iel)*rcodcl(ifac,itk,1) 280 endif 281 rcodcl(ifac,ien,1) = rcodcl(ifac,ien,1) & 282 + 0.5d0*(vel(1,iel)**2+vel(2,iel)**2+vel(3,iel)**2) & 283 + w5(iel) 284 ! w5 contains epsilon_sup 285 286 ! fluxes in grad(epsilon_sup and cinetic energy) have to be zero 287 ! since they are already accounted for in the energy diffusion term 288 ifbet(ifac) = 1 289 290 ! Dirichlet condition on the temperature for gradient reconstruction 291 ! used only in post-processing (typically Nusselt computation) 292 icodcl(ifac,itk) = 1 293 294 ! imposed flux 295 elseif(icodcl(ifac,itk).eq.3) then 296 297 ! zero flux on energy 298 icodcl(ifac,ien) = 3 299 rcodcl(ifac,ien,3) = rcodcl(ifac,itk,3) 300 301 ! fluxes in grad(epsilon_sup and cinetic energy) have to be zero 302 ! since they are already accounted for in the energy diffusion term 303 ifbet(ifac) = 1 304 305 ! zero flux for the possible temperature reconstruction 306 icodcl(ifac,itk) = 3 307 rcodcl(ifac,itk,3) = 0.d0 308 309 endif 310 311!=============================================================================== 312! 3. Treatment of all inlet/outlet boundary faces and thermo step 313!=============================================================================== 314 315!=============================================================================== 316! 3.1 Imposed Inlet/outlet (for example: supersonic inlet) 317!=============================================================================== 318 319 elseif ( itypfb(ifac).eq.iesicf ) then 320 321 ! we have 322 ! - velocity, 323 ! - 2 variables among P, rho, T, E (but not the couple (T,E)), 324 ! - turbulence variables 325 ! - scalars 326 327 ! we look for the variable to be initialized 328 ! (if a zero value has been given, it is not adapted, so it will 329 ! be considered as not initialized and the computation will stop 330 ! displaying an error message. The boundary density may 331 ! be pre-initialized to the cell density also, so is tested last. 332 333 iel = ifabor(ifac) 334 drom = abs(crom(iel) - brom(ifac)) 335 336 niv = 0 337 iccfth = 10000 338 if (rcodcl(ifac,ipr,1).lt.rinfin*0.5d0) then 339 iccfth = 2*iccfth 340 niv = niv + 1 341 endif 342 if (rcodcl(ifac,itk,1).lt.rinfin*0.5d0) then 343 iccfth = 5*iccfth 344 niv = niv + 1 345 endif 346 if (rcodcl(ifac,ien,1).lt.rinfin*0.5d0) then 347 iccfth = 7*iccfth 348 niv = niv + 1 349 endif 350 351 if (brom(ifac).gt.0.d0 .and. (niv.lt.2 .or. drom.gt.epzero)) then 352 iccfth = 3*iccfth 353 niv = niv + 1 354 endif 355 356 if (niv .ne. 2) then 357 write(nfecra,1000) iccfth 358 call csexit (1) 359 endif 360 iccfth = iccfth + 900 361 362 ! missing thermo variables among P, rho, T, E are computed 363 bc_en(ifac) = rcodcl(ifac,ien,1) 364 bc_pr(ifac) = rcodcl(ifac,ipr,1) 365 bc_tk(ifac) = rcodcl(ifac,itk,1) 366 bc_vel(1,ifac) = rcodcl(ifac,iu,1) 367 bc_vel(2,ifac) = rcodcl(ifac,iv,1) 368 bc_vel(3,ifac) = rcodcl(ifac,iw,1) 369 370 call cs_cf_thermo(iccfth, ifac-1, bc_en, bc_pr, bc_tk, bc_vel) 371 372!=============================================================================== 373! 3.2 Outlet with imposed pressure 374!=============================================================================== 375 376 elseif ( itypfb(ifac).eq.isopcf ) then 377 378 ! If no value was given for P or if its value is negative, the computation 379 ! stops (a negative value could be possible, but in most cases it would be 380 ! an error). 381 if(rcodcl(ifac,ipr,1).lt.-rinfin*0.5d0) then 382 write(nfecra,1100) 383 call csexit (1) 384 endif 385 386 bc_en(ifac) = rcodcl(ifac,ien,1) 387 bc_pr(ifac) = rcodcl(ifac,ipr,1) 388 bc_tk(ifac) = rcodcl(ifac,itk,1) 389 bc_vel(1,ifac) = rcodcl(ifac,iu,1) 390 bc_vel(2,ifac) = rcodcl(ifac,iv,1) 391 bc_vel(3,ifac) = rcodcl(ifac,iw,1) 392 393 call cs_cf_thermo_subsonic_outlet_bc(bc_en, bc_pr, bc_vel, ifac-1) 394 395!=============================================================================== 396! 3.3 Inlet with Ptot, Htot imposed (reservoir boundary conditions) 397!=============================================================================== 398 399 elseif ( itypfb(ifac).eq.iephcf ) then 400 401 ! If values for Ptot and Htot were not given, the computation stops. 402 403 ! rcodcl(ifac,isca(ienerg),1) contains the boundary total enthalpy values 404 ! prescribed by the user 405 406 if(rcodcl(ifac,ipr ,1).lt.-rinfin*0.5d0.or. & 407 rcodcl(ifac,isca(ienerg) ,1).lt.-rinfin*0.5d0) then 408 write(nfecra,1200) 409 call csexit (1) 410 endif 411 412 bc_en(ifac) = rcodcl(ifac,ien,1) 413 bc_pr(ifac) = rcodcl(ifac,ipr,1) 414 bc_tk(ifac) = rcodcl(ifac,itk,1) 415 bc_vel(1,ifac) = rcodcl(ifac,iu,1) 416 bc_vel(2,ifac) = rcodcl(ifac,iv,1) 417 bc_vel(3,ifac) = rcodcl(ifac,iw,1) 418 419 call cs_cf_thermo_ph_inlet_bc(bc_en, bc_pr, bc_vel, ifac-1) 420 421!=============================================================================== 422! 3.4 Inlet with imposed rho*U and rho*U*H 423!=============================================================================== 424 425 elseif ( itypfb(ifac).eq.ieqhcf ) then 426 427 ! TODO to be implemented 428 write(nfecra,1301) 429 call csexit (1) 430 431 ! On utilise un scenario dans lequel on a un 2-contact et une 432 ! 3-détente entrant dans le domaine. On détermine les conditions 433 ! sur l'interface selon la thermo et on passe dans Rusanov 434 ! ensuite pour lisser. 435 436 ! Si rho et u ne sont pas donnés, erreur 437 if(rcodcl(ifac,irunh,1).lt.-rinfin*0.5d0) then 438 write(nfecra,1300) 439 call csexit (1) 440 endif 441 442 endif ! end of test on boundary condition types 443 444!=============================================================================== 445! 4. Complete the treatment for inlets and outlets: 446! - boundary convective fluxes computation (analytical or Rusanov) if needed 447! - B.C. code (Dirichlet or Neumann) 448! - Dirichlet values 449!=============================================================================== 450 451 if (itypfb(ifac).eq.iesicf.or. & 452 itypfb(ifac).eq.isspcf.or. & 453 itypfb(ifac).eq.iephcf.or. & 454 itypfb(ifac).eq.isopcf.or. & 455 itypfb(ifac).eq.ieqhcf) then 456 457!=============================================================================== 458! 4.1 Boundary convective fluxes computation (analytical or Rusanov) if needed 459! (gamma should already have been computed if Rusanov fluxes are computed) 460!=============================================================================== 461 462 ! Rusanov fluxes are computed only for the imposed inlet for stability 463 ! reasons. 464 if (itypfb(ifac).eq.iesicf) then 465 466 ! Dirichlet for velocity and pressure are computed in order to 467 ! impose the Rusanov fluxes in mass, momentum and energy balance. 468 call cfrusb(ifac, bc_en, bc_pr, bc_vel) 469 470 ! For the other types of inlets/outlets (subsonic outlet, QH inlet, 471 ! PH inlet), analytical fluxes are computed 472 elseif (itypfb(ifac).ne.isspcf) then 473 474 ! the pressure part of the boundary analytical flux is not added here, 475 ! but set through the pressure gradient boundary conditions (Dirichlet) 476 call cffana(ifac, bc_en, bc_pr, bc_vel) 477 478 endif 479 480!=============================================================================== 481! 4.2 Copy of boundary values into the Dirichlet values array 482!=============================================================================== 483 484 if (itypfb(ifac).ne.isspcf) then 485 rcodcl(ifac,ien,1) = bc_en(ifac) 486 rcodcl(ifac,ipr,1) = bc_pr(ifac) 487 rcodcl(ifac,itk,1) = bc_tk(ifac) 488 rcodcl(ifac,iu,1) = bc_vel(1,ifac) 489 rcodcl(ifac,iv,1) = bc_vel(2,ifac) 490 rcodcl(ifac,iw,1) = bc_vel(3,ifac) 491 if (ippmod(icompf).eq.2) then ! FIXME fill bc_frac... 492 rcodcl(ifac,isca(ifracv),1) = bc_fracv(ifac) 493 rcodcl(ifac,isca(ifracm),1) = bc_fracm(ifac) 494 rcodcl(ifac,isca(ifrace),1) = bc_frace(ifac) 495 endif 496 else ! supersonic outlet 497 rcodcl(ifac,ien,3) = 0.d0 498 rcodcl(ifac,ipr,3) = 0.d0 499 rcodcl(ifac,itk,3) = 0.d0 500 rcodcl(ifac,iu,3) = 0.d0 501 rcodcl(ifac,iv,3) = 0.d0 502 rcodcl(ifac,iw,3) = 0.d0 503 if (ippmod(icompf).eq.2) then 504 rcodcl(ifac,isca(ifracv),3) = 0.d0 505 rcodcl(ifac,isca(ifracm),3) = 0.d0 506 rcodcl(ifac,isca(ifrace),3) = 0.d0 507 endif 508 endif 509 510!=============================================================================== 511! 4.3 Boundary conditions codes (Dirichlet or Neumann) 512!=============================================================================== 513 514! P : Neumann but pressure part of momentum flux is imposed 515! as a Dirichlet BC for the pressure gradient (code 13) 516! rho, U, E, T : Dirichlet 517! k, R, eps, scal : Dirichlet/Neumann depending on the flux mass value 518 519 if (itypfb(ifac).ne.isspcf) then 520 ! Pressure : - Dirichlet for the gradient computation, allowing to have the 521 ! pressure part of the convective flux at the boundary 522 ! - Homogeneous Neumann for the diffusion 523 icodcl(ifac,ipr) = 13 524 ! velocity 525 icodcl(ifac,iu) = 1 526 icodcl(ifac,iv) = 1 527 icodcl(ifac,iw) = 1 528 ! total energy 529 icodcl(ifac,ien) = 1 530 ! temperature 531 icodcl(ifac,itk) = 1 532 ! mixture fractions 533 if (ippmod(icompf).eq.2) then 534 icodcl(ifac,isca(ifracv)) = 1 535 icodcl(ifac,isca(ifracm)) = 1 536 icodcl(ifac,isca(ifrace)) = 1 537 endif 538 else ! supersonic outlet 539 icodcl(ifac,ipr) = 3 540 icodcl(ifac,iu) = 3 541 icodcl(ifac,iv) = 3 542 icodcl(ifac,iw) = 3 543 icodcl(ifac,ien) = 3 544 icodcl(ifac,itk) = 3 545 ! mixture fractions 546 if (ippmod(icompf).eq.2) then 547 icodcl(ifac,isca(ifracv)) = 3 548 icodcl(ifac,isca(ifracm)) = 3 549 icodcl(ifac,isca(ifrace)) = 3 550 endif 551 endif 552 553!------------------------------------------------------------------------------- 554! Turbulence and passive scalars: Dirichlet / Neumann depending on the mass flux 555!------------------------------------------------------------------------------- 556 557 ! Dirichlet or homogeneous Neumann 558 ! A Dirichlet is chosen if the mass flux is ingoing and if the user provided 559 ! a value in rcodcl(ifac,ivar,1) 560 561 bmasfl = brom(ifac) & 562 *( bc_vel(1,ifac)*suffbo(1,ifac) & 563 + bc_vel(2,ifac)*suffbo(2,ifac) & 564 + bc_vel(3,ifac)*suffbo(3,ifac) ) 565 566 if (itypfb(ifac).ne.isspcf.and.bmasfl.ge.0.d0) then 567 if(itytur.eq.2) then 568 icodcl(ifac,ik ) = 3 569 icodcl(ifac,iep) = 3 570 elseif(itytur.eq.3) then 571 icodcl(ifac,ir11) = 3 572 icodcl(ifac,ir22) = 3 573 icodcl(ifac,ir33) = 3 574 icodcl(ifac,ir12) = 3 575 icodcl(ifac,ir13) = 3 576 icodcl(ifac,ir23) = 3 577 icodcl(ifac,iep ) = 3 578 elseif(iturb.eq.50) then 579 icodcl(ifac,ik ) = 3 580 icodcl(ifac,iep ) = 3 581 icodcl(ifac,iphi) = 3 582 icodcl(ifac,ifb ) = 3 583 elseif(iturb.eq.60) then 584 icodcl(ifac,ik ) = 3 585 icodcl(ifac,iomg) = 3 586 elseif(iturb.eq.70) then 587 icodcl(ifac,inusa) = 3 588 endif 589 if (nscaus.gt.0) then 590 do ii = 1, nscaus 591 icodcl(ifac,isca(ii)) = 3 592 enddo 593 endif 594 if (nscasp.gt.0) then 595 do ii = 1, nscasp 596 icodcl(ifac,iscasp(ii)) = 3 597 enddo 598 endif 599 else 600 if(itytur.eq.2) then 601 if(rcodcl(ifac,ik ,1).lt.rinfin*0.5d0 .and. & 602 rcodcl(ifac,iep,1).lt.rinfin*0.5d0) then 603 icodcl(ifac,ik ) = 1 604 icodcl(ifac,iep) = 1 605 else 606 icodcl(ifac,ik ) = 3 607 icodcl(ifac,iep) = 3 608 endif 609 elseif(itytur.eq.3) then 610 if(rcodcl(ifac,ir11,1).lt.rinfin*0.5d0 .and. & 611 rcodcl(ifac,ir22,1).lt.rinfin*0.5d0 .and. & 612 rcodcl(ifac,ir33,1).lt.rinfin*0.5d0 .and. & 613 rcodcl(ifac,ir12,1).lt.rinfin*0.5d0 .and. & 614 rcodcl(ifac,ir13,1).lt.rinfin*0.5d0 .and. & 615 rcodcl(ifac,ir23,1).lt.rinfin*0.5d0 .and. & 616 rcodcl(ifac,iep ,1).lt.rinfin*0.5d0) then 617 icodcl(ifac,ir11) = 1 618 icodcl(ifac,ir22) = 1 619 icodcl(ifac,ir33) = 1 620 icodcl(ifac,ir12) = 1 621 icodcl(ifac,ir13) = 1 622 icodcl(ifac,ir23) = 1 623 icodcl(ifac,iep ) = 1 624 else 625 icodcl(ifac,ir11) = 3 626 icodcl(ifac,ir22) = 3 627 icodcl(ifac,ir33) = 3 628 icodcl(ifac,ir12) = 3 629 icodcl(ifac,ir13) = 3 630 icodcl(ifac,ir23) = 3 631 icodcl(ifac,iep ) = 3 632 endif 633 elseif(iturb.eq.50) then 634 if(rcodcl(ifac,ik ,1).lt.rinfin*0.5d0 .and. & 635 rcodcl(ifac,iep ,1).lt.rinfin*0.5d0 .and. & 636 rcodcl(ifac,iphi,1).lt.rinfin*0.5d0 .and. & 637 rcodcl(ifac,ifb ,1).lt.rinfin*0.5d0) then 638 icodcl(ifac,ik ) = 1 639 icodcl(ifac,iep ) = 1 640 icodcl(ifac,iphi) = 1 641 icodcl(ifac,ifb ) = 1 642 else 643 icodcl(ifac,ik ) = 3 644 icodcl(ifac,iep ) = 3 645 icodcl(ifac,iphi) = 3 646 icodcl(ifac,ifb ) = 3 647 endif 648 elseif(iturb.eq.60) then 649 if(rcodcl(ifac,ik ,1).lt.rinfin*0.5d0 .and. & 650 rcodcl(ifac,iomg,1).lt.rinfin*0.5d0) then 651 icodcl(ifac,ik ) = 1 652 icodcl(ifac,iomg) = 1 653 else 654 icodcl(ifac,ik ) = 3 655 icodcl(ifac,iomg) = 3 656 endif 657 elseif(iturb.eq.70) then 658 if(rcodcl(ifac,inusa,1).gt.0.d0) then 659 icodcl(ifac,inusa) = 1 660 else 661 icodcl(ifac,inusa) = 3 662 endif 663 endif 664 if (nscaus.gt.0) then 665 do ii = 1, nscaus 666 if(rcodcl(ifac,isca(ii),1).lt.rinfin*0.5d0) then 667 icodcl(ifac,isca(ii)) = 1 668 else 669 icodcl(ifac,isca(ii)) = 3 670 endif 671 enddo 672 endif 673 if (nscasp.gt.0) then 674 do ii = 1, nscasp 675 if(rcodcl(ifac,iscasp(ii),1).lt.rinfin*0.5d0) then 676 icodcl(ifac,iscasp(ii)) = 1 677 else 678 icodcl(ifac,iscasp(ii)) = 3 679 endif 680 enddo 681 endif 682 endif 683 684 endif ! end of test on inlet/outlet faces 685 686enddo ! end of loop on boundary faces 687 688! Free memory 689deallocate(w5) 690deallocate(w7, wbfb, wbfa) 691deallocate(bc_en, bc_pr, bc_tk, bc_fracv, bc_fracm, bc_frace, bc_vel) 692 693!---- 694! FORMATS 695!---- 696 697 1000 format( & 698'@ ',/,& 699'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 700'@ ',/,& 701'@ @@ WARNING : Error during execution, ',/,& 702'@ ========= ',/,& 703'@ two and only two independant variables among ',/,& 704'@ P, rho, T and E have to be imposed at boundaries of type',/,& 705'@ iesicf in uscfcl (iccfth = ',i10,'). ',/,& 706'@ ',/,& 707'@ The computation will stop. ',/,& 708'@ ',/,& 709'@ Check the boundary condition definitions. ',/,& 710'@ ',/,& 711'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 712'@ ',/) 713 1100 format( & 714'@ ',/,& 715'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 716'@ ',/,& 717'@ @@ WARNING : Error during execution, ',/,& 718'@ ========= ',/,& 719'@ The pressure was not provided at outlet with pressure ',/,& 720'@ imposed. ',/,& 721'@ ',/,& 722'@ The computation will stop. ',/,& 723'@ ',/,& 724'@ Check the boundary conditions in ',/,& 725'@ cs_user_boundary_conditions ',/,& 726'@ ',/,& 727'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 728'@ ',/) 729 1200 format( & 730'@ ',/,& 731'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 732'@ ',/,& 733'@ @@ WARNING : Error during execution, ',/,& 734'@ ========= ',/,& 735'@ The total pressure or total enthalpy were not provided ',/,& 736'@ at inlet with total pressure and total enthalpy imposed.',/,& 737'@ ',/,& 738'@ The computation will stop. ',/,& 739'@ ',/,& 740'@ Check the boundary conditions in ',/,& 741'@ cs_user_boundary_conditions ',/,& 742'@ ',/,& 743'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 744'@ ',/) 745 1300 format( & 746'@ ',/,& 747'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 748'@ ',/,& 749'@ @@ WARNING : Error during execution, ',/,& 750'@ ========= ',/,& 751'@ The mass or enthalpy flow rate were not provided ',/,& 752'@ at inlet with mass and enthalpy flow rate imposed. ',/,& 753'@ ',/,& 754'@ The computation will stop. ',/,& 755'@ ',/,& 756'@ Check the boundary conditions in ',/,& 757'@ cs_user_boundary_conditions ',/,& 758'@ ',/,& 759'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 760'@ ',/) 761 1301 format( & 762'@ ',/,& 763'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 764'@ ',/,& 765'@ @@ WARNING : Error during execution, ',/,& 766'@ ========= ',/,& 767'@ Inlet with mass and enthalpy flow rate not provided. ',/,& 768'@ ',/,& 769'@ The computation will stop. ',/,& 770'@ ',/,& 771'@ Check the boundary conditions in ',/,& 772'@ cs_user_boundary_conditions ',/,& 773'@ ',/,& 774'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 775'@ ',/) 776 777!---- 778! END 779!---- 780 781return 782end subroutine 783