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 tridim.f90 24!> \brief Resolution of incompressible Navier Stokes and scalar transport 25!> equations for a time step. 26!> 27!------------------------------------------------------------------------------ 28 29!------------------------------------------------------------------------------ 30! Arguments 31!------------------------------------------------------------------------------ 32! mode name role 33!------------------------------------------------------------------------------ 34!> \param[in] itrale ALE iteration number 35!> \param[in] nvar total number of variables 36!> \param[in] nscal total number of scalars 37!> \param[in] dt time step (per cell) 38!______________________________________________________________________________ 39 40subroutine tridim & 41 ( itrale , & 42 nvar , nscal , & 43 dt ) 44 45!=============================================================================== 46! Module files 47!=============================================================================== 48 49use paramx 50use numvar 51use optcal 52use entsor 53use cstphy 54use cstnum 55use pointe 56use albase 57use alstru 58use alaste 59use parall 60use period 61use ppppar 62use ppthch 63use ppincl 64use cpincl 65use coincl 66use atincl 67use ctincl 68use atsoil 69use lagran 70use radiat 71use cplsat 72use ppcpfu 73use cs_fuel_incl 74use mesh 75use field 76use rotation 77use turbomachinery 78use darcy_module 79use cs_f_interfaces 80use cs_c_bindings 81use cs_tagmr, only: rob, condb, cpb, hext, text 82use cs_tagms, only: t_metal, tmet0 83use cs_nz_tagmr 84use cs_nz_condensation 85use turbomachinery 86use cdomod 87 88! les " use pp* " ne servent que pour recuperer le pointeur IIZFPP 89 90!=============================================================================== 91 92implicit none 93 94! Arguments 95 96integer itrale 97integer nvar , nscal 98 99double precision, pointer, dimension(:) :: dt 100 101! Local variables 102 103logical must_return 104 105integer iel , ifac , ivar , iscal , iappel, n_fans 106integer iok , nfld , f_id , f_dim , f_type 107integer nbccou 108integer ntrela 109integer icmst 110integer st_id 111 112integer isvhb, iz 113integer ii 114integer iterns, inslst, icvrge 115integer italim, itrfin, itrfup, ineefl 116integer ielpdc, iflmas, iflmab 117integer kcpsyr, icpsyr 118integer key_buoyant_id, is_buoyant_fld, st_prv_id 119 120double precision cvcst 121double precision xxp0, xyp0, xzp0 122double precision relaxk, relaxe, relaxw, relaxn 123double precision hdls(6) 124 125double precision, save :: tpar, tmet 126 127integer ipass 128data ipass /0/ 129save ipass 130 131integer, pointer, dimension(:,:) :: icodcl 132integer, allocatable, dimension(:) :: isostd 133 134double precision, pointer, dimension(:) :: xprale 135double precision, pointer, dimension(:,:) :: cofale 136double precision, pointer, dimension(:,:) :: dttens 137double precision, pointer, dimension(:,:,:) :: rcodcl 138double precision, pointer, dimension(:) :: hbord, theipb 139double precision, pointer, dimension(:) :: visvdr 140double precision, allocatable, dimension(:) :: prdv2f 141double precision, allocatable, dimension(:) :: mass_source 142double precision, dimension(:), pointer :: brom, crom, cpro_rho_mass 143 144double precision, pointer, dimension(:,:) :: frcxt => null() 145double precision, pointer, dimension(:,:) :: trava 146double precision, dimension(:,:), pointer :: vel 147double precision, dimension(:,:), pointer :: cvar_vec 148double precision, dimension(:), pointer :: cvar_sca 149double precision, dimension(:), pointer :: cvar_pr 150double precision, dimension(:), pointer :: cvar_k, cvara_k, cvar_ep, cvara_ep 151double precision, dimension(:), pointer :: cvar_omg, cvara_omg 152double precision, dimension(:), pointer :: cvar_nusa, cvara_nusa 153double precision, dimension(:), pointer :: cpro_prtot 154double precision, dimension(:), pointer :: cvar_scalt, cvar_totwt 155 156! Darcy 157integer mbrom 158double precision, dimension(:), pointer :: cpro_delay, cpro_capacity, cpro_sat 159double precision, dimension(:), pointer :: cproa_delay, cproa_capacity 160double precision, dimension(:), pointer :: cproa_sat 161double precision, dimension(:), pointer :: i_mass_flux, b_mass_flux 162 163double precision, dimension(:), pointer :: coefap, cofafp, cofbfp 164double precision, dimension(:), pointer :: cpro_scal_st, cproa_scal_st 165 166type(gwf_soilwater_partition) :: sorption_scal 167 168type(var_cal_opt) :: vcopt, vcopt_u, vcopt_p 169 170!=============================================================================== 171! Interfaces 172!=============================================================================== 173 174interface 175 176 subroutine condli & 177 ( nvar , nscal , iterns , & 178 isvhb , & 179 itrale , italim , itrfin , ineefl , itrfup , & 180 cofale , xprale , & 181 icodcl , isostd , & 182 dt , rcodcl , & 183 visvdr , hbord , theipb ) 184 185 use mesh, only: nfac, nfabor 186 187 implicit none 188 189 integer nvar, nscal, iterns, isvhb 190 integer itrale , italim , itrfin , ineefl , itrfup 191 192 double precision, pointer, dimension(:) :: xprale 193 double precision, pointer, dimension(:,:) :: cofale 194 integer, pointer, dimension(:,:) :: icodcl 195 integer, dimension(nfabor+1) :: isostd 196 double precision, pointer, dimension(:) :: dt 197 double precision, pointer, dimension(:,:,:) :: rcodcl 198 double precision, pointer, dimension(:) :: visvdr, hbord, theipb 199 200 end subroutine condli 201 202 subroutine navstv & 203 ( nvar , nscal , iterns , icvrge , itrale , & 204 isostd , & 205 dt , & 206 frcxt , & 207 trava ) 208 209 use mesh, only: nfabor 210 211 implicit none 212 213 integer nvar , nscal , iterns , icvrge , itrale 214 215 integer isostd(nfabor+1) 216 217 double precision, pointer, dimension(:) :: dt 218 double precision, pointer, dimension(:,:) :: frcxt 219 double precision, pointer, dimension(:,:) :: trava 220 221 end subroutine navstv 222 223 subroutine richards(icvrge, dt) 224 225 implicit none 226 227 integer icvrge 228 double precision, pointer, dimension(:) :: dt 229 230 end subroutine richards 231 232 subroutine strdep(itrale , italim , itrfin, nvar, dt, cofale, xprale) 233 234 implicit none 235 236 integer :: itrale , italim , itrfin, nvar 237 double precision, dimension(:) :: dt 238 double precision, pointer, dimension(:) :: xprale 239 double precision, pointer, dimension(:,:) :: cofale 240 241 end subroutine strdep 242 243 subroutine cs_lagr_head_losses(n_hl_cells, cell_ids, bc_type, cku) & 244 bind(C, name='cs_lagr_head_losses') 245 use, intrinsic :: iso_c_binding 246 implicit none 247 integer(c_int), value :: n_hl_cells 248 integer(c_int), dimension(*), intent(in) :: cell_ids 249 integer(c_int), dimension(*) :: bc_type 250 real(kind=c_double), dimension(*) :: cku 251 end subroutine cs_lagr_head_losses 252 253 subroutine cs_syr_coupling_send_boundary(h_wall, t_wall) & 254 bind(C, name = 'cs_syr_coupling_send_boundary') 255 256 use, intrinsic :: iso_c_binding 257 implicit none 258 real(kind=c_double), dimension(*), intent(in) :: h_wall 259 real(kind=c_double), dimension(*), intent(inout) :: t_wall 260 261 end subroutine cs_syr_coupling_send_boundary 262 263 subroutine cs_turbulence_ke & 264 (ncesmp, icetsm, itypsm, dt, smacel, prdv2f) & 265 bind(C, name='cs_turbulence_ke') 266 use, intrinsic :: iso_c_binding 267 implicit none 268 integer(c_int), value :: ncesmp 269 integer(c_int), dimension(*), intent(in) :: icetsm, itypsm 270 real(kind=c_double), dimension(*) :: dt, smacel 271 real(kind=c_double), dimension(*), intent(in) :: prdv2f 272 end subroutine cs_turbulence_ke 273 274 subroutine cs_turbulence_kw & 275 (ncesmp, icetsm, itypsm, dt, smacel) & 276 bind(C, name='cs_turbulence_kw') 277 use, intrinsic :: iso_c_binding 278 implicit none 279 integer(c_int), value :: ncesmp 280 integer(c_int), dimension(*), intent(in) :: icetsm, itypsm 281 real(kind=c_double), dimension(*) :: dt, smacel 282 end subroutine cs_turbulence_kw 283 284 subroutine cs_turbulence_sa & 285 (ncesmp, icetsm, itypsm, dt, smacel, itypfb) & 286 bind(C, name='cs_turbulence_sa') 287 use, intrinsic :: iso_c_binding 288 implicit none 289 integer(c_int), value :: ncesmp 290 integer(c_int), dimension(*), intent(in) :: icetsm, itypsm 291 real(kind=c_double), dimension(*) :: dt, smacel 292 integer(kind=c_int), dimension(*), intent(in) :: itypfb 293 end subroutine cs_turbulence_sa 294 295 subroutine cs_turbulence_v2f & 296 (ncesmp, icetsm, itypsm, dt, smacel, prdv2f) & 297 bind(C, name='cs_turbulence_v2f') 298 use, intrinsic :: iso_c_binding 299 implicit none 300 integer(c_int), value :: ncesmp 301 integer(c_int), dimension(*), intent(in) :: icetsm, itypsm 302 real(kind=c_double), dimension(*) :: dt, smacel, prdv2f 303 end subroutine cs_turbulence_v2f 304 305 subroutine cs_volume_mass_injection_eval & 306 (nvar, ncesmp, itypsm, smacel) & 307 bind(C, name='cs_volume_mass_injection_eval') 308 use, intrinsic :: iso_c_binding 309 implicit none 310 integer(c_int), value :: nvar, ncesmp 311 integer(c_int), dimension(*), intent(in) :: itypsm 312 real(kind=c_double), dimension(*) :: smacel 313 end subroutine cs_volume_mass_injection_eval 314 315end interface 316 317!=============================================================================== 318 319! Map field arrays 320call field_get_val_v(ivarfl(iu), vel) 321 322! Number of fields 323call field_get_n_fields(nfld) 324 325call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt_u) 326call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt_p) 327 328!=============================================================================== 329! 1. Initialisation 330!=============================================================================== 331 332allocate(isostd(nfabor+1)) 333 334must_return = .false. 335 336if (vcopt_u%iwarni.ge.1) then 337 write(nfecra,1000) 338endif 339 340ipass = ipass + 1 341 342! --- Indicateur de stockage d'un scalaire et de son coef 343! d'echange associe. 344! Pour le moment, on stocke uniquement dans le cas couplage SYRTHES. 345! ISVHB donne le numero du scalaire (on suppose qu'il n'y en a qu'un). 346! Dans le cas ou on a un couplage avec le module thermique 1D en paroi, 347! on utilise le meme scalaire que celui qui sert a Syrthes (s'il y a 348! couplage Syrthes), sinon on stocke le scalaire thermique. 349 350call field_get_key_id("syrthes_coupling", kcpsyr) 351 352nbccou = cs_syr_coupling_n_couplings() 353isvhb = 0 354if (nbccou .ge. 1) then 355 do iscal = 1, nscal 356 call field_get_key_int(ivarfl(isca(iscal)), kcpsyr, icpsyr) 357 if(icpsyr.eq.1) then 358 isvhb = iscal 359 endif 360 enddo 361endif 362 363if ((nfpt1t.gt.0).and.(nbccou.le.0)) then 364 isvhb = iscalt 365endif 366 367call field_get_val_s(ivarfl(ipr), cvar_pr) 368 369if (iphydr.eq.1) then 370 call field_get_val_v_by_name('volume_forces', frcxt) 371else 372 frcxt => rvoid2 373endif 374 375! Compute z ground everywhere 376if (ippmod(iatmos).ne.-1) then 377 call cs_atmo_z_ground_compute() 378endif 379 380!=============================================================================== 381! 2. AU DEBUT DU CALCUL ON REINITIALISE LA PRESSION 382!=============================================================================== 383 384if (ippmod(idarcy).eq.-1) then 385 386! On le fait sur ntinit pas de temps, car souvent, le champ de flux de masse 387! initial n'est pas a divergence nulle (CL incluses) et l'obtention 388! d'un flux a divergence nulle coherent avec la contrainte stationnaire 389! peut prendre quelques pas de temps. 390! Remarquer que la pression est rattrapee dans l'etape de stokes. 391! On ne le fait pas dans le cas de la prise en compte de la pression 392! hydrostatique, ni dans le cas du compressible 393 394 if ( ntcabs.le.ntinit .and. isuite.eq.0 & 395 .and. (iphydr.eq.0) & 396 .and. ippmod(icompf).lt.0 & 397 .and. idilat.le.1) then 398 399 if(vcopt_p%iwarni.ge.2) then 400 write(nfecra,2000) ntcabs 401 endif 402 call field_get_val_s(iprtot, cpro_prtot) 403 xxp0 = xyzp0(1) 404 xyp0 = xyzp0(2) 405 xzp0 = xyzp0(3) 406 do iel = 1, ncel 407 cvar_pr(iel) = pred0 408 cpro_prtot(iel) = p0 + ro0*( gx*(xyzcen(1,iel)-xxp0) & 409 + gy*(xyzcen(2,iel)-xyp0) & 410 + gz*(xyzcen(3,iel)-xzp0)) 411 enddo 412 endif 413 414endif 415 416 2000 format( & 417 ' REINITIALISATION DE LA PRESSION A L''ITERATION ',I10) 418 419!=============================================================================== 420! 3. COMMUNICATIONS 421!=============================================================================== 422 423! Halo synchronization (only variables require this) 424 425if (irangp.ge.0 .or. iperio.eq.1) then 426 427 do f_id = 0, nfld-1 428 429 call field_get_dim(f_id, f_dim) 430 call field_get_type(f_id, f_type) 431 432 ! Is the field of type FIELD_VARIABLE? 433 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 434 ! Is this field not managed by CDO? 435 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 436 437 if (f_dim.eq.1) then 438 439 call field_get_val_s(f_id, cvar_sca) 440 call synsce (cvar_sca) 441 442 else if (f_dim.eq.3) then 443 444 call field_get_val_v(f_id, cvar_vec) 445 call synvie(cvar_vec) 446 447 else if (f_dim.eq.6) then 448 449 call field_get_val_v(f_id, cvar_vec) 450 call syntis(cvar_vec) 451 else 452 call csexit(1) 453 endif 454 455 endif 456 endif 457 enddo 458 459endif 460 461!=============================================================================== 462! 4. POUR IPHYDR ON DOIT COMMUNIQUER FRCXT AU PREMIER PASSAGE 463! (FRCXT SERT DANS TYPECL) 464! If icalhy=1, rho must be synchronized before copying to previous values 465!=============================================================================== 466 467if (ipass.eq.1) then 468 469! --- Communication de FRCXT 470 if (iphydr.eq.1) then 471 472 if (irangp.ge.0 .or. iperio.eq.1) then 473 call synvin (frcxt) 474 endif 475 476 endif 477 478! --- Communication de RHO 479 if (icalhy.eq.1 .or. idilat.eq.3) then 480 481 call field_get_val_s(icrom, crom) 482 if (irangp.ge.0 .or. iperio.eq.1) then 483 call synsce (crom) 484 endif 485 486 endif 487 488endif 489 490!=============================================================================== 491! 5. Temporal update of previous values (mass flux, density, ...) 492!=============================================================================== 493! We exit before SCHTMP as otherwise for 2nd order in time the value of the 494! mass flux at the previous time is overwritten by the value at the current 495! time step. When ntmabs = 0, there is no issue since all mass fluxes are 0. 496 497if (ntmabs.eq.ntpabs .and. isuite.eq.1) return 498 499! If itrale=0, we are initializing ALE, we do not touch the mas flux either. 500 501if (itrale.gt.0) then 502 iappel = 1 503 call schtmp(nscal, iappel) 504endif 505 506!=============================================================================== 507! 6. MISE A JOUR DE LA LOCALISATION DES INTERFACES DE COUPLAGE CS/CS 508!=============================================================================== 509 510! Localisation des interfaces de couplage via la librairie FVM 511 512! On fait cette mise a jour des interfaces de localisation juste apres 513! les changements de geometries dus : 514! - soit a la methode ALE (en fin de pas de temps precedent) 515! - soit a un deplacement impose (cf ci-dessus) 516 517if (nbrcpl.gt.0) call cscloc 518 519!=============================================================================== 520! 7. CALCUL DES PROPRIETES PHYSIQUES VARIABLES 521! SOIT VARIABLES AU COURS DU TEMPS 522! SOIT VARIABLES LORS D'UNE REPRISE DE CALCUL 523! (VISCOSITES ET MASSE VOLUMIQUE) 524!=============================================================================== 525 526if (vcopt_u%iwarni.ge.1) then 527 write(nfecra,1010) 528endif 529 530! Disable solid cells in fluid_solid mode 531if (fluid_solid) call cs_porous_model_set_has_disable_flag(1) 532iterns = -1 533call phyvar(nvar, nscal, iterns, dt) 534 535if (itrale.gt.0) then 536 iappel = 2 537 call schtmp(nscal, iappel) 538endif 539 540 541! REMPLISSAGE DES COEFS DE PDC 542! ON Y PASSE MEME S'IL N'Y A PAS DE PDC SUR LE PROC COURANT AU CAS OU 543! UN UTILISATEUR DECIDERAIT D'AVOIR UN COEFF DE PDC DEPENDANT DE 544! LA VITESSE MOYENNE OU MAX. 545 546if (ncpdct.gt.0) then 547 548 call cs_head_losses_compute(ckupdc) 549 550 if (iflow .eq.1) then 551 call cs_lagr_head_losses(ncepdc, icepdc, itypfb, ckupdc) 552 endif 553 554endif 555 556! Evaluate mass source term coefficients 557! (called on all ranks in case user calls global operations). 558 559if (nctsmt.gt.0) then 560 561 call cs_volume_mass_injection_eval(nvar, ncetsm, itypsm, smacel) 562 563 call cs_user_mass_source_terms(nvar, nscal, ncepdc, ncetsm, 3, & 564 icepdc, icetsm, itypsm, izctsm, & 565 dt, ckupdc, smacel) 566 567 if (ippmod(iaeros).gt.0) then 568 569 allocate(mass_source(ncelet)) 570 571 ! Cooling tower model evaporation mass exchange term 572 call cs_ctwr_bulk_mass_source_term(p0, molmass_rat, mass_source) 573 574 do ii = 1, ncetsm 575 iel = icetsm(ii) 576 smacel(ii, ipr) = smacel(ii, ipr) + mass_source(iel) 577 enddo 578 579 deallocate(mass_source) 580 endif 581 582endif 583 584!------------------------------------------------------------------------ 585!-- Fill the condensation arrays spcond for the sink term of condensation 586!-- and hpcond the thermal exchange coefficient associated to the phase 587!-- change (gas phase to liquid phase) 588!------------------------------------------------------------------------ 589 590if (nftcdt.gt.0) then 591 592 iappel = 3 593 594 ! Condensation source terms arrays initialized 595 do ii = 1, nfbpcd 596 do ivar = 1, nvar 597 itypcd(ii,ivar) = 0 598 spcond(ii,ivar) = 0.d0 599 hpcond(ii) = 0.d0 600 enddo 601 enddo 602 603 call cs_user_boundary_mass_source_terms & 604( nvar , nscal , & 605 nfbpcd , iappel , & 606 ifbpcd , itypcd , izftcd , & 607 spcond , tpar) 608 609 if (nzones.eq.1) then 610 do ii = 1, nfbpcd 611 iz = izzftcd(ii) 612 zrob (iz) = rob 613 zcondb(iz) = condb 614 zcpb (iz) = cpb 615 zhext (iz) = hext 616 ztext (iz) = text 617 ztpar (iz) = tpar 618 enddo 619 endif 620 621 ! Empiric laws used by COPAIN condensation model to 622 ! the computation of the condensation source term and 623 ! exchange coefficient of the heat transfer imposed 624 ! as boundary condition. 625 626 call condensation_copain_model & 627( nvar , nfbpcd , ifbpcd , izzftcd , & 628 spcond , hpcond ) 629 630endif 631 632!---------------------------------------------------------- 633!-- Fill the condensation arrays (svcond) for the sink term 634!-- of condensation and source term type (itypst) of each 635!-- variable solved associated to the metal structures 636!-- condensation modelling. 637!---------------------------------------------------------- 638 639if (icondv.eq.0) then 640 641 !-- Condensation source terms arrays initialized 642 do iel = 1, ncelet 643 ltmast(iel) = 0 644 do ivar = 1, nvar 645 itypst(iel, ivar) = 0 646 svcond(iel, ivar) = 0.d0 647 enddo 648 flxmst(iel) = 0.d0 649 enddo 650 651 call cs_user_metal_structures_source_terms & 652( nvar , nscal , & 653 ncmast , ltmast, & 654 itypst , izmast , & 655 svcond , tmet) 656 657 ! Condensation model to compute the sink source term 658 ! (svcond) and the heat transfer flux (flxmst) imposed 659 ! in the cells associated to the metal structures 660 ! volume where this phenomenon occurs. 661 662 call metal_structures_copain_model & 663( ncmast , ltmast , & 664 tmet , & 665 svcond(:, ipr) , flxmst ) 666 667 ! array initialization if the metal structures 668 ! condensation model is coupled with 669 ! a 0-D thermal model 670 ! FIXME add restart file later 671 if (itagms.eq.1) then 672 do icmst = 1, ncmast 673 iel = ltmast(icmst) 674 t_metal(iel,1) = tmet0 675 t_metal(iel,2) = tmet0 676 enddo 677 endif 678 679endif 680 681!=============================================================================== 682! 7.bis Current to previous for variables and GWF module 683!=============================================================================== 684 685do f_id = 0, nfld - 1 686 call field_get_type(f_id, f_type) 687 ! Is the field of type FIELD_VARIABLE? 688 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 689 ! Is this field not managed by CDO ? 690 if (iand(f_type, FIELD_CDO)/=FIELD_CDO) then 691 692 call field_current_to_previous(f_id) 693 694 ! For buoyant scalar with source termes, current to previous for them 695 call field_get_key_id("is_buoyant", key_buoyant_id) 696 call field_get_key_int(f_id, key_buoyant_id, is_buoyant_fld) 697 call field_get_key_int(f_id, kstprv, st_prv_id) 698 if (is_buoyant_fld.eq.1.and.st_prv_id.ge.0.and.itrale.gt.1) then 699 call field_get_key_int(f_id, kst, st_id) 700 call field_get_val_s(st_id, cpro_scal_st) 701 call field_get_key_int(f_id, kstprv, st_id) 702 call field_get_val_s(st_id, cproa_scal_st) 703 do iel = 1, ncel 704 cproa_scal_st(iel) = cpro_scal_st(iel) 705 enddo 706 endif 707 708 endif 709 endif 710enddo 711 712 713 714if (ippmod(idarcy).eq.1) then 715 716 ! Index of the corresponding field 717 call field_get_val_prev_s_by_name('capacity', cproa_capacity) 718 call field_get_val_prev_s_by_name('saturation', cproa_sat) 719 call field_get_val_s_by_name('capacity', cpro_capacity) 720 call field_get_val_s_by_name('saturation', cpro_sat) 721 722 do iel = 1, ncel 723 cproa_capacity(iel) = cpro_capacity(iel) 724 cproa_sat(iel) = cpro_sat(iel) 725 enddo 726 727 do ii = 1, nscal 728 ivar = ivarfl(isca(ii)) 729 call field_get_key_struct_gwf_soilwater_partition(ivarfl(isca(ii)), & 730 sorption_scal) 731 call field_get_val_s(sorption_scal%idel, cpro_delay) 732 call field_get_val_prev_s(sorption_scal%idel, cproa_delay) 733 do iel = 1, ncel 734 cproa_delay(iel) = cpro_delay(iel) 735 enddo 736 enddo 737 738endif 739 740!=============================================================================== 741! 8. Compute time step if variable 742!=============================================================================== 743 744if (vcopt_u%iwarni.ge.1) then 745 write(nfecra,1020) 746endif 747 748call dttvar & 749 ( nvar , nscal , ncepdc , ncetsm , & 750 vcopt_u%iwarni , & 751 icepdc , icetsm , itypsm , & 752 dt , & 753 ckupdc , smacel ) 754 755if (nbaste.gt.0.and.itrale.gt.nalinf) then 756 ntrela = ntcabs - ntpabs 757 call astpdt(dt) 758endif 759 760! Compute the pseudo tensorial time step if needed for the pressure solving 761 762if (iand(vcopt_p%idften, ANISOTROPIC_DIFFUSION).ne.0 & 763 .and.(ippmod(idarcy).eq.-1)) then 764 765 call field_get_val_v(idtten, dttens) 766 767 do iel = 1, ncel 768 dttens(1, iel) = dt(iel) 769 dttens(2, iel) = dt(iel) 770 dttens(3, iel) = dt(iel) 771 dttens(4, iel) = 0.d0 772 dttens(5, iel) = 0.d0 773 dttens(6, iel) = 0.d0 774 enddo 775 776 do ielpdc = 1, ncepdc 777 iel = icepdc(ielpdc) 778 779 ! dttens = (1/dt + Kpdc)^-1 780 hdls(1) = ckupdc(1, ielpdc) + 1.d0/dt(iel) 781 hdls(2) = ckupdc(2, ielpdc) + 1.d0/dt(iel) 782 hdls(3) = ckupdc(3, ielpdc) + 1.d0/dt(iel) 783 hdls(4) = ckupdc(4, ielpdc) 784 hdls(5) = ckupdc(5, ielpdc) 785 hdls(6) = ckupdc(6, ielpdc) 786 787 call symmetric_matrix_inverse(hdls, dttens(:, iel)) 788 enddo 789 790 if (irangp.ge.0) then 791 call syntis(dttens) 792 endif 793 794endif 795 796!=============================================================================== 797! RECALAGE DE LA PRESSION Pth ET MASSE VOLUMIQUE rho 798! POUR L'ALGORITHME A MASSE VOLUMIQUE VARIABLE. 799!=============================================================================== 800 801if (idilat.eq.3.or.ipthrm.eq.1) then 802 call pthrbm & 803 ( nvar , ncetsm , nfbpcd , ncmast, & 804 dt , smacel , spcond , svcond ) 805 806endif 807 808!=============================================================================== 809! 9. CHARGEMENT ET TRADUCTION DES CONDITIONS AUX LIMITES 810!=============================================================================== 811 812if(vcopt_u%iwarni.ge.1) then 813 write(nfecra,1030) 814endif 815 816! --- Methode ALE : debut de boucle d'implicitation du deplacement des 817! structures. itrfin=0 indique qu'on a besoin de refaire une iteration 818! pour Syrthes, T1D ou rayonnement. 819italim = 1 820itrfin = 1 821ineefl = 0 822if (iale.ge.1 .and. nalimx.gt.1 .and. itrale.gt.nalinf) then 823! On reserve certains tableaux pour permettre le retour a l'etat 824! initial en fin d'iteration ALE 825! - mass flux: save at the first call of schtmp 826! - conditions aux limites de gradient de P et U (car on a un appel 827! a GDRCEL pour les non orthogonalites pour calculer les CL reelles) 828! -> n'est peut-etre pas reellement necessaire 829! - la pression initiale (car RTPA est aussi ecrase dans le cas 830! ou NTERUP>1) -> on pourrait optimiser en ne reservant que si 831! necessaire ... 832 allocate(cofale(nfabor,11)) 833 allocate(xprale(ncelet)) 834 ineefl = 1 835 836 if (nbccou.gt.0 .or. nfpt1t.gt.0 .or. iirayo.gt.0) itrfin = 0 837 838else 839 cofale => null() 840 xprale => null() 841endif 842 843icodcl => null() 844rcodcl => null() 845 846300 continue 847 848hbord => null() 849theipb => null() 850visvdr => null() 851 852! --- Boucle sur navstv pour couplage vitesse/pression 853! on s'arrete a NTERUP ou quand on a converge 854! ITRFUP=0 indique qu'on a besoin de refaire une iteration 855! pour Syrthes, T1D ou rayonnement. 856itrfup = 1 857 858if (nterup.gt.1) then 859 allocate(trava(ndim,ncelet)) 860else 861 trava => rvoid2 862endif 863 864if (nterup.gt.1.or.isno2t.gt.0) then 865 if (nbccou.gt.0 .or. nfpt1t.gt.0 .or. iirayo.gt.0) itrfup = 0 866endif 867 868! Allocate temporary arrays for boundary conditions 869if (italim .eq. 1) then 870 allocate(icodcl(nfabor,nvar)) 871 allocate(rcodcl(nfabor,nvar,3)) 872endif 873if (isvhb.gt.0) then 874 allocate(hbord(nfabor)) 875endif 876! Boundary value of the thermal scalar in I' 877if (iscalt.gt.0) then 878 allocate(theipb(nfabor)) 879endif 880if (itytur.eq.4 .and. idries.eq.1) then 881 allocate(visvdr(ncelet)) 882endif 883 884icvrge = 0 885inslst = 0 886 887! Darcy : in case of a steady flow, we resolve Richards only once, 888! at the first time step. 889if (ippmod(idarcy).eq.1) then 890 if ((darcy_unsteady.eq.0).and.(ntcabs.gt.1)) goto 100 891endif 892 893iterns = 1 894do while (iterns.le.nterup) 895 896 ! Calls user BCs and computes BC coefficients 897 call condli & 898 (nvar , nscal , iterns , & 899 isvhb , & 900 itrale , italim , itrfin , ineefl , itrfup , & 901 cofale , xprale , & 902 icodcl , isostd , & 903 dt , rcodcl , & 904 visvdr , hbord , theipb ) 905 906 if (nftcdt.gt.0) then 907 ! Coefficient exchange of the enthalpy scalar 908 ivar = isca(iscalt) 909 call field_get_coefa_s(ivarfl(ivar) , coefap) 910 call field_get_coefaf_s(ivarfl(ivar), cofafp) 911 call field_get_coefbf_s(ivarfl(ivar), cofbfp) 912 913 ! Pass the heat transfer computed by the Empiric laws 914 ! of the COPAIN condensation to impose the heat transfer 915 ! at the wall due to condensation for the enthalpy scalar. 916 do ii = 1, nfbpcd 917 918 ifac= ifbpcd(ii) 919 iel = ifabor(ifac) 920 921 ! Enthalpy Boundary condition associated 922 ! to the heat transfer due to condensation. 923 cofafp(ifac) = -hpcond(ii)*coefap(ifac) 924 cofbfp(ifac) = hpcond(ii) 925 926 enddo 927 928 endif 929 930! ============================================== 931! Appel de l'interface sol-atmosphere 932! ============================================== 933 934 if (ippmod(iatmos).eq.2.and.iatsoil.eq.1.and.nfmodsol.gt.0) then 935 call field_get_val_s(icrom, crom) 936 call field_get_val_s(ivarfl(isca(iscalt)), cvar_scalt) 937 call field_get_val_s(ivarfl(isca(iymw)), cvar_totwt) 938 call solvar(cvar_scalt , cvar_totwt , & 939 crom , dt , & 940 rcodcl ) 941 endif 942 943 ! UNE FOIS LES COEFFICIENTS CALCULES, ON PEUT EN DEDUIRE PLUS 944 ! FACILEMENT (I.E. SANS RECALCULS INUTILES) LES TERMES A 945 ! ENVOYER POUR LES COUPLAGES AUX BORDS (TYPE SYRTHES) 946 947 ! En compressible et si on couple ave l'energie 948 ! on recupere le Cv de la phase couplee 949 950 if (itherm .eq. 3) then 951 952 if(icv.ge.0) then 953 cvcst = 0.d0 954 else 955 cvcst = cv0 956 endif 957 else 958 cvcst = 0.d0 959 endif 960 961 ! On envoie le tout vers SYRTHES, en distinguant CP 962 ! constant ou variable 963 if (itrfin.eq.1 .and. itrfup.eq.1) then 964 965 if (isvhb.gt.0) then 966 call cs_syr_coupling_send_boundary(hbord, theipb) 967 endif 968 969 if (iscalt.gt.0 .and. nfpt1t.gt.0) then 970 call cou1do(cvcst, hbord, theipb) 971 972 if (iirayo.ge.1) call cou1di(nfabor, iscalt, icodcl, rcodcl) 973 974 endif 975 976 ! 1-D thermal model coupling with condensation 977 ! on a surface region 978 if (nftcdt.gt.0.and.nztag1d.eq.1) then 979 call cs_tagmro & 980 ( nfbpcd , ifbpcd , izzftcd , & 981 dt ) 982 endif 983 984 ! 0-D thermal model coupling with condensation 985 ! on a volume region associated to metal structures 986 if (icondv.eq.0.and.itagms.eq.1) then 987 call cs_metal_structures_tag & 988 ( ncmast , ltmast , & 989 dt ) 990 endif 991 992 endif 993 994 ! ON N'A PLUS BESOIN DE ISVHB OU ISVHT (POUR HBORD ET TBORD) 995 ! A PARTIR D'ICI 996 997 ! Compute wall distance 998 ! TODO it has to be moved before phyvar, for that bc types have to be known 999 ! (itypfb) 1000 1001 ! (Nouvel algorithme. L'ancien est dans condli) 1002 ! In ALE, this computation is done only for the first step. 1003 if (italim.eq.1) then 1004 1005 ! Pour le moment, on suppose que l'on peut se contenter de faire 1006 ! cela au premier passage, sauf avec les maillages mobiles. Attention donc 1007 ! aux conditions aux limites variables (faces qui deviennent des parois ou 1008 ! parois qui deviennent autre chose) 1009 1010 ! Wall distance is computed if: 1011 ! - it has to be updated 1012 ! - we need it 1013 ! In case there is no wall, distance is a big value. 1014 if (imajdy.eq.0 .and. ineedy.eq.1) then 1015 1016 if (abs(icdpar).eq.1) then 1017 call distpr(itypfb) 1018 ! Deprecated algorithm 1019 else if (abs(icdpar).eq.2) then 1020 call distpr2(itypfb) 1021 endif 1022 ! Wall distance is not updated except if ALE is switched on 1023 if (iale.eq.0) imajdy = 1 1024 endif 1025 1026 endif 1027 1028 ! Compute y+ if needed 1029 ! and Van Driest "amortissement" 1030 if (itytur.eq.4 .and. idries.eq.1) then 1031 call distyp(itypfb, visvdr) 1032 endif 1033 1034!=============================================================================== 1035! 10. DANS LE CAS "zero pas de temps" EN "NON SUITE" DE CALCUL 1036! ON SORT ICI 1037!=============================================================================== 1038 1039 if (ntmabs.eq.ntpabs .and. isuite.eq.0) must_return = .true. 1040 1041 if (iilagr.eq.3) must_return = .true. 1042 1043!=============================================================================== 1044! 11. RESOLUTION DE LA VITESSE DE MAILLAGE EN ALE 1045!=============================================================================== 1046 1047 if (iale.ge.1) then 1048 1049 ! Otherwise it is done in navstv.f90 1050 if (itrale.eq.0) then 1051 1052 call cs_ale_solve_mesh_velocity(iterns, impale, ialtyb) 1053 must_return = .true. 1054 1055 endif 1056 1057 endif 1058 1059!=============================================================================== 1060! Return now if required 1061!=============================================================================== 1062 1063 if (must_return) then 1064 1065 if (associated(hbord)) deallocate(hbord) 1066 if (associated(theipb)) deallocate(theipb) 1067 if (associated(visvdr)) deallocate(visvdr) 1068 1069 if (nterup.gt.1) then 1070 deallocate(trava) 1071 endif 1072 1073 deallocate(icodcl, rcodcl) 1074 1075 deallocate(isostd) 1076 1077 return 1078 1079 endif 1080 1081!=============================================================================== 1082 1083!=============================================================================== 1084! 11. CALCUL A CHAMP DE VITESSE NON FIGE : 1085! ON RESOUT VITESSE ET TURBULENCE 1086! ON SUPPOSE QUE TOUTES LES PHASES SONT FIGEES OU AUCUNE 1087!=============================================================================== 1088 1089 ! En cas de champ de vitesse fige, on ne boucle pas sur U/P 1090 if (iccvfg.eq.0) then 1091 1092!=============================================================================== 1093! 12. Solve momentum and mass equation 1094!=============================================================================== 1095 1096 ! In case of buoyancy, scalars and momentum are coupled 1097 if (n_buoyant_scal.ge.1) then 1098 1099 if(vcopt_u%iwarni.ge.1) then 1100 write(nfecra,1060) 1101 endif 1102 1103 ! Enable solid cells in fluid_solid mode 1104 if (fluid_solid) call cs_porous_model_set_has_disable_flag(0) 1105 1106 ! Update buoyant scalar(s) 1107 call scalai(nvar, nscal , iterns , dt) 1108 1109 ! Diffusion terms for weakly compressible algorithm 1110 if (idilat.ge.4) then 1111 call diffst(nscal, iterns) 1112 endif 1113 1114 ! Update the density and turbulent viscosity 1115 !------------------------------------------- 1116 1117 ! Disable solid cells in fluid_solid mode 1118 if (fluid_solid) call cs_porous_model_set_has_disable_flag(1) 1119 call phyvar(nvar, nscal, iterns, dt) 1120 1121 ! Correct the scalar to ensure scalar conservation 1122 call field_get_key_id("is_buoyant", key_buoyant_id) 1123 call field_get_val_s(icrom,crom) 1124 call field_get_id_try("density_mass",f_id) 1125 if (f_id.ge.0) then 1126 call field_get_val_s(f_id, cpro_rho_mass) 1127 do iscal = 1, nscal 1128 ivar = isca(iscal) 1129 call field_get_key_int(ivarfl(ivar), key_buoyant_id, is_buoyant_fld) 1130 if (is_buoyant_fld.eq.1) then 1131 call field_get_val_s(ivarfl(ivar),cvar_sca) 1132 do iel = 1, ncel 1133 cvar_sca(iel) = cvar_sca(iel)*cpro_rho_mass(iel)/crom(iel) 1134 enddo 1135 endif 1136 enddo 1137 endif 1138 endif 1139 1140 if (vcopt_u%iwarni.ge.1) then 1141 write(nfecra,1040) iterns 1142 endif 1143 1144 ! Coupled solving of the velocity components 1145 1146 if (ippmod(idarcy).eq.-1) then 1147 1148 call navstv & 1149 ( nvar , nscal , iterns , icvrge , itrale , & 1150 isostd , & 1151 dt , & 1152 frcxt , & 1153 trava ) 1154 1155 ! Update local pointer arrays for transient turbomachinery computations 1156 if (iturbo.eq.2) then 1157 call field_get_val_s(ivarfl(ipr), cvar_pr) 1158 endif 1159 1160 else 1161 1162 call richards (icvrge, dt) 1163 1164 call uidapp & 1165 ( darcy_anisotropic_permeability, & 1166 darcy_anisotropic_dispersion, & 1167 darcy_gravity, & 1168 darcy_gravity_x, darcy_gravity_y, darcy_gravity_z, & 1169 darcy_unsaturated) 1170 1171 ! Darcy : update data specific to underground flow 1172 mbrom = 0 1173 call usphyv(nvar, nscal, mbrom, dt) 1174 1175 ! C version 1176 call user_physical_properties() 1177 1178 if (darcy_unsteady.eq.0) then 1179 1180 do iel = 1, ncel 1181 cproa_capacity(iel) = cpro_capacity(iel) 1182 cproa_sat(iel) = cpro_sat(iel) 1183 enddo 1184 1185 do ii = 1, nscal 1186 call field_get_key_struct_gwf_soilwater_partition(ivarfl(isca(ii)), & 1187 sorption_scal) 1188 call field_get_val_s(sorption_scal%idel, cpro_delay) 1189 call field_get_val_prev_s(sorption_scal%idel, cproa_delay) 1190 do iel = 1, ncel 1191 cproa_delay(iel) = cpro_delay(iel) 1192 enddo 1193 enddo 1194 1195 endif 1196 1197 endif 1198 1199 if (istmpf.eq.2.and.itpcol.eq.1) then 1200 iappel = 3 1201 call schtmp(nscal, iappel) 1202 endif 1203 1204 ! Si c'est la derniere iteration : INSLST = 1 1205 if ((icvrge.eq.1).or.(iterns.eq.nterup)) then 1206 1207 ! Si on a besoin de refaire une nouvelle iteration pour SYRTHES, 1208 ! rayonnement, paroi thermique 1D... 1209 ! ET que l'on est a la derniere iteration en ALE ! 1210 1211 ! ...alors, on remet a zero les indicateurs de convergence 1212 if (itrfup.eq.0.and.itrfin.eq.1) then 1213 itrfup = 1 1214 icvrge = 0 1215 iterns = iterns - 1 1216 1217 ! ...sinon, on termine 1218 else 1219 inslst = 1 1220 endif 1221 1222 ! For explicit mass flux 1223 if (istmpf.eq.0.and.inslst.eq.0) then 1224 iappel = 3 1225 call schtmp(nscal, iappel) 1226 endif 1227 1228 if (inslst.eq.1) goto 100 1229 1230 endif 1231 1232 endif ! Fin si calcul sur champ de vitesse figee 1233 1234 iterns = iterns + 1 1235 1236enddo 1237 1238100 continue 1239 1240!=============================================================================== 1241! Compute Courant and Fourier number for log 1242!=============================================================================== 1243 1244if (vcopt_u%iwarni.ge.1) then 1245 write(nfecra,1021) 1246endif 1247 1248call cs_compute_courant_fourier() 1249 1250! DARCY : the hydraulic head, identified with the pressure, 1251! has been updated by the call to Richards. 1252! As diffusion of scalars depends on hydraulic head in the 1253! general case, in order to compute the exact 1254! values of the boundary faces coefficients, we have to 1255! call boundary conditions routine again. 1256! Moreover, we need an update of the boundary 1257! conditions in the cases where they vary in time. 1258 1259if (ippmod(idarcy).eq.1) then 1260 1261 ! Calls user BCs and computes BC coefficients 1262 call condli & 1263 (nvar , nscal , iterns , & 1264 isvhb , & 1265 itrale , italim , itrfin , ineefl , itrfup , & 1266 cofale , xprale , & 1267 icodcl , isostd , & 1268 dt , rcodcl , & 1269 visvdr , hbord , theipb ) 1270 1271endif 1272 1273! Free memory 1274if (associated(hbord)) deallocate(hbord) 1275if (associated(theipb)) deallocate(theipb) 1276if (associated(visvdr)) deallocate(visvdr) 1277 1278if (nterup.gt.1) then 1279 deallocate(trava) 1280endif 1281 1282! Calcul sur champ de vitesse NON fige SUITE (a cause de la boucle U/P) 1283if (iccvfg.eq.0) then 1284 1285!=============================================================================== 1286! 13. DEPLACEMENT DES STRUCTURES EN ALE ET TEST DE BOUCLAGE IMPLICITE 1287!=============================================================================== 1288 1289 if (nbstru.gt.0.or.nbaste.gt.0) then 1290 1291 call strdep(itrale, italim, itrfin, nvar, dt, cofale, xprale) 1292 1293 ! On boucle eventuellement sur de deplacement des structures 1294 if (itrfin.ne.-1) then 1295 italim = italim + 1 1296 goto 300 1297 endif 1298 1299 ! Free memory 1300 if (associated(cofale)) then 1301 deallocate(cofale) 1302 deallocate(xprale) 1303 endif 1304 1305 endif 1306 1307 ! On ne passe dans SCHTMP que si ISTMPF.EQ.0 (explicite) 1308 if (istmpf.eq.0) then 1309 iappel = 4 1310 call schtmp(nscal, iappel) 1311 endif 1312 1313!=============================================================================== 1314! 14. RESOLUTION TURBULENCE 1315!=============================================================================== 1316 1317 iok = 0 1318 if(vcopt_u%iwarni.ge.1) then 1319 if( itytur.eq.2 .or. itytur.eq.3 & 1320 .or. itytur.eq.5 .or. iturb.eq.60 ) then 1321 iok = 1 1322 endif 1323 if(iok.eq.1) then 1324 write(nfecra,1050) 1325 endif 1326 endif 1327 1328 if ((itytur.eq.2) .or. (itytur.eq.5)) then 1329 1330 ! Reserve array to avoid recomputing production in cs_turbulence_v2f 1331 if (itytur.eq.5) then 1332 allocate(prdv2f(ncelet)) 1333 endif 1334 1335 call cs_turbulence_ke(ncetsm, icetsm, itypsm, dt, smacel, prdv2f) 1336 1337 if (itytur.eq.5) then 1338 1339 call cs_turbulence_v2f(ncetsm, icetsm, itypsm, dt, smacel, prdv2f) 1340 1341 ! Free memory 1342 deallocate(prdv2f) 1343 1344 endif 1345 1346 call field_get_val_s(ivarfl(ik), cvar_k) 1347 call field_get_val_prev_s(ivarfl(ik), cvara_k) 1348 call field_get_val_s(ivarfl(iep), cvar_ep) 1349 call field_get_val_prev_s(ivarfl(iep), cvara_ep) 1350 1351 ! RELAXATION DE K ET EPSILON SI IKECOU=0 EN INSTATIONNAIRE 1352 if (ikecou.eq.0 .and. idtvar.ge.0) then 1353 call field_get_key_struct_var_cal_opt(ivarfl(ik), vcopt) 1354 relaxk = vcopt%relaxv 1355 call field_get_key_struct_var_cal_opt(ivarfl(iep), vcopt) 1356 relaxe = vcopt%relaxv 1357 do iel = 1,ncel 1358 cvar_k(iel) = relaxk*cvar_k(iel) + (1.d0-relaxk)*cvara_k(iel) 1359 cvar_ep(iel) = relaxe*cvar_ep(iel) + (1.d0-relaxe)*cvara_ep(iel) 1360 enddo 1361 endif 1362 1363 else if(itytur.eq.3) then 1364 1365 ! Compute Alpha for EBRSM 1366 if (iturb.eq.32) then 1367 1368 call resalp(ivarfl(ial), xcl) 1369 1370 endif 1371 1372 call turrij & 1373 ( nvar , nscal , & 1374 ncepdc , ncetsm , & 1375 icepdc , icetsm , itypsm , & 1376 dt , & 1377 tslagr , & 1378 ckupdc , smacel ) 1379 1380 else if (iturb.eq.60) then 1381 1382 call cs_turbulence_kw(ncetsm, icetsm, itypsm, dt, smacel) 1383 1384 call field_get_val_s(ivarfl(ik), cvar_k) 1385 call field_get_val_prev_s(ivarfl(ik), cvara_k) 1386 call field_get_val_s(ivarfl(iomg), cvar_omg) 1387 call field_get_val_prev_s(ivarfl(iomg), cvara_omg) 1388 1389 ! RELAXATION DE K ET OMEGA SI IKECOU=0 1390 if (ikecou.eq.0 .and. idtvar.ge.0) then 1391 call field_get_key_struct_var_cal_opt(ivarfl(ik), vcopt) 1392 relaxk = vcopt%relaxv 1393 call field_get_key_struct_var_cal_opt(ivarfl(iomg), vcopt) 1394 relaxw = vcopt%relaxv 1395 do iel = 1,ncel 1396 cvar_k(iel) = relaxk*cvar_k(iel) + (1.d0-relaxk)*cvara_k(iel) 1397 cvar_omg(iel) = relaxw*cvar_omg(iel) + (1.d0-relaxw)*cvara_omg(iel) 1398 enddo 1399 end if 1400 1401 else if (iturb.eq.70) then 1402 1403 call cs_turbulence_sa(ncetsm, icetsm, itypsm, dt, smacel, itypfb) 1404 1405 call field_get_val_s(ivarfl(inusa), cvar_nusa) 1406 call field_get_val_prev_s(ivarfl(inusa), cvara_nusa) 1407 1408 ! RELAXATION DE NUSA 1409 if (idtvar.ge.0) then 1410 call field_get_key_struct_var_cal_opt(ivarfl(inusa), vcopt) 1411 relaxn = vcopt%relaxv 1412 do iel = 1,ncel 1413 cvar_nusa(iel) = relaxn*cvar_nusa(iel)+(1.d0-relaxn)*cvara_nusa(iel) 1414 enddo 1415 endif 1416 1417 endif 1418 1419endif ! Fin si calcul sur champ de vitesse fige SUITE 1420 1421! Re enable solid cells in fluid_solid mode 1422if (fluid_solid) call cs_porous_model_set_has_disable_flag(0) 1423 1424!=============================================================================== 1425! 15. RESOLUTION DES SCALAIRES 1426!=============================================================================== 1427 1428if (nscal.ge.1 .and. iirayo.gt.0) then 1429 1430 if (vcopt_u%iwarni.ge.1) then 1431 write(nfecra,1070) 1432 endif 1433 1434 ! Call the 1D radiative model 1435 ! Compute the divergence of the ir and solar radiative fluxes: 1436 if (ippmod(iatmos).ge.1.and.iatra1.ge.1) then 1437 call atr1vf() 1438 endif 1439 1440 call cs_rad_transfer_solve(itypfb, cp2fol, cp2ch, ichcor) 1441endif 1442 1443if (nscal.ge.1) then 1444 1445 if (vcopt_u%iwarni.ge.1) then 1446 write(nfecra,1060) 1447 endif 1448 1449 ! Update non-buoyant scalar(s) 1450 iterns = -1 1451 call scalai(nvar, nscal, iterns, dt) 1452 1453 ! Diffusion terms for weakly compressible algorithm 1454 if (idilat.ge.4) then 1455 call diffst(nscal, iterns) 1456 endif 1457 1458endif 1459 1460! Free memory 1461deallocate(icodcl, rcodcl) 1462 1463deallocate(isostd) 1464 1465!=============================================================================== 1466! 16. TRAITEMENT DU FLUX DE MASSE, DE LA VISCOSITE, 1467! DE LA MASSE VOLUMIQUE ET DE LA CHALEUR SPECIFIQUE POUR 1468! UN THETA SCHEMA 1469!=============================================================================== 1470 1471iappel = 5 1472call schtmp(nscal, iappel) 1473 1474!=============================================================================== 1475! Update flow through fans 1476!=============================================================================== 1477 1478n_fans = cs_fan_n_fans() 1479if (n_fans .gt. 0) then 1480 call field_get_key_int(ivarfl(iu), kimasf, iflmas) 1481 call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 1482 call field_get_val_s(iflmas, i_mass_flux) 1483 call field_get_val_s(iflmab, b_mass_flux) 1484 call field_get_val_s(icrom, crom) 1485 call field_get_val_s(ibrom, brom) 1486 call debvtl(i_mass_flux, b_mass_flux, crom, brom) 1487endif 1488 1489!=============================================================================== 1490 1491!-------- 1492! Formats 1493!-------- 1494 1495 1000 format(/, & 1496' ------------------------------------------------------------',/,& 1497 /,/,& 1498' INITIALISATIONS' ,/,& 1499' ===============' ,/) 1500 1010 format(/, & 1501' ------------------------------------------------------------',/,& 1502 /,/,& 1503' COMPUTATION OF PHYSICAL QUANTITIES' ,/,& 1504' ==================================' ,/) 1505 1020 format(/, & 1506' ------------------------------------------------------------',/,& 1507 /,/,& 1508' COMPUTATION OF CFL, FOURIER AND VARIABLE DT' ,/,& 1509' ===========================================' ,/) 1510 1511 1021 format(/, & 1512' ------------------------------------------------------------',/,& 1513 /,/,& 1514' COMPUTATION OF CFL AND FOURIER',/,& 1515' ==============================',/) 1516 1517 1030 format(/, & 1518' ------------------------------------------------------------',/,& 1519 /,/,& 1520' SETTING UP THE BOUNDARY CONDITIONS' ,/,& 1521' ==================================' ,/) 1522 1040 format(/, & 1523' ------------------------------------------------------------',/,& 1524 /,/,& 1525' SOLVING NAVIER-STOKES EQUATIONS (sub iter: ',i3,')' ,/,& 1526' ===============================' ,/) 1527 1050 format(/, & 1528' ------------------------------------------------------------',/,& 1529 /,/,& 1530' SOLVING TURBULENT VARIABLES EQUATIONS' ,/,& 1531' =====================================' ,/) 1532 1060 format(/, & 1533' ------------------------------------------------------------',/,& 1534 /,/,& 1535' SOLVING ENERGY AND SCALARS EQUATIONS' ,/,& 1536' ====================================' ,/) 1537 1070 format(/, & 1538 '------------------------------------------------------------',/,& 1539 /,/,& 1540 ' SOLVING THERMAL RADIATIVE TRANSFER' ,/,& 1541' ==================================' ,/) 1542 1543!---- 1544! End 1545!---- 1546 1547end subroutine 1548