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!=============================================================================== 24! Function : 25! -------- 26 27!> \file driflu.f90 28!> 29!> \brief Compute the modified convective flux for scalars with a drift. 30!> 31!------------------------------------------------------------------------------- 32 33!------------------------------------------------------------------------------- 34! Arguments 35!______________________________________________________________________________. 36! mode name role ! 37!______________________________________________________________________________! 38!> \param[in] iflid index of the current drift scalar field 39!> \param[in] dt time step (per cell) 40!> \param[in,out] imasfl scalar mass flux at interior face centers 41!> \param[in,out] bmasfl scalar mass flux at boundary face centers 42!> \param[in,out] divflu divergence of drift flux 43!______________________________________________________________________________ 44 45subroutine driflu & 46( iflid , & 47 dt , & 48 imasfl , bmasfl , & 49 divflu ) 50 51!=============================================================================== 52! Module files 53!=============================================================================== 54 55use paramx 56use dimens, only: nvar, ndimfb 57use numvar 58use entsor 59use optcal 60use cstphy 61use cstnum 62use ppppar 63use ppthch 64use pointe 65use field 66use mesh 67use parall 68use period 69use cpincl 70use cs_coal_incl 71use ppincl 72use cs_c_bindings 73 74!=============================================================================== 75 76implicit none 77 78! Arguments 79 80integer iflid 81 82double precision dt(ncelet) 83double precision imasfl(nfac), bmasfl(nfabor) 84double precision divflu(ncelet) 85 86! Local variables 87 88integer ivar 89integer ifac , iel 90integer init , inc , iccocg 91integer ifcvsl, iflmas, iflmab 92integer imrgrp, nswrgp, imligp, iwarnp 93integer iconvp, idiffp, imvisp 94integer isou , jsou 95integer f_id , f_id0 96integer iflmb0, iphydp, ivisep, itypfl 97integer keysca, iscal, keydri, iscdri, icvflb 98integer keyccl 99integer icla, jcla, jvar 100integer ivoid(1) 101integer ielup, id_x1, id_vdp_i, imasac, id_pro, id_drift 102 103double precision epsrgp, climgp, extrap 104double precision visls_0 105double precision thetap 106double precision rhovdt 107double precision omegaa 108double precision rho 109 110double precision rvoid(1) 111 112character(len=80) :: fname 113 114double precision, dimension(:), allocatable :: w1, viscce, rtrace 115double precision, dimension(:), allocatable :: coefap, coefbp 116double precision, dimension(:), allocatable :: cofafp, cofbfp 117double precision, dimension(:,:), allocatable :: coefa1 118double precision, dimension(:,:,:), allocatable :: coefb1 119double precision, dimension(:,:), allocatable :: dudt 120double precision, dimension(:), allocatable :: viscf, viscb 121double precision, dimension(:), allocatable :: flumas, flumab 122 123double precision, dimension(:), pointer :: cpro_taup 124double precision, dimension(:), pointer :: cpro_taufpt 125double precision, dimension(:,:), pointer :: cpro_drift 126double precision, dimension(:,:), pointer :: coefav, cofafv 127double precision, dimension(:,:,:), pointer :: coefbv, cofbfv 128double precision, dimension(:), pointer :: imasfl_mix, bmasfl_mix 129double precision, dimension(:,:), pointer :: vdp_i 130double precision, dimension(:), pointer :: x2 131double precision, dimension(:), pointer :: imasfl_gas, bmasfl_gas 132double precision, dimension(:), pointer :: cpro_x1, bpro_x1 133double precision, dimension(:), pointer :: brom, crom 134double precision, dimension(:,:), pointer :: vel, vela 135double precision, dimension(:), pointer :: cvar_k 136double precision, dimension(:,:), pointer :: cvar_rij 137double precision, dimension(:), pointer :: visct, cpro_viscls 138double precision, dimension(:), pointer :: cvara_var 139 140type(var_cal_opt) :: vcopt, vcopt_u 141type(var_cal_opt), target :: vcopt_loc 142type(var_cal_opt), pointer :: p_k_value 143type(c_ptr) :: c_k_value 144 145!=============================================================================== 146 147!=============================================================================== 148! 0. Key words for fields 149!=============================================================================== 150 151! Key id for scalar id 152call field_get_key_id("scalar_id", keysca) 153call field_get_key_int(iflid, keysca, iscal) 154 155ivar = isca(iscal) 156f_id0 = -1 157 158call field_get_val_prev_s(ivarfl(ivar), cvara_var) 159 160! Eventual scalar class: 161! 0: scalar belonging to no class 162! -1: deduced continuous (gas) class 163! >0: particle class 164call field_get_key_id("scalar_class", keyccl) 165call field_get_key_int(ivarfl(ivar), keyccl, icla) 166 167! Key id for drift scalar 168call field_get_key_id("drift_scalar_model", keydri) 169call field_get_key_int(iflid, keydri, iscdri) 170 171! Massic fraction of gas 172call field_get_id_try("x_c", id_x1) 173if (id_x1.ne.-1) then 174 call field_get_val_s(id_x1, cpro_x1) 175 176 ! Mass fraction of the gas at the boundary 177 call field_get_val_s_by_name("b_x_c", bpro_x1) 178 179 ! Get the mass flux of the continuous phase (gas) 180 ! that is the negative scalar class 181 do jvar = 1, nvar 182 call field_get_key_int(ivarfl(jvar), keyccl, jcla) 183 if (jcla.eq.-1) then 184 call field_get_key_int(ivarfl(jvar), kimasf, iflmas) 185 call field_get_key_int(ivarfl(jvar), kbmasf, iflmab) 186 endif 187 enddo 188 call field_get_val_s(iflmas, imasfl_gas) 189 ! Pointer to the Boundary mass flux 190 call field_get_val_s(iflmab, bmasfl_gas) 191 192endif 193 194call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 195call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt_u) 196 197! Pointers to the mass fluxes of the mix (based on mix velocity) 198call field_get_key_int(ivarfl(iu), kimasf, iflmas) 199call field_get_key_int(ivarfl(iu), kbmasf, iflmab) 200call field_get_val_s(iflmas, imasfl_mix) 201call field_get_val_s(iflmab, bmasfl_mix) 202 203! Map field arrays 204call field_get_val_v(ivarfl(iu), vel) 205call field_get_val_prev_v(ivarfl(iu), vela) 206 207!=============================================================================== 208! 1. Initialization 209!=============================================================================== 210 211! --- Physical properties 212call field_get_val_s(icrom, crom) 213call field_get_val_s(ibrom, brom) 214call field_get_val_s(ivisct, visct) 215 216if (itytur.eq.3) then 217 call field_get_val_v(ivarfl(irij), cvar_rij) 218elseif (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) then 219 call field_get_val_s(ivarfl(ik), cvar_k) 220endif 221 222! --- Brownian diffusivity 223call field_get_key_int (ivarfl(isca(iscal)), kivisl, ifcvsl) 224if (ifcvsl.ge.0) then 225 call field_get_val_s(ifcvsl, cpro_viscls) 226endif 227 228! Name of the scalar 229call field_get_name(ivarfl(ivar), fname) 230 231! Vector containing all the additional convective terms 232allocate(dudt(3, ncelet)) 233allocate(w1(ncelet), viscce(ncelet)) 234allocate(coefap(ndimfb), coefbp(ndimfb)) 235allocate(cofafp(ndimfb), cofbfp(ndimfb)) 236allocate(coefa1(3, ndimfb), coefb1(3, 3, ndimfb)) 237allocate(viscf(nfac), viscb(nfabor)) 238allocate(flumas(nfac), flumab(nfabor)) 239 240if (btest(iscdri, DRIFT_SCALAR_ADD_DRIFT_FLUX)) then 241 ! Index of the corresponding relaxation time (cpro_taup) 242 call field_get_id_try('drift_tau_'//trim(fname), id_pro) 243 if (id_pro.ne.-1) call field_get_val_s(id_pro, cpro_taup) 244 245 ! Index of the corresponding relaxation time (cpro_taup) 246 call field_get_id_try('drift_vel_'//trim(fname), id_drift) 247 if (id_drift.ne.-1) call field_get_val_v(id_drift, cpro_drift) 248 249 ! Index of the corresponding interaction time particle--eddies (cpro_taufpt) 250 if (btest(iscdri, DRIFT_SCALAR_TURBOPHORESIS)) then 251 call field_get_id('drift_turb_tau_'//trim(fname), f_id) 252 call field_get_val_s(f_id, cpro_taufpt) 253 endif 254 255 ! Initialization of the convection flux for the current particle class 256 do ifac = 1, nfac 257 viscf(ifac) = 0.d0 258 flumas(ifac) = 0.d0 259 enddo 260 261 do ifac = 1, nfabor 262 viscb(ifac) = 0.d0 263 flumab(ifac) = 0.d0 264 enddo 265 266 ! Initialization of the gas "class" convective flux by the 267 ! first particle "class": 268 ! it is initialized by the mass flux of the bulk 269 if (icla.eq.1.and.id_x1.ne.-1) then 270 do ifac = 1, nfac 271 imasfl_gas(ifac) = imasfl_mix(ifac) 272 enddo 273 274 do ifac = 1, nfabor 275 bmasfl_gas(ifac) = bmasfl_mix(ifac) 276 enddo 277 endif 278 279!=============================================================================== 280! 2. initialize the additional convective flux with the gravity term 281!=============================================================================== 282 283 ! Test if a deviation velocity of particles class exists 284 285 if (icla.ge.1) then 286 287 write(fname,'(a,i2.2)')'vd_p_' ,icla 288 call field_get_id_try(fname, id_vdp_i) 289 290 if (id_vdp_i.ne.-1) then 291 call field_get_val_v(id_vdp_i, vdp_i) 292 293 do iel = 1, ncel 294 295 rho = crom(iel) 296 cpro_drift(1, iel) = rho*vdp_i(1, iel) 297 cpro_drift(2, iel) = rho*vdp_i(2, iel) 298 cpro_drift(3, iel) = rho*vdp_i(3, iel) 299 300 enddo 301 endif 302 303 else if (icla.ge.0.and.id_pro.ne.-1.and.id_drift.ne.-1) then 304 305 do iel = 1, ncel 306 307 rho = crom(iel) 308 cpro_drift(1, iel) = rho*cpro_taup(iel)*gx 309 cpro_drift(2, iel) = rho*cpro_taup(iel)*gy 310 cpro_drift(3, iel) = rho*cpro_taup(iel)*gz 311 312 enddo 313 314 endif 315 316!=============================================================================== 317! 3. Computation of the turbophoresis and the thermophoresis terms 318!=============================================================================== 319 320 ! Initialized to 0 321 do iel = 1, ncel 322 viscce(iel) = 0.d0 323 enddo 324 325 if (btest(iscdri, DRIFT_SCALAR_TURBOPHORESIS).and.iturb.ne.0) then 326 327 ! The diagonal part is easy to implicit (Grad (K) . n = (K_j - K_i)/IJ) 328 ! Compute the K=1/3*trace(R) coefficient (diffusion of Zaichik) 329 330 if (itytur.eq.3) then 331 332 allocate(rtrace(ncel)) 333 do iel = 1, ncel 334 rtrace(iel) = cvar_rij(1,iel) + cvar_rij(2,iel) + cvar_rij(3,iel) 335 enddo 336 do iel = 1, ncel 337 ! Correction by Omega 338 omegaa = cpro_taup(iel)/cpro_taufpt(iel) 339 ! FIXME: use idifft or not? 340 viscce(iel) = 1.d0/3.d0 * cpro_taup(iel)/(1.d0+omegaa)*rtrace(iel) 341 enddo 342 deallocate(rtrace) 343 344 elseif (itytur.eq.2 .or. itytur.eq.5 .or. iturb.eq.60) then 345 346 do iel = 1, ncel 347 ! Correction by Omega 348 omegaa = cpro_taup(iel)/cpro_taufpt(iel) 349 viscce(iel) = 2.d0/3.d0*cpro_taup(iel)/(1.d0+omegaa)*cvar_k(iel) 350 enddo 351 352 endif 353 354 endif 355 356 357 if (btest(iscdri, DRIFT_SCALAR_THERMOPHORESIS)) then 358 359 ! cpro_viscls(iel): contains the Brownian motion 360 !----------------------------------------------- 361 362 if (ifcvsl.ge.0) then 363 364 do iel = 1, ncel 365 viscce(iel) = viscce(iel) + cpro_viscls(iel)/crom(iel) 366 enddo 367 368 else 369 370 call field_get_key_double(ivarfl(ivar), kvisl0, visls_0) 371 372 do iel = 1, ncel 373 viscce(iel) = viscce(iel) + visls_0/crom(iel) 374 enddo 375 376 endif 377 378 endif 379 380 if (btest(iscdri, DRIFT_SCALAR_TURBOPHORESIS).or. & 381 btest(iscdri, DRIFT_SCALAR_THERMOPHORESIS)) then 382 383 iphydp = 0 384 inc = 1 385 iccocg = 1 386 imrgrp = vcopt%imrgra 387 nswrgp = vcopt%nswrgr 388 imligp = vcopt%imligr 389 iwarnp = vcopt%iwarni 390 epsrgp = vcopt%epsrgr 391 climgp = vcopt%climgr 392 imvisp = vcopt%imvisf 393 extrap = 0 394 395 ! Face diffusivity of rho to compute rho*(Grad K . n)_face 396 do iel = 1, ncel 397 w1(iel) = crom(iel) 398 enddo 399 400 if (irangp.ge.0.or.iperio.eq.1) then 401 call synsca(w1) 402 endif 403 404 call viscfa & 405 ( imvisp , & 406 w1 , & 407 viscf , viscb ) 408 409 ! Homogeneous Neumann BC 410 do ifac = 1, nfabor 411 ! BCs for gradients 412 coefap(ifac) = 0.d0 413 coefbp(ifac) = 1.d0 414 ! BCs for fluxes 415 cofafp(ifac) = 0.d0 416 cofbfp(ifac) = 0.d0 417 enddo 418 419 init = 0 420 421 ! The computed convective flux has the dimension of rho*velocity 422 call itrmas & 423 ( f_id0 , init , inc , imrgrp , iccocg , nswrgp , imligp , iphydp , & 424 0 , iwarnp , & 425 epsrgp , climgp , extrap , & 426 rvoid , & 427 viscce , & 428 coefap , coefbp , & 429 cofafp , cofbfp , & 430 viscf , viscb , & 431 w1 , & 432 flumas , flumab ) 433 434 ! TODO add extradiagonal part 435 endif 436 437!=============================================================================== 438! 4. Centrifugal force (particular derivative Du/Dt) 439!=============================================================================== 440 441 if (btest(iscdri, DRIFT_SCALAR_CENTRIFUGALFORCE)) then 442 443 do iel = 1, ncel 444 445 rhovdt = crom(iel)*volume(iel)/dt(iel) 446 447 dudt(1,iel) = -rhovdt*(vel(1,iel)-vela(1,iel)) 448 dudt(2,iel) = -rhovdt*(vel(2,iel)-vela(2,iel)) 449 dudt(3,iel) = -rhovdt*(vel(3,iel)-vela(3,iel)) 450 enddo 451 452 iconvp = 1 453 idiffp = 0 454 inc = 1 455 ivisep = 0 456 icvflb = 0 457 458 ! Reset viscf and viscb 459 do ifac = 1, nfac 460 viscf(ifac) = 0.d0 461 enddo 462 463 do ifac = 1, nfabor 464 viscb(ifac) = 0.d0 465 enddo 466 467 ! Get Boundary conditions of the velocity 468 call field_get_coefa_v (ivarfl(iu), coefav) 469 call field_get_coefb_v (ivarfl(iu), coefbv) 470 call field_get_coefaf_v(ivarfl(iu), cofafv) 471 call field_get_coefbf_v(ivarfl(iu), cofbfv) 472 473 ! The added convective scalar mass flux is: 474 ! (thetap*Y_\face-imasac*Y_\celli)*mf. 475 ! When building the implicit part of the rhs, one 476 ! has to impose 1 on mass accumulation. 477 imasac = 1 478 479 vcopt_loc = vcopt_u 480 481 vcopt_loc%iconv = iconvp 482 vcopt_loc%istat = -1 483 vcopt_loc%idiff = idiffp 484 vcopt_loc%idifft = -1 485 vcopt_loc%iswdyn = -1 486 vcopt_loc%nswrsm = -1 487 vcopt_loc%iwgrec = 0 488 vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 489 vcopt_loc%epsilo = -1 490 vcopt_loc%epsrsm = -1 491 492 p_k_value => vcopt_loc 493 c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 494 495 ! Warning: cs_balance_vector adds "-( grad(u) . rho u)" 496 call cs_balance_vector & 497 ( idtvar , ivarfl(iu) , imasac , inc , ivisep , & 498 c_k_value , vel , vel , coefav , coefbv , cofafv , & 499 cofbfv , imasfl_mix , bmasfl_mix , viscf , viscb , rvoid , & 500 rvoid , rvoid , rvoid , rvoid , & 501 icvflb , ivoid , dudt ) 502 503 do iel = 1, ncel 504 cpro_drift(1, iel) = cpro_drift(1, iel) & 505 + cpro_taup(iel)*dudt(1, iel)/volume(iel) 506 cpro_drift(2, iel) = cpro_drift(2, iel) & 507 + cpro_taup(iel)*dudt(2, iel)/volume(iel) 508 cpro_drift(3, iel) = cpro_drift(3, iel) & 509 + cpro_taup(iel)*dudt(3, iel)/volume(iel) 510 enddo 511 512 endif 513 514!=============================================================================== 515! 5. Electrophoresis term 516!=============================================================================== 517 518 if (btest(iscdri, DRIFT_SCALAR_ELECTROPHORESIS)) then 519 520 !TODO 521 call csexit(1) 522 523 endif 524 525!=============================================================================== 526! 6. Finalization of the mass flux of the current class 527!=============================================================================== 528 529 ! For all scalar with a drift excpted the gas phase which is deduced 530 ! And for those whom mass flux is imposed elsewhere 531 if (icla.ge.0.and.(.not.btest(iscdri, DRIFT_SCALAR_IMPOSED_MASS_FLUX))) then 532 533 ! Homogeneous Neumann at the boundary 534 if (btest(iscdri, DRIFT_SCALAR_ZERO_BNDY_FLUX)) then 535 do ifac = 1, nfabor 536 do isou = 1, 3 537 coefa1(isou, ifac) = 0.d0 538 do jsou = 1, 3 539 coefb1(isou, jsou, ifac) = 0.d0 540 enddo 541 enddo 542 enddo 543 else if (btest(iscdri, DRIFT_SCALAR_ZERO_BNDY_FLUX_AT_WALLS)) then 544 do ifac = 1, nfabor 545 do isou = 1, 3 546 coefa1(isou, ifac) = 0.d0 547 do jsou = 1, 3 548 coefb1(isou, jsou, ifac) = 0.d0 549 enddo 550 if (itypfb(ifac).ne.iparoi .and. itypfb(ifac).ne.iparug) then 551 coefb1(isou, isou, ifac) = 1.d0 552 endif 553 enddo 554 enddo 555 else 556 do ifac = 1, nfabor 557 do isou = 1, 3 558 coefa1(isou, ifac) = 0.d0 559 do jsou = 1, 3 560 coefb1(isou, jsou, ifac) = 0.d0 561 enddo 562 coefb1(isou, isou, ifac) = 1.d0 563 enddo 564 enddo 565 endif 566 567 init = 0 568 inc = 1 569 iflmb0 = 0 570 itypfl = 0 ! drift has already been multiplied by rho 571 imrgrp = vcopt%imrgra 572 nswrgp = vcopt%nswrgr 573 imligp = vcopt%imligr 574 iwarnp = vcopt%iwarni 575 epsrgp = vcopt%epsrgr 576 climgp = vcopt%climgr 577 extrap = 0 578 579 call inimav & 580 ( f_id0 , itypfl , & 581 iflmb0 , init , inc , imrgrp , nswrgp , imligp , & 582 iwarnp , & 583 epsrgp , climgp , & 584 crom , brom , & 585 cpro_drift , & 586 coefa1 , coefb1 , & 587 flumas , flumab ) 588 589 ! Update the convective flux, exception for the Gas "class" 590 do ifac = 1, nfac 591 imasfl(ifac) = imasfl_mix(ifac) + flumas(ifac) 592 enddo 593 594 do ifac = 1, nfabor 595 bmasfl(ifac) = bmasfl_mix(ifac) + flumab(ifac) 596 enddo 597 598 599!=============================================================================== 600! 7. Deduce the convective flux of the gas "class" by removing the flux 601! of the current particle "class": 602! (rho x1 V1)_ij = (rho Vs)_ij - sum_classes (rho x2 V2)_ij 603!=============================================================================== 604 605 if (icla.ge.1) then 606 607 write(fname,'(a,i2.2)')'x_p_' ,icla 608 call field_get_id_try(fname, f_id) 609 610 if (f_id.ne.-1) then 611 call field_get_val_s(f_id, x2) 612 613 do ifac = 1, nfac 614 615 ! Upwind value of x2 at the face, consistent with the 616 ! other transport equations 617 ielup = ifacel(2,ifac) 618 if (imasfl(ifac).ge.0.d0) ielup = ifacel(1,ifac) 619 620 imasfl_gas(ifac) = imasfl_gas(ifac) - x2(ielup)*imasfl(ifac) 621 enddo 622 623 do ifac = 1, nfabor 624 ! TODO Upwind value of x2 at the face, consistent with the 625 ! other transport equations 626 !if (bmasfl(ifac).ge.0.d0) 627 ielup = ifabor(ifac) 628 bmasfl_gas(ifac) = bmasfl_gas(ifac) - x2(ielup)*bmasfl(ifac) 629 enddo 630 endif 631 endif 632 633 ! Finalize the convective flux of the gas "class" by scaling by x1 634 ! (rho x1 V1)_ij = (rho Vs)_ij - sum_classes (rho x2 V2)_ij 635 ! Warning, x1 at the face must be computed so that it is consistent 636 ! with an upwind scheme on (rho V1) 637 else if (icla.eq.-1.and.id_x1.ne.-1) then 638 639 do ifac = 1, nfac 640 641 ! Upwind value of x2 at the face, consistent with the 642 ! other transport equations 643 ielup = ifacel(2,ifac) 644 if (imasfl_gas(ifac).ge.0.d0) ielup = ifacel(1,ifac) 645 646 imasfl_gas(ifac) = imasfl_gas(ifac)/cpro_x1(ielup) 647 enddo 648 649 do ifac = 1, nfabor 650 ! Upwind value of x1 at the face, consistent with the 651 ! other transport equations 652 ielup = ifabor(ifac) 653 if (bmasfl_gas(ifac).lt.0.d0) then 654 bmasfl_gas(ifac) = bmasfl_gas(ifac)/bpro_x1(ifac) 655 else 656 bmasfl_gas(ifac) = bmasfl_gas(ifac)/cpro_x1(ielup) 657 endif 658 enddo 659 endif 660 661endif 662 663!=============================================================================== 664! 8. Mass aggregation term of the additional part "div(rho(u_p-u_f))" 665!=============================================================================== 666 667init = 1 668iconvp = vcopt%iconv 669thetap = vcopt%thetav!FIXME not multiplied by theta? 670 671! recompute the difference between mixture and the class 672if (btest(iscdri, DRIFT_SCALAR_IMPOSED_MASS_FLUX)) then 673 do ifac = 1, nfac 674 flumas(ifac) = - imasfl_mix(ifac) 675 enddo 676 677 do ifac = 1, nfabor 678 flumab(ifac) = - bmasfl_mix(ifac) 679 enddo 680else 681 do ifac = 1, nfac 682 flumas(ifac) = imasfl(ifac) - imasfl_mix(ifac) 683 enddo 684 685 do ifac = 1, nfabor 686 flumab(ifac) = bmasfl(ifac) - bmasfl_mix(ifac) 687 enddo 688endif 689 690call divmas(init, flumas, flumab, divflu) 691 692! Free memory 693deallocate(viscce) 694deallocate(dudt) 695deallocate(w1) 696deallocate(viscf, viscb) 697deallocate(flumas, flumab) 698deallocate(coefap, coefbp) 699deallocate(cofafp, cofbfp) 700deallocate(coefa1, coefb1) 701 702!-------- 703! Formats 704!-------- 705 706!---- 707! End 708!---- 709 710return 711end subroutine 712