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