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 typecl.f90 24!> \brief Handle boundary condition type code (\ref itypfb). 25!> 26!------------------------------------------------------------------------------ 27 28!------------------------------------------------------------------------------ 29! Arguments 30!------------------------------------------------------------------------------ 31! mode name role 32!------------------------------------------------------------------------------ 33!> \param[in] nvar total number of variables 34!> \param[in] nscal total number of scalars 35!> \param[in] init partial treatment (before time loop) if true 36!> \param[in,out] itypfb boundary face types 37!> \param[out] itrifb tab d'indirection pour tri des faces 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[out] isostd standard output indicator 53!> + reference face number 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 typecl & 71 ( nvar , nscal , init , & 72 itypfb , itrifb , icodcl , isostd , & 73 rcodcl ) 74 75!=============================================================================== 76! Module files 77!=============================================================================== 78 79use paramx 80use numvar 81use optcal 82use cstnum 83use cstphy 84use dimens, only: ndimfb 85use pointe, only: b_head_loss 86use lagran, only: iilagr 87use entsor 88use parall 89use ppppar 90use ppthch 91use ppincl 92use cplsat 93use mesh 94use field 95use field_operator 96use cs_c_bindings 97 98!=============================================================================== 99 100implicit none 101 102! Arguments 103 104integer nvar , nscal 105logical init 106 107integer icodcl(ndimfb,nvar) 108integer itypfb(ndimfb) , itrifb(ndimfb) 109integer isostd(ndimfb+1) 110 111double precision rcodcl(ndimfb,nvar,3) 112 113! Local variables 114 115integer ifac, ivar, iel 116integer iok, inc, iccocg, iprev, ideb, ifin, inb, isum 117integer ityp, ii, jj, iflmab 118integer ifadir 119integer iut , ivt , iwt, ialt, iscal 120integer keyvar, keycpl 121integer iivar, icpl 122integer f_id, i_dim, f_type, nfld, f_dim, f_id_yplus, f_id_z_ground 123integer modntl 124integer kturt, turb_flux_model, turb_flux_model_type 125 126double precision pref 127double precision flumbf, flumty(ntypmx) 128double precision d0, d0min 129double precision xyzref(3) 130 131character(len=80) :: fname 132 133double precision, allocatable, dimension(:) :: pripb 134double precision, allocatable, dimension(:,:) :: grad 135double precision, dimension(:), pointer :: bmasfl 136 137double precision, pointer, dimension(:,:) :: frcxt 138double precision, dimension(:), pointer :: cvara_pr 139double precision, dimension(:), pointer :: cpro_prtot 140double precision, dimension(1,1), target :: rvoid2 141 142integer, save :: irangd 143 144integer ipass 145data ipass /0/ 146save ipass 147 148type(var_cal_opt) :: vcopt 149 150!=============================================================================== 151 152!=============================================================================== 153! 1. Initialization 154!=============================================================================== 155 156! Allocate temporary arrays 157allocate(pripb(ndimfb)) 158 159! Initialize variables to avoid compiler warnings 160pripb = 0.d0 161pref = 0.d0 162 163call field_get_val_prev_s(ivarfl(ipr), cvara_pr) 164 165call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 166 167call field_get_key_id("variable_id", keyvar) 168 169call field_get_key_id('turbulent_flux_model', kturt) 170 171call field_get_id_try("wall_yplus", f_id_yplus) 172call field_get_id_try("z_ground", f_id_z_ground) 173 174! Number of fields 175call field_get_n_fields(nfld) 176 177!=============================================================================== 178! 2. Check consistency of types given in cs_user_boundary_conditions 179!=============================================================================== 180 181iok = 0 182 183do ifac = 1, nfabor 184 ityp = itypfb(ifac) 185 if(ityp.le.0.or.ityp.gt.ntypmx) then 186 itypfb(ifac) = 0 187 iok = iok + 1 188 endif 189enddo 190 191if (irangp.ge.0) call parcmx(iok) 192 193if (iok.ne.0) then 194 call boundary_conditions_error(itypfb) 195endif 196 197!=============================================================================== 198! 3. Sort boundary faces 199!=============================================================================== 200 201 202! Count faces of each type (temporarily in ifinty) 203 204do ii = 1, ntypmx 205 ifinty(ii) = 0 206enddo 207 208do ifac = 1, nfabor 209 ityp = itypfb(ifac) 210 ifinty(ityp) = ifinty(ityp) + 1 211enddo 212 213 214! Set start of each group of faces in itrifb (sorted by type): idebty 215 216do ii = 1, ntypmx 217 idebty(ii) = 1 218enddo 219 220do ii = 1, ntypmx-1 221 do jj = ii+1, ntypmx 222 idebty(jj) = idebty(jj) + ifinty(ii) 223 enddo 224enddo 225 226! Sort faces in itrifb and use the opportunity to correctly set ifinty 227 228do ii = 1, ntypmx 229 ifinty(ii) = idebty(ii)-1 230enddo 231 232do ifac = 1, nfabor 233 ityp = itypfb(ifac) 234 ifin = ifinty(ityp)+1 235 itrifb(ifin) = ifac 236 ifinty(ityp) = ifin 237enddo 238 239! Basic check 240 241iok = 0 242do ii = 1, ntypmx-1 243 if(ifinty(ii).ge.idebty(ii+1)) then 244 if (iok.eq.0) iok = ii 245 endif 246enddo 247if (irangp.ge.0) call parcmx(iok) 248if (iok.gt.0) then 249 ii = iok 250 if(ifinty(ii).ge.idebty(ii+1)) then 251 write(nfecra,2020) (ifinty(jj),jj=1,ntypmx) 252 write(nfecra,2030) (idebty(jj),jj=1,ntypmx) 253 write(nfecra,2040) (itypfb(jj),jj=1,nfabor) 254 write(nfecra,2098) ii,ifinty(ii),ii+1,idebty(ii+1) 255 else 256 write(nfecra,2099) ii,ii+1 257 endif 258 call csexit (1) 259endif 260 261iok = 0 262isum = 0 263do ii = 1, ntypmx 264 isum = isum + ifinty(ii) - idebty(ii) + 1 265enddo 266if (irangp.ge.0) call parcpt (isum) 267if(isum.ne.nfbrgb) then 268 write(nfecra,3099) isum, nfbrgb 269 iok = iok + 1 270endif 271if (iok.ne.0) then 272 call csexit (1) 273 !========== 274endif 275 276 277! ---> On ecrit les types de faces avec la borne inf et sup et le nb 278! pour chaque type de face trouve (tjrs pour les types par defaut) 279 280if(ipass.eq.0.or.vcopt%iwarni.ge.2) then 281 282 ipass = 1 283 284 write(nfecra,6010) 285 286 write(nfecra,6011) 287 288 if ( ippmod(icompf).lt.0 ) then 289 290 ii = ientre 291 inb = ifinty(ii)-idebty(ii)+1 292 if (irangp.ge.0) call parcpt (inb) 293 write(nfecra,6020) 'Inlet ', ii, inb 294 ii = iparoi 295 inb = ifinty(ii)-idebty(ii)+1 296 if (irangp.ge.0) call parcpt (inb) 297 write(nfecra,6020) 'Smooth wall ', ii, inb 298 ii = iparug 299 inb = ifinty(ii)-idebty(ii)+1 300 if (irangp.ge.0) call parcpt (inb) 301 write(nfecra,6020) 'Rough wall ', ii, inb 302 ii = isymet 303 inb = ifinty(ii)-idebty(ii)+1 304 if (irangp.ge.0) call parcpt (inb) 305 write(nfecra,6020) 'Symmetry ', ii, inb 306 ii = isolib 307 inb = ifinty(ii)-idebty(ii)+1 308 if (irangp.ge.0) call parcpt (inb) 309 write(nfecra,6020) 'Free outlet ', ii, inb 310 ii = ifrent 311 inb = ifinty(ii)-idebty(ii)+1 312 if (irangp.ge.0) call parcpt (inb) 313 write(nfecra,6020) 'Free inlet ', ii, inb 314 315 ! Presence of free entrance face(s) 316 if (inb.gt.0) then 317 iifren = 1 318 else 319 iifren = 0 320 endif 321 if (irangp.ge.0) call parcmx(iifren) 322 323 if (iifren.eq.1) then 324 if (.not.(allocated(b_head_loss))) allocate(b_head_loss(ndimfb)) 325 endif 326 327 ii = i_convective_inlet 328 inb = ifinty(ii)-idebty(ii)+1 329 if (irangp.ge.0) call parcpt (inb) 330 write(nfecra,6020) 'Convective inlet ', ii, inb 331 332 if (nbrcpl.ge.1) then 333 if (ifaccp.eq.0) then 334 ii = icscpl 335 else 336 ii = icscpd 337 endif 338 inb = ifinty(ii)-idebty(ii)+1 339 if (irangp.ge.0) call parcpt (inb) 340 write(nfecra,6020) 'Sat/Sat coupling ', ii, inb 341 endif 342 343 ii = ifresf 344 inb = ifinty(ii)-idebty(ii)+1 345 if (irangp.ge.0) call parcpt (inb) 346 write(nfecra,6020) 'Free surface ', ii, inb 347 348 ii = iindef 349 inb = ifinty(ii)-idebty(ii)+1 350 if (irangp.ge.0) call parcpt (inb) 351 write(nfecra,6020) 'Undefined ', ii, inb 352 353 do ii = 1, ntypmx 354 if (ii.ne.ientre .and. & 355 ii.ne.i_convective_inlet .and. & 356 ii.ne.iparoi .and. & 357 ii.ne.iparug .and. & 358 ii.ne.isymet .and. & 359 ii.ne.isolib .and. & 360 ii.ne.ifrent .and. & 361 ii.ne.ifresf .and. & 362 ii.ne.icscpl .and. & 363 ii.ne.icscpd .and. & 364 ii.ne.iindef ) then 365 inb = ifinty(ii)-idebty(ii)+1 366 if (irangp.ge.0) call parcpt (inb) 367 if(inb.gt.0) then 368 write(nfecra,6020) 'User type ', ii, inb 369 endif 370 endif 371 enddo 372 373 else 374 375 ii = ieqhcf 376 inb = ifinty(ii)-idebty(ii)+1 377 if (irangp.ge.0) call parcpt (inb) 378 write(nfecra,6020) 'Sub. enth. inlet ', ii, inb 379 380 ii = iephcf 381 inb = ifinty(ii)-idebty(ii)+1 382 if (irangp.ge.0) call parcpt (inb) 383 write(nfecra,6020) 'Ptot, Htot ', ii, inb 384 385 ii = iesicf 386 inb = ifinty(ii)-idebty(ii)+1 387 if (irangp.ge.0) call parcpt (inb) 388 write(nfecra,6020) 'Imp inlet/outlet ', ii, inb 389 390 ii = isopcf 391 inb = ifinty(ii)-idebty(ii)+1 392 if (irangp.ge.0) call parcpt (inb) 393 write(nfecra,6020) 'Subsonic outlet ', ii, inb 394 395 ii = isspcf 396 inb = ifinty(ii)-idebty(ii)+1 397 if (irangp.ge.0) call parcpt (inb) 398 write(nfecra,6020) 'Supersonic outlet', ii, inb 399 400 ii = iparoi 401 inb = ifinty(ii)-idebty(ii)+1 402 if (irangp.ge.0) call parcpt (inb) 403 write(nfecra,6020) 'Smooth wall ', ii, inb 404 405 ii = iparug 406 inb = ifinty(ii)-idebty(ii)+1 407 if (irangp.ge.0) call parcpt (inb) 408 write(nfecra,6020) 'Rough wall ', ii, inb 409 410 ii = isymet 411 inb = ifinty(ii)-idebty(ii)+1 412 if (irangp.ge.0) call parcpt (inb) 413 write(nfecra,6020) 'Symmetry ', ii, inb 414 415 ii = iindef 416 inb = ifinty(ii)-idebty(ii)+1 417 if (irangp.ge.0) call parcpt (inb) 418 write(nfecra,6020) 'Undefined ', ii, inb 419 420 do ii = 1, ntypmx 421 if (ii.ne.iesicf .and. & 422 ii.ne.isspcf .and. & 423 ii.ne.ieqhcf .and. & 424 ii.ne.iephcf .and. & 425 ii.ne.isopcf .and. & 426 ii.ne.iparoi .and. & 427 ii.ne.iparug .and. & 428 ii.ne.isymet .and. & 429 ii.ne.iindef ) then 430 inb = ifinty(ii)-idebty(ii)+1 431 if (irangp.ge.0) call parcpt (inb) 432 if(inb.gt.0) then 433 write(nfecra,6020) 'User type ',ii, inb 434 endif 435 endif 436 enddo 437 438 endif 439 440 write(nfecra,6030) 441 442endif 443 444!================================================================================ 445! 4. rcodcl(., .,1) has been initialized as rinfin so as to check what 446! the user has modified. Those not modified are reset to zero here. 447! isolib, ifrent and ientre are handled later. 448!================================================================================ 449 450do ivar = 1, nvar 451 do ifac = 1, nfabor 452 if ((itypfb(ifac) .ne. isolib) .and. & 453 (itypfb(ifac) .ne. ifrent) .and. & 454 (itypfb(ifac) .ne. ifresf) .and. & 455 (itypfb(ifac) .ne. i_convective_inlet).and. & 456 (itypfb(ifac) .ne. ientre) .and. & 457 (rcodcl(ifac,ivar,1) .gt. rinfin*0.5d0)) then 458 rcodcl(ifac,ivar,1) = 0.d0 459 endif 460 enddo 461enddo 462 463!=============================================================================== 464! 5. Compute pressure at boundary (in pripb(*)) 465! (if we need it, that is if there are outlet boudary faces). 466!=============================================================================== 467 468! ifrslb = closest free standard outlet face to xyzp0 (icodcl not modified) 469! (or closest free inlet) 470! itbslb = max of ifrslb on all ranks, standard outlet face presence indicator 471 472! Even when the user has not chosen xyzp0 (and it is thus at the 473! origin), we choose the face whose center is closest to it, so 474! as to be mesh numbering (and partitioning) independent. 475 476if (ntcabs.eq.ntpabs+1) then 477 478 d0min = rinfin 479 480 ideb = idebty(isolib) 481 ifin = ifinty(isolib) 482 483 do ii = ideb, ifin 484 ifac = itrifb(ii) 485 if (icodcl(ifac,ipr).eq.0) then 486 d0 = (cdgfbo(1,ifac)-xyzp0(1))**2 & 487 + (cdgfbo(2,ifac)-xyzp0(2))**2 & 488 + (cdgfbo(3,ifac)-xyzp0(3))**2 489 if (d0.lt.d0min) then 490 ifrslb = ifac 491 d0min = d0 492 endif 493 endif 494 enddo 495 496 ideb = idebty(ifrent) 497 ifin = ifinty(ifrent) 498 499 do ii = ideb, ifin 500 ifac = itrifb(ii) 501 if (icodcl(ifac,ipr).eq.0) then 502 d0 = (cdgfbo(1,ifac)-xyzp0(1))**2 & 503 + (cdgfbo(2,ifac)-xyzp0(2))**2 & 504 + (cdgfbo(3,ifac)-xyzp0(3))**2 505 if (d0.lt.d0min) then 506 ifrslb = ifac 507 d0min = d0 508 endif 509 endif 510 enddo 511 512 ! If we have free outlet faces, irangd and itbslb will 513 ! contain respectively the rank having the boundary face whose 514 ! center is closest to xyzp0, and the local number of that face 515 ! on that rank (also equal to ifrslb on that rank). 516 ! If we do not have free outlet faces, than itbslb = 0 517 ! (as it was initialized that way on all ranks). 518 519 itbslb = ifrslb 520 irangd = irangp 521 if (irangp.ge.0) then 522 call parfpt(itbslb, irangd, d0min) 523 endif 524 525endif 526 527if (iphydr.eq.1.or.iifren.eq.1) then 528 529! En cas de prise en compte de la pression hydrostatique, 530! on remplit le tableau ISOSTD 531! 0 -> pas une face de sortie standard (i.e. pas sortie ou sortie avec CL 532! de pression modifiee) 533! 1 -> face de sortie libre avec CL de pression automatique. 534! le numero de la face de reference est stocke dans ISOSTD(NFABOR+1) 535! qui est d'abord initialise a -1 (i.e. pas de face de sortie std) 536 isostd(nfabor+1) = -1 537 do ifac = 1,nfabor 538 isostd(ifac) = 0 539 if ((itypfb(ifac).eq.isolib.or.itypfb(ifac).eq.ifrent).and. & 540 (icodcl(ifac,ipr).eq.0)) then 541 isostd(ifac) = 1 542 endif 543 enddo 544endif 545 546! ---> Reference pressure (unique, even if there are multiple outlets) 547! In case we account for the hydrostatic pressure, we search for the 548! reference face. 549 550! Determine a unique P I' pressure in parallel 551! if there are free outlet faces, we have determined that the rank 552! with the outlet face closest to xyzp0 is irangd. 553 554! We also retrieve the coordinates of the reference point, so as to 555! calculate pref later on. 556 557if (itbslb.gt.0) then 558 559 ! If irangd is the local rank, we assign PI' to xyzref(4) 560 ! (this is always the case in serial mode) 561 562 if (irangp.eq.irangd) then 563 xyzref(1) = cdgfbo(1,ifrslb) 564 xyzref(2) = cdgfbo(2,ifrslb) 565 xyzref(3) = cdgfbo(3,ifrslb) 566 if (iphydr.eq.1.or.iifren.eq.1) isostd(nfabor+1) = ifrslb 567 endif 568 569 ! Broadcast PI' and pressure reference 570 ! from irangd to all other ranks. 571 if (irangp.ge.0) then 572 inb = 3 573 call parbcr(irangd, inb, xyzref) 574 endif 575 576 ! If the user has not specified anything, we set ixyzp0 to 2 so as 577 ! to update the reference point. 578 579 if (ixyzp0.eq.-1) ixyzp0 = 2 580 581elseif (ixyzp0.lt.0.and.ntcabs.eq.ntpabs+1) then 582 583 ! If there are no outlet faces, we search for possible Dirichlets 584 ! specified by the user so as to locate the reference point. 585 ! As before, we chose the face closest to xyzp0 so as to 586 ! be mesh numbering (and partitioning) independent. 587 588 d0min = rinfin 589 590 ifadir = -1 591 do ifac = 1, nfabor 592 if (abs(icodcl(ifac,ipr)).eq.1) then 593 d0 = (cdgfbo(1,ifac)-xyzp0(1))**2 & 594 + (cdgfbo(2,ifac)-xyzp0(2))**2 & 595 + (cdgfbo(3,ifac)-xyzp0(3))**2 596 if (d0.lt.d0min) then 597 ifadir = ifac 598 d0min = d0 599 endif 600 endif 601 enddo 602 603 irangd = irangp 604 if (irangp.ge.0) call parfpt(ifadir, irangd, d0min) 605 606 if (ifadir.gt.0) then 607 608 ! on met ixyzp0 a 2 pour mettre a jour le point de reference 609 ixyzp0 = 2 610 611 if (irangp.eq.irangd) then 612 xyzref(1) = cdgfbo(1,ifadir) 613 xyzref(2) = cdgfbo(2,ifadir) 614 xyzref(3) = cdgfbo(3,ifadir) 615 endif 616 617 ! Broadcast xyzref from irangd to all other ranks. 618 if (irangp.ge.0) then 619 inb = 3 620 call parbcr(irangd, inb, xyzref) 621 endif 622 623 endif 624 625endif 626 627! Si le point de reference n'a pas ete specifie par l'utilisateur 628! on le change et on decale alors COEFU s'il y a des sorties. 629! La pression totale est aussi decalee (c'est a priori 630! inutile sauf si l'utilisateur l'utilise dans ustsnv par exemple) 631 632if (ixyzp0.eq.2) then 633 ixyzp0 = 1 634 if (ippmod(icompf).lt.0.and.ippmod(idarcy).lt.0) then 635 call field_get_val_s(iprtot, cpro_prtot) 636 do iel = 1, ncelet 637 cpro_prtot(iel) = cpro_prtot(iel) - ro0*( gx*(xyzref(1) - xyzp0(1)) & 638 + gy*(xyzref(2) - xyzp0(2)) & 639 + gz*(xyzref(3) - xyzp0(3))) 640 enddo 641 endif 642 643 xyzp0(1) = xyzref(1) 644 xyzp0(2) = xyzref(2) 645 xyzp0(3) = xyzref(3) 646 647 if (itbslb.gt.0) then 648 write(nfecra,8000) xyzp0(1), xyzp0(2), xyzp0(3) 649 else 650 write(nfecra,8001) xyzp0(1), xyzp0(2), xyzp0(3) 651 endif 652 653elseif (ixyzp0.eq.-1) then 654 ! Il n'y a pas de sorties ni de Dirichlet et l'utilisateur n'a 655 ! rien specifie -> on met IXYZP0 a 0 pour ne plus y toucher, tout 656 ! en differenciant du cas =1 qui necessitera une ecriture en suite 657 ixyzp0 = 0 658endif 659 660! No need of computing pressure gradient for frozen field computations 661if (itbslb.gt.0.and.iilagr.ne.3) then 662 663 if (iphydr.eq.1) then 664 call field_get_val_v_by_name('volume_forces', frcxt) 665 else 666 frcxt => rvoid2 667 endif 668 669 ! Allocate a work array for the gradient calculation 670 allocate(grad(3,ncelet)) 671 672 iprev = 1 673 inc = 1 674 iccocg = 1 675 676 call grdpor(inc) 677 678 call field_gradient_potential(ivarfl(ipr), iprev, 0, inc, & 679 iccocg, iphydr, & 680 frcxt, grad) 681 682 ! Put in pripb the value at F of the 683 ! total pressure, computed from P* shifted from the ro0*g.(X-Xref) 684 685 do ifac = 1, nfabor 686 ii = ifabor(ifac) 687 pripb(ifac) = cvara_pr(ii) & 688 + (cdgfbo(1,ifac)-xyzcen(1,ii))*grad(1,ii) & 689 + (cdgfbo(2,ifac)-xyzcen(2,ii))*grad(2,ii) & 690 + (cdgfbo(3,ifac)-xyzcen(3,ii))*grad(3,ii) 691 enddo 692 693 ! Free memory 694 deallocate(grad) 695 696 if (irangp.eq.irangd) pref = pripb(ifrslb) 697 if (irangp.ge.0) then 698 call parall_bcast_r(irangd, pref) 699 endif 700 701endif 702 703!=============================================================================== 704! 6. Convert to rcodcl and icodcl 705! (if this has not already been set by the user) 706 707! First, process variables for which a specific treatement is done 708! (pressure, velocity, ...) 709!=============================================================================== 710 711 712! Translate total pressure P_tot given by the user into solved pressure P 713! P = P_tot - p0 - rho0 ( g . (z-z0)) 714! 715! If icodcl = -1, then the BC Dirichlet value is given in solved pressure P, 716! no need of transformation from P_tot to P 717 718ivar = ipr 719if (ippmod(icompf).lt.0.and.ippmod(idarcy).lt.0) then 720 do ifac = 1, nfabor 721 if (icodcl(ifac,ivar).eq.-1) then 722 icodcl(ifac,ivar) = 1 723 else if (icodcl(ifac,ivar).ne.0) then 724 rcodcl(ifac,ivar,1) = rcodcl(ifac,ivar,1) & 725 - ro0*( gx*(cdgfbo(1,ifac) - xyzp0(1)) & 726 + gy*(cdgfbo(2,ifac) - xyzp0(2)) & 727 + gz*(cdgfbo(3,ifac) - xyzp0(3))) & 728 - p0 729 endif 730 enddo 731endif 732 733! 6.1.1 Inlet 734! =========== 735 736! ---> La pression a un traitement Neumann, le reste Dirichlet 737! sera traite plus tard. 738 739ideb = idebty(ientre) 740ifin = ifinty(ientre) 741 742ivar = ipr 743do ii = ideb, ifin 744 ifac = itrifb(ii) 745 if (icodcl(ifac,ivar).eq.0) then 746 icodcl(ifac,ivar) = 3 747 rcodcl(ifac,ivar,1) = 0.d0 748 rcodcl(ifac,ivar,2) = rinfin 749 rcodcl(ifac,ivar,3) = 0.d0 750 endif 751enddo 752 753! 6.1.2 Convective Inlet 754! ====================== 755 756! ---> La pression a un traitement Neumann, le reste Dirichlet 757! sera traite plus tard. 758 759ideb = idebty(i_convective_inlet) 760ifin = ifinty(i_convective_inlet) 761 762ivar = ipr 763do ii = ideb, ifin 764 ifac = itrifb(ii) 765 if (icodcl(ifac,ivar).eq.0) then 766 icodcl(ifac,ivar) = 3 767 rcodcl(ifac,ivar,1) = 0.d0 768 rcodcl(ifac,ivar,2) = rinfin 769 rcodcl(ifac,ivar,3) = 0.d0 770 endif 771enddo 772 773 774! 6.2 SORTIE (entree-sortie libre) (ISOLIB) 775! =================== 776 777! ---> La pression a un traitement Dirichlet, les vitesses 9 778! (le reste Neumann, ou Dirichlet si donnee utilisateur, 779! sera traite plus tard) 780 781! ---> Entree/Sortie libre 782 783ideb = idebty(isolib) 784ifin = ifinty(isolib) 785 786do ivar = 1, nvar 787 if (ivar.eq.ipr) then 788 do ii = ideb, ifin 789 ifac = itrifb(ii) 790 if(icodcl(ifac,ivar).eq.0) then 791 icodcl(ifac,ivar) = 1 792 rcodcl(ifac,ivar,1) = pripb(ifac) - pref 793 rcodcl(ifac,ivar,2) = rinfin 794 rcodcl(ifac,ivar,3) = 0.d0 795 endif 796 enddo 797 elseif(ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 798 do ii = ideb, ifin 799 ifac = itrifb(ii) 800 if(icodcl(ifac,ivar).eq.0) then 801 icodcl(ifac,ivar) = 9 802 rcodcl(ifac,ivar,1) = 0.d0 803 rcodcl(ifac,ivar,2) = rinfin 804 rcodcl(ifac,ivar,3) = 0.d0 805 endif 806 enddo 807 endif 808enddo 809 810! ---> Free entrance (Bernoulli relation), std free outlet 811! (a specific treatment is performed on the increment of pressure) 812 813ideb = idebty(ifrent) 814ifin = ifinty(ifrent) 815 816do ivar = 1, nvar 817 if (ivar.eq.ipr) then 818 do ii = ideb, ifin 819 ifac = itrifb(ii) 820 iel = ifabor(ifac) 821 822 if (icodcl(ifac,ivar).eq.0) then 823 824 ! If the user has given a value of boundary head loss 825 if (rcodcl(ifac,ivar,2).le.0.5d0*rinfin) then 826 b_head_loss(ifac) = rcodcl(ifac,ivar,2) 827 else 828 b_head_loss(ifac) = 0.d0 829 endif 830 831 ! Std outlet 832 icodcl(ifac,ivar) = 1 833 rcodcl(ifac,ivar,1) = pripb(ifac) - pref 834 rcodcl(ifac,ivar,2) = rinfin 835 rcodcl(ifac,ivar,3) = 0.d0 836 837 endif 838 enddo 839 elseif(ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 840 do ii = ideb, ifin 841 ifac = itrifb(ii) 842 ! Homogeneous Neumann 843 if(icodcl(ifac,ivar).eq.0) then 844 icodcl(ifac,ivar) = 3 845 rcodcl(ifac,ivar,1) = 0.d0 846 rcodcl(ifac,ivar,2) = rinfin 847 rcodcl(ifac,ivar,3) = 0.d0 848 endif 849 enddo 850 endif 851enddo 852 853 854! ---> Free surface: Dirichlet on the pressure 855 856ideb = idebty(ifresf) 857ifin = ifinty(ifresf) 858 859do ivar = 1, nvar 860 if (ivar.eq.ipr) then 861 do ii = ideb, ifin 862 ifac = itrifb(ii) 863 if(icodcl(ifac,ivar).eq.0) then 864 icodcl(ifac,ivar) = 1 865 rcodcl(ifac,ivar,1) = - ro0*( gx*(cdgfbo(1,ifac) - xyzp0(1)) & 866 + gy*(cdgfbo(2,ifac) - xyzp0(2)) & 867 + gz*(cdgfbo(3,ifac) - xyzp0(3))) 868 rcodcl(ifac,ivar,2) = rinfin 869 rcodcl(ifac,ivar,3) = 0.d0 870 endif 871 enddo 872 elseif(ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 873 do ii = ideb, ifin 874 ifac = itrifb(ii) 875 ! Homogeneous Neumann 876 if(icodcl(ifac,ivar).eq.0) then 877 icodcl(ifac,ivar) = 3 878 rcodcl(ifac,ivar,1) = 0.d0 879 rcodcl(ifac,ivar,2) = rinfin 880 rcodcl(ifac,ivar,3) = 0.d0 881 endif 882 enddo 883 endif 884enddo 885 886! Free memory 887deallocate(pripb) 888 889! Symmetry 890! ========= 891 892! ---> Vectors and tensors have a special treatment. 893! Scalars have an homogeneous Neumann. 894 895ideb = idebty(isymet) 896ifin = ifinty(isymet) 897 898do f_id = 0, nfld - 1 899 900 call field_get_type(f_id, f_type) 901 902 ! Is the field of type FIELD_VARIABLE? 903 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 904 ! Is this field not managed by CDO 905 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 906 907 call field_get_dim (f_id, f_dim) 908 call field_get_key_int(f_id, keyvar, ivar) 909 910 ! Loop over faces 911 do ii = ideb, ifin 912 ifac = itrifb(ii) 913 ! Special treatment for uncoupled version of Rij models, 914 ! will be removed with the coupled solver 915 if (icodcl(ifac,ivar).eq.0.and.irijco.eq.0.and. & 916 (ivar.eq.ir11.or.ivar.eq.ir22.or.ivar.eq.ir33.or. & 917 ivar.eq.ir12.or.ivar.eq.ir13.or.ivar.eq.ir23)) then 918 icodcl(ifac,ivar) = 4 919 rcodcl(ifac,ivar,1) = 0.d0 920 rcodcl(ifac,ivar,2) = rinfin 921 rcodcl(ifac,ivar,3) = 0.d0 922 endif 923 924 ! Homogeneous Neumann on scalars 925 if (f_dim.eq.1.and.icodcl(ifac,ivar).eq.0) then 926 icodcl(ifac,ivar) = 3 927 rcodcl(ifac,ivar,1) = 0.d0 928 rcodcl(ifac,ivar,2) = rinfin 929 rcodcl(ifac,ivar,3) = 0.d0 930 931 ! Symmetry BC if nothing is set by the user on vector and tensors 932 else if (icodcl(ifac,ivar).eq.0) then 933 do i_dim = 0, f_dim - 1 934 icodcl(ifac,ivar + i_dim) = 4 935 rcodcl(ifac,ivar + i_dim, 1) = 0.d0 936 rcodcl(ifac,ivar + i_dim, 2) = rinfin 937 rcodcl(ifac,ivar + i_dim, 3) = 0.d0 938 enddo 939 endif 940 enddo 941 endif 942 endif 943enddo 944 945! 6.4 Smooth wall 946! =============== 947 948! ---> Velocity and turbulent quantities have icodcl = 5 949! Turbulent fluxes of scalars has 0 Dirichlet if scalars 950! have a Dirichlet, otherwise treated in clptur. 951! Other quantities are treated afterwards (Homogeneous Neumann) 952 953ideb = idebty(iparoi) 954ifin = ifinty(iparoi) 955 956do ivar = 1, nvar 957 if (ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 958 do ii = ideb, ifin 959 ifac = itrifb(ii) 960 if (icodcl(ifac,ivar).eq.0) then 961 icodcl(ifac,ivar) = 5 962 ! rcodcl(ifac,ivar,1) = User 963 rcodcl(ifac,ivar,2) = rinfin 964 ! rcodcl(ifac,ivar,3) unused value, let it as default 965 endif 966 enddo 967 elseif ((itytur.eq.2 & 968 .and.(ivar.eq.ik.or.ivar.eq.iep)) & 969 .or.(itytur.eq.3 & 970 .and.(ivar.eq.ir11.or.ivar.eq.ir22.or.ivar.eq.ir33.or. & 971 ivar.eq.ir12.or.ivar.eq.ir13.or.ivar.eq.ir23.or. & 972 ivar.eq.iep.or.ivar.eq.ial)) & 973 .or.(iturb.eq.50 & 974 .and.(ivar.eq.ik.or.ivar.eq.iep.or.ivar.eq.iphi.or. & 975 ivar.eq.ifb)) & 976 .or.(iturb.eq.51 & 977 .and.(ivar.eq.ik.or.ivar.eq.iep.or.ivar.eq.iphi.or. & 978 ivar.eq.ial)) & 979 .or.(iturb.eq.60 & 980 .and.(ivar.eq.ik.or.ivar.eq.iomg)) & 981 .or.(iturb.eq.70 & 982 .and.(ivar.eq.inusa))) then 983 do ii = ideb, ifin 984 ifac = itrifb(ii) 985 if(icodcl(ifac,ivar).eq.0) then 986 icodcl(ifac,ivar) = 5 987 rcodcl(ifac,ivar,1) = 0.d0 988 rcodcl(ifac,ivar,2) = rinfin 989 rcodcl(ifac,ivar,3) = 0.d0 990 endif 991 enddo 992 993 ! Homogeneous Neumann on the pressure 994 elseif(ivar.eq.ipr) then 995 do ii = ideb, ifin 996 ifac = itrifb(ii) 997 if (icodcl(ifac,ivar).eq.0) then 998 icodcl(ifac,ivar) = 3 999 rcodcl(ifac,ivar,1) = 0.d0 1000 rcodcl(ifac,ivar,2) = rinfin 1001 rcodcl(ifac,ivar,3) = 0.d0 1002 endif 1003 enddo 1004 endif 1005enddo 1006 1007! Turbulent fluxes 1008do iscal = 1, nscal 1009 call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 1010 turb_flux_model_type = turb_flux_model / 10 1011 1012 if (turb_flux_model_type.eq.3) then 1013 ! Name of the scalar ivar 1014 call field_get_name(ivarfl(isca(iscal)), fname) 1015 1016 ! Index of the corresponding turbulent flux 1017 call field_get_id(trim(fname)//'_turbulent_flux', f_id) 1018 1019 ! Set pointer values of turbulent fluxes in icodcl 1020 call field_get_key_int(f_id, keyvar, iut) 1021 ivt = iut + 1 1022 iwt = iut + 2 1023 1024 do ii = ideb, ifin 1025 ifac = itrifb(ii) 1026 if (icodcl(ifac, iut).eq.0) then 1027 icodcl(ifac,iut) = 5 1028 icodcl(ifac,ivt) = 5 1029 icodcl(ifac,iwt) = 5 1030 rcodcl(ifac,iut,1) = 0.d0 1031 rcodcl(ifac,ivt,1) = 0.d0 1032 rcodcl(ifac,iwt,1) = 0.d0 1033 endif 1034 enddo 1035 endif 1036 1037 ! EB-GGDH/AFM/DFM alpha boundary conditions 1038 if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then 1039 ! Name of the scalar ivar 1040 call field_get_name(ivarfl(isca(iscal)), fname) 1041 1042 ! Index of the corresponding turbulent flux 1043 call field_get_id(trim(fname)//'_alpha', f_id) 1044 1045 ! Set pointer values of turbulent fluxes in icodcl 1046 call field_get_key_int(f_id, keyvar, ialt) 1047 1048 do ii = ideb, ifin 1049 ifac = itrifb(ii) 1050 if (icodcl(ifac, ialt).eq.0) then 1051 icodcl(ifac,ialt) = 5 1052 rcodcl(ifac,ialt,1) = 0.d0 1053 endif 1054 enddo 1055 endif 1056 1057! End loop over scalars 1058enddo 1059 1060 1061! 6.5 PAROI RUGUEUSE 1062! ================== 1063 1064! ---> La vitesse et les grandeurs turbulentes ont le code 6 1065! la rugosite est stockee dans rcodcl(..,..,3) 1066! le reste Neumann sera traite plus tard (idem paroi lisse) 1067 1068ideb = idebty(iparug) 1069ifin = ifinty(iparug) 1070 1071do ivar = 1, nvar 1072 if ( ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 1073 do ii = ideb, ifin 1074 ifac = itrifb(ii) 1075 if(icodcl(ifac,ivar).eq.0) then 1076 icodcl(ifac,ivar) = 6 1077 ! rcodcl(ifac,ivar,1) = Utilisateur 1078 rcodcl(ifac,ivar,2) = rinfin 1079 ! rcodcl(ifac,ivar,3) = Utilisateur 1080 endif 1081 enddo 1082 elseif ((itytur.eq.2 & 1083 .and.(ivar.eq.ik.or.ivar.eq.iep)) & 1084 .or.(itytur.eq.3 & 1085 .and.(ivar.eq.ir11.or.ivar.eq.ir22.or.ivar.eq.ir33.or. & 1086 ivar.eq.ir12.or.ivar.eq.ir13.or.ivar.eq.ir23.or. & 1087 ivar.eq.iep.or.ivar.eq.ial)) & 1088 .or.(iturb.eq.50 & 1089 .and.(ivar.eq.ik.or.ivar.eq.iep.or.ivar.eq.iphi.or. & 1090 ivar.eq.ifb)) & 1091 .or.(iturb.eq.51 & 1092 .and.(ivar.eq.ik.or.ivar.eq.iep.or.ivar.eq.iphi.or. & 1093 ivar.eq.ial)) & 1094 .or.(iturb.eq.60 & 1095 .and.(ivar.eq.ik.or.ivar.eq.iomg)) & 1096 .or.(iturb.eq.70 & 1097 .and.(ivar.eq.inusa))) then 1098 do ii = ideb, ifin 1099 ifac = itrifb(ii) 1100 if(icodcl(ifac,ivar).eq.0) then 1101 icodcl(ifac,ivar) = 6 1102 rcodcl(ifac,ivar,1) = 0.d0 1103 rcodcl(ifac,ivar,2) = rinfin 1104 rcodcl(ifac,ivar,3) = 0.d0 1105 endif 1106 enddo 1107 elseif(ivar.eq.ipr) then 1108 do ii = ideb, ifin 1109 ifac = itrifb(ii) 1110 if(icodcl(ifac,ivar).eq.0) then 1111 icodcl(ifac,ivar) = 3 1112 rcodcl(ifac,ivar,1) = 0.d0 1113 rcodcl(ifac,ivar,2) = rinfin 1114 rcodcl(ifac,ivar,3) = 0.d0 1115 endif 1116 enddo 1117 endif 1118enddo 1119 1120! Turbulent fluxes 1121do iscal = 1, nscal 1122 call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 1123 turb_flux_model_type = turb_flux_model / 10 1124 1125 if (turb_flux_model_type.eq.3) then 1126 ! Name of the scalar ivar 1127 call field_get_name(ivarfl(isca(iscal)), fname) 1128 1129 ! Index of the corresponding turbulent flux 1130 call field_get_id(trim(fname)//'_turbulent_flux', f_id) 1131 1132 ! Set pointer values of turbulent fluxes in icodcl 1133 call field_get_key_int(f_id, keyvar, iut) 1134 ivt = iut + 1 1135 iwt = iut + 2 1136 1137 do ii = ideb, ifin 1138 ifac = itrifb(ii) 1139 if (icodcl(ifac, iut).eq.0) then 1140 icodcl(ifac,iut) = 6 1141 icodcl(ifac,ivt) = 6 1142 icodcl(ifac,iwt) = 6 1143 rcodcl(ifac,iut,1) = 0.d0 1144 rcodcl(ifac,ivt,1) = 0.d0 1145 rcodcl(ifac,iwt,1) = 0.d0 1146 endif 1147 enddo 1148 endif 1149 1150 if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then 1151 ! Name of the scalar ivar 1152 call field_get_name(ivarfl(isca(iscal)), fname) 1153 1154 ! Index of the corresponding turbulent flux 1155 call field_get_id(trim(fname)//'_alpha', f_id) 1156 1157 ! Set pointer values of turbulent fluxes in icodcl 1158 call field_get_key_int(f_id, keyvar, ialt) 1159 1160 do ii = ideb, ifin 1161 ifac = itrifb(ii) 1162 if (icodcl(ifac, ialt).eq.0) then 1163 icodcl(ifac,ialt) = 6 1164 rcodcl(ifac,ialt,1) = 0.d0 1165 endif 1166 enddo 1167 endif 1168 1169! End loop over scalars 1170enddo 1171 1172! When called before time loop, some values are not yet available. 1173if (init .eqv. .true.) return 1174 1175!=============================================================================== 1176! 6.bis CONVERSION EN RCODCL ICODCL 1177! (SI CE DERNIER N'A PAS DEJA ETE RENSEIGNE PAR L'UTILISATEUR) 1178 1179! MAINTENANT LES VARIABLES POUR LESQUELLES IL N'EXISTE PAS DE 1180! TRAITEMENT PARTICULIER (HORS PRESSION, VITESSE ...) 1181!=============================================================================== 1182 1183! 6.1.1 Inlet bis 1184! =============== 1185 1186! ---> Pressure is already treated (with a Neumann BC) 1187! Dirichlet BC for the velocity. 1188! Dirichlet BC for scalars if the user give a value, otherwise homogeneous 1189! Neumann if the mass flux is outgoing. 1190 1191ideb = idebty(ientre) 1192ifin = ifinty(ientre) 1193 1194iok = 0 1195iivar = 0 1196do ivar = 1, nvar 1197 1198 ! Convective mass flux of the variable 1199 call field_get_key_int(ivarfl(ivar), kbmasf, iflmab) 1200 call field_get_val_s(iflmab, bmasfl) 1201 1202 do ii = ideb, ifin 1203 ifac = itrifb(ii) 1204 if(icodcl(ifac,ivar).eq.0) then 1205 1206 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1207 1208 ! For not convected variables, if nothing is defined: homogeneous Neumann 1209 if (vcopt%iconv .eq. 0 .and. rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1210 icodcl(ifac,ivar) = 3 1211 rcodcl(ifac,ivar,1) = 0.d0 1212 rcodcl(ifac,ivar,2) = rinfin 1213 rcodcl(ifac,ivar,3) = 0.d0 1214 else if (ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 1215 if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1216 itypfb(ifac) = - abs(itypfb(ifac)) 1217 if (iok.eq.0.or.iok.eq.2) then 1218 iok = iok + 1 1219 endif 1220 else 1221 icodcl(ifac,ivar) = 1 1222 ! rcodcl(ifac,ivar,1) = given by the user 1223 rcodcl(ifac,ivar,2) = rinfin 1224 rcodcl(ifac,ivar,3) = 0.d0 1225 endif 1226 1227 else if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1228 1229 flumbf = bmasfl(ifac) 1230 ! Outgoing flux or yplus or z_ground 1231 if (flumbf.ge.-epzero.or. ivarfl(ivar).eq.f_id_yplus & 1232 .or.ivarfl(ivar).eq.f_id_z_ground) then 1233 icodcl(ifac,ivar) = 3 1234 rcodcl(ifac,ivar,1) = 0.d0 1235 rcodcl(ifac,ivar,2) = rinfin 1236 rcodcl(ifac,ivar,3) = 0.d0 1237 else 1238 itypfb(ifac) = - abs(itypfb(ifac)) 1239 if (iok.lt.2) then 1240 iok = iok + 2 1241 iivar = max(iivar, ivar) 1242 endif 1243 endif 1244 else 1245 icodcl(ifac,ivar) = 1 1246 ! rcodcl(ifac,ivar,1) = given by the user 1247 rcodcl(ifac,ivar,2) = rinfin 1248 rcodcl(ifac,ivar,3) = 0.d0 1249 endif 1250 1251 endif 1252 enddo 1253enddo 1254 1255if (irangp.ge.0) then 1256 call parcmx(iok) 1257 call parcmx(iivar) 1258endif 1259if (iok.gt.0) then 1260 if (iok.eq.1 .or. iok.eq.3) write(nfecra,6060) 1261 if (iok.eq.2 .or. iok.eq.3) then 1262 call field_get_name(ivarfl(iivar), fname) 1263 write(nfecra,6070) trim(fname) 1264 endif 1265 call boundary_conditions_error(itypfb) 1266endif 1267 1268! 6.1.1 Convective Inlet bis 1269! ========================== 1270 1271! ---> Pressure is already treated (with a Neumann BC) 1272! Dirichlet BC for the velocity. 1273! Dirichlet BC for scalars if the user give a value, otherwise homogeneous 1274! Neumann if the mass flux is outgoing. 1275 1276ideb = idebty(i_convective_inlet) 1277ifin = ifinty(i_convective_inlet) 1278 1279iok = 0 1280iivar = 0 1281do ivar = 1, nvar 1282 1283 ! Convective mass flux of the variable 1284 call field_get_key_int(ivarfl(ivar), kbmasf, iflmab) 1285 call field_get_val_s(iflmab, bmasfl) 1286 1287 do ii = ideb, ifin 1288 ifac = itrifb(ii) 1289 if(icodcl(ifac,ivar).eq.0) then 1290 1291 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1292 1293 ! For not convected variables, if nothing is defined: homogeneous Neumann 1294 if (vcopt%iconv .eq. 0 .and. rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1295 icodcl(ifac,ivar) = 3 1296 rcodcl(ifac,ivar,1) = 0.d0 1297 rcodcl(ifac,ivar,2) = rinfin 1298 rcodcl(ifac,ivar,3) = 0.d0 1299 1300 else if (ivar.eq.iu.or.ivar.eq.iv.or.ivar.eq.iw) then 1301 if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1302 itypfb(ifac) = - abs(itypfb(ifac)) 1303 if (iok.eq.0.or.iok.eq.2) iok = iok + 1 1304 else 1305 icodcl(ifac,ivar) = 13 1306 ! rcodcl(ifac,ivar,1) = User 1307 rcodcl(ifac,ivar,2) = rinfin 1308 rcodcl(ifac,ivar,3) = 0.d0 1309 endif 1310 1311 else if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1312 1313 flumbf = bmasfl(ifac) 1314 ! Outgoing flux or yplus 1315 if (flumbf.ge.-epzero.or. ivarfl(ivar).eq.f_id_yplus & 1316 .or.ivarfl(ivar).eq.f_id_z_ground) then 1317 icodcl(ifac,ivar) = 3 1318 rcodcl(ifac,ivar,1) = 0.d0 1319 rcodcl(ifac,ivar,2) = rinfin 1320 rcodcl(ifac,ivar,3) = 0.d0 1321 else 1322 itypfb(ifac) = - abs(itypfb(ifac)) 1323 if (iok.lt.2) then 1324 iok = iok + 2 1325 iivar = max(iivar, ivar) 1326 endif 1327 endif 1328 else 1329 icodcl(ifac,ivar) = 13 1330 ! rcodcl(ifac,ivar,1) = User 1331 rcodcl(ifac,ivar,2) = rinfin 1332 rcodcl(ifac,ivar,3) = 0.d0 1333 endif 1334 1335 endif 1336 enddo 1337enddo 1338 1339if (irangp.ge.0) then 1340 call parcmx(iok) 1341 call parcmx(iivar) 1342endif 1343if (iok.gt.0) then 1344 if (iok.eq.1 .or. iok.eq.3) write(nfecra,6060) 1345 if (iok.eq.2 .or. iok.eq.3) then 1346 call field_get_name(ivarfl(iivar), fname) 1347 write(nfecra,6070) trim(fname) 1348 endif 1349 call boundary_conditions_error(itypfb) 1350endif 1351 1352 1353! 6.2.a SORTIE (entree sortie libre) 1354! ================================== 1355 1356! ---> La pression a un traitement Dirichlet, les vitesses 9 ont ete 1357! traites plus haut. 1358! Le reste Dirichlet si l'utilisateur fournit une donnee 1359! (flux de masse entrant ou sortant). 1360! S'il n'y a pas de donnee utilisateur, on utilise un Neumann homogene 1361! (flux entrant et sortant) 1362 1363 1364! ---> Sortie ISOLIB 1365 1366ideb = idebty(isolib) 1367ifin = ifinty(isolib) 1368 1369do ivar = 1, nvar 1370 do ii = ideb, ifin 1371 ifac = itrifb(ii) 1372 if(icodcl(ifac,ivar).eq.0) then 1373 1374 if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1375 icodcl(ifac,ivar) = 3 1376 rcodcl(ifac,ivar,1) = 0.d0 1377 rcodcl(ifac,ivar,2) = rinfin 1378 rcodcl(ifac,ivar,3) = 0.d0 1379 else 1380 icodcl(ifac,ivar) = 1 1381 ! rcodcl(ifac,ivar,1) = Utilisateur 1382 rcodcl(ifac,ivar,2) = rinfin 1383 rcodcl(ifac,ivar,3) = 0.d0 1384 endif 1385 endif 1386 enddo 1387enddo 1388 1389! 6.2.b Free inlet (outlet) Bernoulli 1390! =================================== 1391 1392! ---> Free entrance (Bernoulli, a dirichlet is needed on the other variables 1393! than velocity and pressure, already treated) Free outlet (homogeneous 1394! Neumann) 1395 1396ideb = idebty(ifrent) 1397ifin = ifinty(ifrent) 1398 1399do ivar = 1, nvar 1400 do ii = ideb, ifin 1401 ifac = itrifb(ii) 1402 if (icodcl(ifac,ivar).eq.0) then 1403 1404 if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1405 icodcl(ifac,ivar) = 3 1406 rcodcl(ifac,ivar,1) = 0.d0 1407 rcodcl(ifac,ivar,2) = rinfin 1408 rcodcl(ifac,ivar,3) = 0.d0 1409 else 1410 icodcl(ifac,ivar) = 1 1411 ! rcodcl(ifac,ivar,1) is given by the user 1412 rcodcl(ifac,ivar,2) = rinfin 1413 rcodcl(ifac,ivar,3) = 0.d0 1414 endif 1415 endif 1416 enddo 1417enddo 1418 1419! 6.2.c Free surface 1420! ================== 1421 1422! ---> Free surface 1423 1424ideb = idebty(ifresf) 1425ifin = ifinty(ifresf) 1426 1427do ivar = 1, nvar 1428 do ii = ideb, ifin 1429 ifac = itrifb(ii) 1430 if(icodcl(ifac,ivar).eq.0) then 1431 1432 if (rcodcl(ifac,ivar,1).gt.rinfin*0.5d0) then 1433 icodcl(ifac,ivar) = 3 1434 rcodcl(ifac,ivar,1) = 0.d0 1435 rcodcl(ifac,ivar,2) = rinfin 1436 rcodcl(ifac,ivar,3) = 0.d0 1437 else 1438 icodcl(ifac,ivar) = 1 1439 ! rcodcl(ifac,ivar,1) is given by the user 1440 rcodcl(ifac,ivar,2) = rinfin 1441 rcodcl(ifac,ivar,3) = 0.d0 1442 endif 1443 endif 1444 enddo 1445enddo 1446 1447! 6.4 PAROI LISSE bis 1448! =============== 1449 1450! ---> La vitesse et les grandeurs turbulentes ont le code 5 1451! traite plus haut 1452! le reste Neumann 1453 1454ideb = idebty(iparoi) 1455ifin = ifinty(iparoi) 1456 1457do f_id = 0, nfld - 1 1458 call field_get_type(f_id, f_type) 1459 ! Is the field of type FIELD_VARIABLE? 1460 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 1461 ! Is this field not managed by CDO 1462 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 1463 call field_get_dim (f_id, f_dim) 1464 call field_get_key_int(f_id, keyvar, ivar) 1465 1466 do ii = ideb, ifin 1467 ifac = itrifb(ii) 1468 if(icodcl(ifac,ivar).eq.0) then 1469 do i_dim = 0, f_dim-1 1470 icodcl(ifac,ivar+i_dim) = 3 1471 rcodcl(ifac,ivar+i_dim,1) = 0.d0 1472 rcodcl(ifac,ivar+i_dim,2) = rinfin 1473 rcodcl(ifac,ivar+i_dim,3) = 0.d0 1474 enddo 1475 endif 1476 enddo 1477 endif 1478 endif 1479enddo 1480 1481! 6.5 PAROI RUGUEUSE bis 1482! ================== 1483 1484! ---> La vitesse et les grandeurs turbulentes ont le code 6 1485! traite plus haut 1486! le reste Neumann 1487 1488ideb = idebty(iparug) 1489ifin = ifinty(iparug) 1490 1491do f_id = 0, nfld - 1 1492 call field_get_type(f_id, f_type) 1493 ! Is the field of type FIELD_VARIABLE? 1494 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 1495 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 1496 call field_get_dim (f_id, f_dim) 1497 call field_get_key_int(f_id, keyvar, ivar) 1498 1499 do ii = ideb, ifin 1500 ifac = itrifb(ii) 1501 if(icodcl(ifac,ivar).eq.0) then 1502 do i_dim = 0, f_dim-1 1503 icodcl(ifac,ivar+i_dim) = 3 1504 rcodcl(ifac,ivar+i_dim,1) = 0.d0 1505 rcodcl(ifac,ivar+i_dim,2) = rinfin 1506 rcodcl(ifac,ivar+i_dim,3) = 0.d0 1507 enddo 1508 endif 1509 enddo 1510 endif 1511 endif 1512enddo 1513 1514!=============================================================================== 1515! Automatic treatment for some variables 1516!=============================================================================== 1517 1518! Put homogeneous Neumann on hydrostatic pressure for iphydr=1 1519! if not modified by the user 1520 1521call field_get_id_try("hydrostatic_pressure", f_id) 1522if (f_id.ge.0) then 1523 call field_get_type(f_id, f_type) 1524 ! Is the field of type FIELD_VARIABLE? 1525 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 1526 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 1527 call field_get_key_int(f_id, keyvar, ivar) 1528 do ifac = 1, nfabor 1529 if(icodcl(ifac,ivar).eq.0) then 1530 icodcl(ifac,ivar) = 3 1531 endif 1532 enddo 1533 endif 1534 endif 1535endif 1536 1537! ensure that for all variables of dimension higher than 1 and which components 1538! are coupled, icodcl is the same for all the components 1539call field_get_key_id('coupled', keycpl) 1540 1541do f_id = 0, nfld - 1 1542 call field_get_type(f_id, f_type) 1543 ! Is the field of type FIELD_VARIABLE? 1544 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 1545 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 1546 call field_get_dim (f_id, f_dim) 1547 call field_get_key_int(f_id, keycpl, icpl) 1548 if (f_dim.gt.1.and.icpl.eq.1) then 1549 call field_get_key_int(f_id, keyvar, ivar) 1550 do ifac = 1, nfabor 1551 do i_dim = 1, f_dim-1 1552 icodcl(ifac,ivar+i_dim) = icodcl(ifac,ivar) 1553 enddo 1554 enddo 1555 endif 1556 endif 1557 endif 1558enddo 1559 1560! Ensure that for all scalars without diffusion 1561! wall values ignore diffusion 1562 1563do iscal = 1, nscal 1564 ivar = isca(iscal) 1565 f_id = ivarfl(ivar) 1566 ! Name of the scalar ivar 1567 call field_get_type(f_id, f_type) 1568 if (iand(f_type, FIELD_CDO) .ne. 0) cycle 1569 call field_get_key_struct_var_cal_opt(f_id, vcopt) 1570 if (vcopt%idiff .eq. 0) then 1571 call field_get_dim (f_id, f_dim) 1572 do i_dim = 0, f_dim-1 1573 do ifac = 1, nfabor 1574 if (icodcl(ifac,ivar+i_dim).eq.5 .or. icodcl(ifac,ivar+i_dim).eq.6) then 1575 icodcl(ifac,ivar+i_dim) = 3 1576 rcodcl(ifac,ivar+i_dim,3) = 0 1577 else if (icodcl(ifac,ivar+i_dim).eq.3) then 1578 rcodcl(ifac,ivar+i_dim,3) = 0 1579 endif 1580 enddo 1581 enddo 1582 endif 1583enddo 1584 1585! Put homogeneous Neumann on wall distance 1586! if not modified by the user 1587 1588call field_get_id_try("wall_distance", f_id) 1589if (f_id.ge.0) then 1590 call field_get_type(f_id, f_type) 1591 ! Is the field of type FIELD_VARIABLE? 1592 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 1593 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 1594 call field_get_key_int(f_id, keyvar, ivar) 1595 do ifac = 1, nfabor 1596 if(icodcl(ifac,ivar).eq.0) then 1597 icodcl(ifac,ivar) = 3 1598 endif 1599 enddo 1600 endif 1601 endif 1602endif 1603 1604!=============================================================================== 1605! 7. RENFORCEMENT DIAGONALE DE LA MATRICE SI AUCUN POINTS DIRICHLET 1606!=============================================================================== 1607! On renforce si ISTAT=0 et si l'option est activee (idircl=1) 1608! Si une de ces conditions est fausse, on force ndircl a valoir 1609! au moins 1 pour ne pas decaler la diagonale. 1610 1611do ivar = 1, nvar 1612 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1613 vcopt%ndircl = 0 1614 if (vcopt%istat.gt.0 .or. vcopt%idircl.eq.0) then 1615 vcopt%ndircl = 1 1616 endif 1617 call field_set_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1618enddo 1619 1620do ivar = 1, nvar 1621 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1622 do ifac = 1, nfabor 1623 if (icodcl(ifac,ivar).eq.1 .or. icodcl(ifac,ivar).eq.5.or. icodcl(ifac,ivar).eq.15) then 1624 vcopt%ndircl = vcopt%ndircl +1 1625 endif 1626 enddo 1627 if (irangp.ge.0) call parcpt (vcopt%ndircl) 1628 call field_set_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 1629enddo 1630 1631!=============================================================================== 1632! 8. ON CALCULE LE FLUX DE MASSE AUX DIFFERENTS TYPES DE FACES 1633! ET ON IMPRIME. 1634 1635! Ca serait utile de faire l'impression dans ECRLIS, mais attention, 1636! on imprime le flux de masse du pas de temps precedent 1637! or dans ECRLIS, on imprime a la fin du pas de temps 1638! d'ou une petite incoherence possible. 1639! D'autre part, ca serait utile de sortir d'autres grandeurs 1640! (flux de chaleur par exemple, bilan de scalaires ...) 1641 1642!=============================================================================== 1643 1644call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 1645 1646! Writings: mass flux if verbosity or when log is on, and at the first two 1647! iterations and the two last iterations. Only the first iteration on navstv. 1648 1649! Always print 2 first iterations and the last 2 iterations 1650if (ntcabs - ntpabs.le.2.or.(ntcabs.ge.ntmabs-1).or.vcopt%iwarni.ge.1) then 1651 modntl = 0 1652else if (ntlist.gt.0) then 1653 modntl = mod(ntcabs,ntlist) 1654 1655! No output 1656else 1657 modntl = 1 1658endif 1659 1660if (modntl.eq.0) then 1661 1662 ! Header 1663 write(nfecra,7010) 1664 1665 ! Convective mass flux of the velocity 1666 call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 1667 call field_get_val_s(iflmab, bmasfl) 1668 1669 do ii = 1, ntypmx 1670 flumty(ii) = 0.d0 1671 enddo 1672 1673 do ii = 1, ntypmx 1674 ideb = idebty(ii) 1675 ifin = ifinty(ii) 1676 do jj = ideb, ifin 1677 ifac = itrifb(jj) 1678 flumty(ii) = flumty(ii) + bmasfl(ifac) 1679 enddo 1680 enddo 1681 1682 1683 write(nfecra,7011) 1684 1685 if (ippmod(icompf).lt.0 ) then 1686 1687 ii = ientre 1688 inb = ifinty(ii)-idebty(ii)+1 1689 if (irangp.ge.0) then 1690 call parcpt (inb) 1691 call parsom (flumty(ii)) 1692 endif 1693 write(nfecra,7020) 'Inlet ',ii,inb,flumty(ii) 1694 ii = iparoi 1695 inb = ifinty(ii)-idebty(ii)+1 1696 if (irangp.ge.0) then 1697 call parcpt (inb) 1698 call parsom (flumty(ii)) 1699 endif 1700 write(nfecra,7020) 'Smooth wall ',ii,inb,flumty(ii) 1701 ii = iparug 1702 inb = ifinty(ii)-idebty(ii)+1 1703 if (irangp.ge.0) then 1704 call parcpt (inb) 1705 call parsom (flumty(ii)) 1706 endif 1707 write(nfecra,7020) 'Rough wall ',ii,inb,flumty(ii) 1708 ii = isymet 1709 inb = ifinty(ii)-idebty(ii)+1 1710 if (irangp.ge.0) then 1711 call parcpt (inb) 1712 call parsom (flumty(ii)) 1713 endif 1714 write(nfecra,7020) 'Symmetry ',ii,inb,flumty(ii) 1715 1716 ii = isolib 1717 inb = ifinty(ii)-idebty(ii)+1 1718 if (irangp.ge.0) then 1719 call parcpt (inb) 1720 call parsom (flumty(ii)) 1721 endif 1722 write(nfecra,7020) 'Free outlet ',ii,inb,flumty(ii) 1723 1724 ii = ifrent 1725 inb = ifinty(ii)-idebty(ii)+1 1726 if (irangp.ge.0) then 1727 call parcpt (inb) 1728 call parsom (flumty(ii)) 1729 endif 1730 write(nfecra,7020) 'Free inlet ',ii,inb,flumty(ii) 1731 ii = i_convective_inlet 1732 inb = ifinty(ii)-idebty(ii)+1 1733 if (irangp.ge.0) then 1734 call parcpt (inb) 1735 call parsom (flumty(ii)) 1736 endif 1737 write(nfecra,7020) 'Convective inlet ',ii,inb,flumty(ii) 1738 1739 ii = ifresf 1740 inb = ifinty(ii)-idebty(ii)+1 1741 if (irangp.ge.0) then 1742 call parcpt (inb) 1743 call parsom (flumty(ii)) 1744 endif 1745 write(nfecra,7020) 'Free surface ',ii,inb,flumty(ii) 1746 1747 if (nbrcpl.ge.1) then 1748 if (ifaccp.eq.0) then 1749 ii = icscpl 1750 else 1751 ii = icscpd 1752 endif 1753 inb = ifinty(ii)-idebty(ii)+1 1754 if (irangp.ge.0) then 1755 call parcpt (inb) 1756 call parsom (flumty(ii)) 1757 endif 1758 write(nfecra,7020) 'Sat/Sat coupling ',ii,inb,flumty(ii) 1759 endif 1760 1761 ii = iindef 1762 inb = ifinty(ii)-idebty(ii)+1 1763 if (irangp.ge.0) then 1764 call parcpt (inb) 1765 call parsom (flumty(ii)) 1766 endif 1767 write(nfecra,7020) 'Undefined ',ii,inb,flumty(ii) 1768 1769 do ii = 1, ntypmx 1770 if ( ii.ne.ientre .and. & 1771 ii.ne.i_convective_inlet .and. & 1772 ii.ne.iparoi .and. & 1773 ii.ne.iparug .and. & 1774 ii.ne.isymet .and. & 1775 ii.ne.isolib .and. & 1776 ii.ne.ifrent .and. & 1777 ii.ne.ifresf .and. & 1778 ii.ne.icscpl .and. & 1779 ii.ne.icscpd .and. & 1780 ii.ne.iindef ) then 1781 inb = ifinty(ii)-idebty(ii)+1 1782 if (irangp.ge.0) then 1783 call parcpt (inb) 1784 call parsom (flumty(ii)) 1785 endif 1786 if(inb.gt.0) then 1787 write(nfecra,7020) 'User type ',ii,inb,flumty(ii) 1788 endif 1789 endif 1790 enddo 1791 1792 else 1793 1794 ii = ieqhcf 1795 inb = ifinty(ii)-idebty(ii)+1 1796 if (irangp.ge.0) then 1797 call parcpt (inb) 1798 call parsom (flumty(ii)) 1799 endif 1800 write(nfecra,7020) 'Sub. enth. inlet ',ii,inb,flumty(ii) 1801 1802 ii = iephcf 1803 inb = ifinty(ii)-idebty(ii)+1 1804 if (irangp.ge.0) then 1805 call parcpt (inb) 1806 call parsom (flumty(ii)) 1807 endif 1808 write(nfecra,7020) 'Ptot, Htot ',ii,inb,flumty(ii) 1809 1810 ii = iesicf 1811 inb = ifinty(ii)-idebty(ii)+1 1812 if (irangp.ge.0) then 1813 call parcpt (inb) 1814 call parsom (flumty(ii)) 1815 endif 1816 write(nfecra,7020) 'Imp inlet/outlet ',ii,inb,flumty(ii) 1817 1818 ii = isopcf 1819 inb = ifinty(ii)-idebty(ii)+1 1820 if (irangp.ge.0) then 1821 call parcpt (inb) 1822 call parsom (flumty(ii)) 1823 endif 1824 write(nfecra,7020) 'Subsonic outlet ',ii,inb,flumty(ii) 1825 1826 ii = isspcf 1827 inb = ifinty(ii)-idebty(ii)+1 1828 if (irangp.ge.0) then 1829 call parcpt (inb) 1830 call parsom (flumty(ii)) 1831 endif 1832 write(nfecra,7020) 'Supersonic outlet',ii,inb,flumty(ii) 1833 1834 ii = iparoi 1835 inb = ifinty(ii)-idebty(ii)+1 1836 if (irangp.ge.0) then 1837 call parcpt (inb) 1838 call parsom (flumty(ii)) 1839 endif 1840 write(nfecra,7020) 'Wall ',ii,inb,flumty(ii) 1841 1842 ii = isymet 1843 inb = ifinty(ii)-idebty(ii)+1 1844 if (irangp.ge.0) then 1845 call parcpt (inb) 1846 call parsom (flumty(ii)) 1847 endif 1848 write(nfecra,7020) 'Symmetry ',ii,inb,flumty(ii) 1849 1850 ii = iindef 1851 inb = ifinty(ii)-idebty(ii)+1 1852 if (irangp.ge.0) then 1853 call parcpt (inb) 1854 call parsom (flumty(ii)) 1855 endif 1856 write(nfecra,7020) 'Undefined ',ii,inb,flumty(ii) 1857 1858 do ii = 1, ntypmx 1859 if (ii.ne.iesicf .and. & 1860 ii.ne.isspcf .and. & 1861 ii.ne.ieqhcf .and. & 1862 ii.ne.iephcf .and. & 1863 ii.ne.isopcf .and. & 1864 ii.ne.iparoi .and. & 1865 ii.ne.isymet .and. & 1866 ii.ne.iindef) then 1867 inb = ifinty(ii)-idebty(ii)+1 1868 if (irangp.ge.0) then 1869 call parcpt (inb) 1870 call parsom (flumty(ii)) 1871 endif 1872 if(inb.gt.0) then 1873 write(nfecra,7020) 'User type ',ii,inb,flumty(ii) 1874 endif 1875 endif 1876 enddo 1877 1878 endif 1879 1880 write(nfecra,7030) 1881 1882endif 1883 1884!=============================================================================== 1885! FORMATS 1886!=============================================================================== 1887 1888 2020 format(/,' IFINTY : ',I10) 1889 2030 format(/,' IDEBTY : ',I10) 1890 2040 format(/,' ITYPFB : ',I10) 1891 2098 format(/, & 1892'@' ,/,& 1893'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1894'@' ,/,& 1895'@ @@ WARNING: ABORT BY BOUNDARY CONDITION CHECK' ,/,& 1896'@ ========' ,/,& 1897'@ PROBLEM WITH ORDERING OF BOUNDARY FACES' ,/,& 1898'@' ,/,& 1899'@ IFINTY(',I10 ,') = ',I10 ,/,& 1900'@ is greater than' ,/,& 1901'@ IDEBTY(',I10 ,') = ',I10 ,/,& 1902'@' ,/,& 1903'@ The calculation will not be run.' ,/,& 1904'@' ,/,& 1905'@ Contact support.' ,/,& 1906'@' ,/,& 1907'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1908'@ ',/) 1909 2099 format(/, & 1910'@' ,/,& 1911'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1912'@' ,/,& 1913'@ @@ WARNING: ABORT BY BOUNDARY CONDITION CHECK' ,/,& 1914'@ ========' ,/,& 1915'@ PROBLEM WITH ORDERING OF BOUNDARY FACES' ,/,& 1916'@ ON A DISTANT RANK.' ,/,& 1917'@' ,/,& 1918'@ IFINTY(',I10 ,') ',/,& 1919'@ is greater than' ,/,& 1920'@ IDEBTY(',I10 ,') ',/,& 1921'@' ,/,& 1922'@ The calculation will not be run.' ,/,& 1923'@' ,/,& 1924'@ Contact support.' ,/,& 1925'@' ,/,& 1926'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1927'@ ',/) 1928 1929 3099 format( & 1930'@' ,/,& 1931'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1932'@' ,/,& 1933'@ @@ WARNING: ABORT BY BOUNDARY CONDITION CHECK' ,/,& 1934'@ ========' ,/,& 1935'@ PROBLEM WITH ORDERING OF BOUNDARY FACES' ,/,& 1936'@' ,/,& 1937'@ number of faces classified by type = ',I10 ,/,& 1938'@ number of boundary faces (NFABOR) = ',I10 ,/,& 1939'@' ,/,& 1940'@ The calculation will not be run.' ,/,& 1941'@' ,/,& 1942'@ Contact support.' ,/,& 1943'@' ,/,& 1944'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1945'@ ',/) 1946 1947 1948 6010 format ( /,/, & 1949 ' ** INFORMATION ON BOUNDARY FACES TYPE',/, & 1950 ' ----------------------------------',/) 1951 6011 format ( & 1952'---------------------------------------------------------------',& 1953'----------', & 1954 /,& 1955'Boundary type Code Nb faces', & 1956 /,& 1957'---------------------------------------------------------------',& 1958'----------') 1959 6020 format ( & 1960 a17,i10,i12) 1961 6030 format( & 1962'---------------------------------------------------------------',& 1963'----------'/) 1964 1965 6060 format( & 1966'@' ,/,& 1967'@' ,/,& 1968'@' ,/,& 1969'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1970'@' ,/,& 1971'@ @@ WARNING: ABORT BY BOUNDARY CONDITION CHECK' ,/,& 1972'@ ========' ,/,& 1973'@ INCORRECT OR INCOMPLETE BOUNDARY CONDITIONS' ,/,& 1974'@' ,/,& 1975'@ At least one boundary face declared as inlet (or' ,/,& 1976'@ outlet) with prescribed velocity for which the' ,/,& 1977'@ velocity value has not been assigned for all' ,/,& 1978'@ components.' ,/,& 1979'@ The calculation will not be run. ',/,& 1980'@' ,/,& 1981'@ Verify the boundary condition definitions in the GUI' ,/,& 1982'@ or in the appropriate user subroutine.' ,/,& 1983'@' ,/,& 1984'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1985'@ ',/) 1986 1987 6070 format( & 1988'@' ,/,& 1989'@' ,/,& 1990'@' ,/,& 1991'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1992'@' ,/,& 1993'@ @@ WARNING: ABORT BY BOUNDARY CONDITION CHECK' ,/,& 1994'@ ========' ,/,& 1995'@ INCORRECT OR INCOMPLETE BOUNDARY CONDITIONS' ,/,& 1996'@' ,/,& 1997'@ At least one boundary face declared as inlet (or' ,/,& 1998'@ outlet) with prescribed velocity with an entering' ,/,& 1999'@ flow for which the value of ',a,' has not been' ,/,& 2000'@ specified (Dirichlet condition).' ,/,& 2001'@ The calculation will not be run. ',/,& 2002'@' ,/,& 2003'@ Verify the boundary condition definitions in the GUI' ,/,& 2004'@ or in the appropriate user subroutine.' ,/,& 2005'@' ,/,& 2006'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 2007'@ ',/) 2008 2009 2010 7010 format ( /,/, & 2011 ' ** BOUNDARY MASS FLOW INFORMATION',/, & 2012 ' ------------------------------',/) 2013 7011 format ( & 2014'---------------------------------------------------------------',& 2015 /,& 2016'Boundary type Code Nb faces Mass flow' , & 2017 /,& 2018'---------------------------------------------------------------') 2019 7020 format ( & 2020 a17,i10,i12,6x,e18.9) 2021 7030 format( & 2022'---------------------------------------------------------------',& 2023 /) 2024 2025 8000 format(/, & 2026'Boundary faces with free inlet/outlet detected' ,/,& 2027'Update of reference point for total pressure' ,/,& 2028' XYZP0 = ',E14.5,E14.5,E14.5 ,/) 2029 8001 format(/, & 2030'Boundary faces with pressure Dirichlet condition detected' ,/,& 2031'Update of reference point for total pressure' ,/,& 2032' XYZP0 = ',E14.5,E14.5,E14.5 ,/) 2033 2034 2035return 2036end subroutine 2037