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 covofv.f90
28!>
29!> \brief This subroutine performs the solving the convection/diffusion
30!> equation (with eventually source terms and/or drift) for a vectorial quantity
31!> over a time step.
32!>
33!-------------------------------------------------------------------------------
34
35!-------------------------------------------------------------------------------
36! Arguments
37!______________________________________________________________________________.
38!  mode           name          role                                           !
39!______________________________________________________________________________!
40!> \param[in]     nvar          total number of variables
41!> \param[in]     nscal         total number of scalars
42!> \param[in]     ncepdp        number of cells with head loss
43!> \param[in]     ncesmp        number of cells with mass source term
44!> \param[in]     iterns        Navier-Stokes iteration number
45!> \param[in]     iscal         scalar number
46!> \param[in]     icepdc        index of cells with head loss
47!> \param[in]     icetsm        index of cells with mass source term
48!> \param[in]     itypsm        type of mass source term for the variables
49!> \param[in]     dt            time step (per cell)
50!> \param[in]     ckupdc        work array for the head loss
51!> \param[in]     smacel        variable value associated to the mass source
52!>                               term (for ivar=ipr, smacel is the mass flux
53!>                               \f$ \Gamma^n \f$)
54!> \param[in]     viscf         visc*surface/dist at internal faces
55!> \param[in]     viscb         visc*surface/dist at boundary faces
56!_______________________________________________________________________________
57
58subroutine covofv &
59 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
60   iterns , iscal  ,                                              &
61   icepdc , icetsm ,                                              &
62   itypsm ,                                                       &
63   dt     ,                                                       &
64   ckupdc , smacel ,                                              &
65   viscf  , viscb  )
66
67!===============================================================================
68
69!===============================================================================
70! Module files
71!===============================================================================
72
73use, intrinsic :: iso_c_binding
74
75use paramx
76use numvar
77use entsor
78use optcal
79use cstphy
80use cstnum
81use ppppar
82use ppthch
83use coincl
84use cpincl
85use cs_fuel_incl
86use ppincl
87use ppcpfu
88use lagran
89use radiat
90use field
91use field_operator
92use mesh
93use parall
94use period
95use cs_f_interfaces
96use atchem
97use darcy_module
98use cs_c_bindings
99
100!===============================================================================
101
102implicit none
103
104! Arguments
105
106integer          nvar   , nscal
107integer          ncepdp , ncesmp
108integer          iterns , iscal
109
110integer          icepdc(ncepdp)
111integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
112
113double precision dt(ncelet)
114double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar)
115double precision viscf(nfac), viscb(nfabor)
116
117! Local variables
118
119logical          lprev
120character(len=80) :: chaine
121integer          ivar
122integer          ifac , iel, isou, jsou
123integer          iiun, ibcl
124integer          ivarsc
125integer          iiscav
126integer          ifcvsl, iflmas, iflmab
127integer          iwarnp
128integer          iconvp, imvisp
129integer          iescap, ivissv
130integer          idftnp
131integer          iflid , st_prv_id, st_id,  keydri, iscdri
132integer          icvflb, f_dim, iflwgr
133integer          key_buoyant_id, is_buoyant_fld
134
135integer          ivoid(1)
136
137double precision sclnor
138double precision thetv , thets , thetap, thetp1
139double precision smbexp(3)
140double precision temp, idifftp
141double precision turb_schmidt, visls_0
142
143double precision rvoid(1)
144
145double precision, allocatable, dimension(:) :: w1
146double precision, allocatable, dimension(:,:) :: smbrv, gavinj
147double precision, allocatable, dimension(:,:,:) :: fimp
148double precision, allocatable, dimension(:,:) :: viscce
149double precision, allocatable, dimension(:,:) :: weighf
150double precision, allocatable, dimension(:) :: weighb
151
152double precision, dimension(:,:), pointer :: visten
153double precision, dimension(:,:), pointer :: cpro_wgrec_v
154double precision, dimension(:), pointer :: cpro_wgrec_s
155double precision, dimension(:), pointer :: imasfl, bmasfl
156double precision, dimension(:), pointer :: crom, croma, pcrom
157double precision, dimension(:,:), pointer :: coefap, cofafp
158double precision, dimension(:,:,:), pointer :: coefbp, cofbfp
159double precision, dimension(:), pointer :: visct
160double precision, dimension(:,:), pointer :: cpro_vect_st, cproa_vect_st
161double precision, dimension(:), pointer :: cpro_viscls
162! Darcy arrays
163double precision, allocatable, dimension(:) :: diverg, divflu
164double precision, dimension(:), pointer :: cpro_delay, cpro_sat
165double precision, dimension(:), pointer :: cproa_delay, cproa_sat
166double precision, dimension(:,:), pointer :: cvar_var, cvara_var
167
168type(var_cal_opt) :: vcopt
169type(var_cal_opt), target :: vcopt_loc
170type(var_cal_opt), pointer :: p_k_value
171type(c_ptr) :: c_k_value
172
173type(gwf_soilwater_partition) :: sorption_scal
174
175!===============================================================================
176
177!===============================================================================
178! 1. Initialization
179!===============================================================================
180
181! --- Variable number
182ivar   = isca(iscal)
183! Index of the field
184iflid = ivarfl(ivar)
185
186call field_get_val_v(iflid, cvar_var)
187call field_get_val_prev_v(iflid, cvara_var)
188
189! Key id for buoyant field (inside the Navier Stokes loop)
190call field_get_key_id("is_buoyant", key_buoyant_id)
191call field_get_key_int(iflid, key_buoyant_id, is_buoyant_fld)
192
193! If the vector is buoyant, it is inside the Navier Stokes loop, and so iterns >=1
194! otherwise it is outside of the loop and iterns = -1.
195if (  (is_buoyant_fld.eq. 1 .and. iterns.eq.-1) &
196  .or.(is_buoyant_fld.eq. 0 .and. iterns.ne.-1)) return
197
198! Key id for drift scalar
199call field_get_key_id("drift_scalar_model", keydri)
200
201! Id of the mass flux
202call field_get_key_int(iflid, kimasf, iflmas) ! interior mass flux
203! Pointer to the internal mass flux
204call field_get_val_s(iflmas, imasfl)
205
206! Id of the mass flux
207call field_get_key_int(iflid, kbmasf, iflmab) ! boundary mass flux
208! Pointer to the Boundary mass flux
209call field_get_val_s(iflmab, bmasfl)
210
211call field_get_key_struct_var_cal_opt(iflid, vcopt)
212
213if (vcopt%iwgrec.eq.1) then
214  ! Id weighting field for gradient
215  call field_get_key_int(iflid, kwgrec, iflwgr)
216  call field_get_dim(iflwgr, f_dim)
217  if (f_dim.gt.1) then
218    call field_get_val_v(iflwgr, cpro_wgrec_v)
219  else
220    call field_get_val_s(iflwgr, cpro_wgrec_s)
221  endif
222endif
223
224! Allocate temporary arrays
225allocate(smbrv(3,ncelet), fimp(3,3,ncelet))
226allocate(weighf(2,nfac))
227allocate(weighb(nfabor))
228
229if (ippmod(idarcy).eq.1) then
230  allocate(diverg(ncelet))
231endif
232
233! --- Numero du scalaire eventuel associe dans le cas fluctuation
234!         et numero de variable de calcul
235iiscav = iscavr(iscal)
236if (iiscav.gt.0.and.iiscav.le.nscal) then
237  ivarsc = isca(iiscav)
238else
239  ivarsc = 0
240endif
241
242! --- Numero des grandeurs physiques
243call field_get_val_s(icrom, crom)
244call field_have_previous(icrom, lprev)
245if (lprev) then
246  call field_get_val_prev_s(icrom, croma)
247endif
248call field_get_val_s(ivisct, visct)
249
250call field_get_key_int (ivarfl(isca(iscal)), kivisl, ifcvsl)
251if (ifcvsl.ge.0) then
252  call field_get_val_s(ifcvsl, cpro_viscls)
253endif
254
255! --- Numero de propriété du terme source si extrapolation
256call field_get_key_int(iflid, kstprv, st_prv_id)
257if (st_prv_id .ge.0) then
258  call field_get_val_v(st_prv_id, cproa_vect_st)
259else
260  cproa_vect_st => null()
261endif
262
263! S pour Source, V pour Variable
264thets  = thetss(iscal)
265thetv  = vcopt%thetav
266
267call field_get_name(iflid, chaine)
268
269if (vcopt%iwarni.ge.1) then
270  write(nfecra,1000) chaine(1:16)
271endif
272
273! Retrieve turbulent Schmidt value for current vector
274call field_get_key_double(ivarfl(isca(iscal)), ksigmas, turb_schmidt)
275
276! Reference diffusivity
277call field_get_key_double(ivarfl(isca(iscal)), kvisl0, visls_0)
278
279!===============================================================================
280! 2. Source terms
281!===============================================================================
282
283! --> Initialization
284
285do iel = 1, ncel
286  do isou = 1, 3
287    do jsou = 1, 3
288      fimp(isou,jsou,iel) = 0.d0
289    enddo
290    smbrv(isou,iel) = 0.d0
291  enddo
292enddo
293
294call ustsvv &
295( nvar   , nscal  , ncepdp , ncesmp ,                            &
296  iscal  ,                                                       &
297  icepdc , icetsm , itypsm ,                                     &
298  dt     ,                                                       &
299  ckupdc , smacel , smbrv  , fimp )
300
301! C version
302call user_source_terms(ivarfl(isca(iscal)), smbrv, fimp)
303
304! Store the source terms for convective limiter
305call field_get_key_int(iflid, kst, st_id)
306if (st_id .ge.0) then
307  call field_get_dim(st_id, f_dim)
308  if (f_dim.ne.3) then
309    call csexit(1)
310  endif
311  call field_get_val_v(st_id, cpro_vect_st)
312
313  do iel = 1, ncel
314    !Fill the scalar source term field
315    do isou = 1, 3
316      cpro_vect_st(isou,iel) = smbrv(isou,iel)
317    enddo
318  end do
319  ! Handle parallelism and periodicity
320  if (irangp.ge.0.or.iperio.eq.1) then
321    call synvin(cpro_vect_st)
322  endif
323end if
324
325! Si on extrapole les TS :
326!   SMBRV recoit -theta TS du pas de temps precedent
327!     (on aurait pu le faire avant ustsvv, mais avec le risque que
328!      l'utilisateur l'ecrase)
329!   SMBRV recoit la partie du terme source qui depend de la variable
330!   A l'ordre 2, on suppose que le ROVSDT fourni par l'utilisateur est <0
331!     on implicite le terme (donc ROVSDT*RTPA va dans SMBRV)
332!   En std, on adapte le traitement au signe de ROVSDT, mais ROVSDT*RTPA va
333!     quand meme dans SMBRV (pas d'autre choix)
334if (st_prv_id .ge. 0) then
335  do iel = 1, ncel
336    ! Stockage temporaire pour economiser un tableau
337    do isou = 1, 3
338      smbexp(isou) = cproa_vect_st(isou,iel)
339
340      ! Terme source utilisateur explicite
341      cproa_vect_st(isou,iel) = smbrv(isou,iel)
342      ! Terme source du pas de temps precedent et
343      ! On suppose -fimp(isou,isou) > 0 : on implicite
344      !    le terme source utilisateur (le reste)
345      smbrv(isou,iel) = fimp(isou,isou,iel)*cvara_var(isou,iel)               &
346                      - thets*smbexp(isou)
347      ! Diagonale
348      fimp(isou,isou,iel) = - thetv*fimp(isou,isou,iel)
349    enddo
350  enddo
351
352! Si on n'extrapole pas les TS :
353else
354  do iel = 1, ncel
355    do isou = 1, 3
356      ! Terme source utilisateur
357      smbrv(isou,iel) = smbrv(isou,iel) &
358                      + fimp(isou,isou,iel)*cvara_var(isou,iel)
359      ! Diagonale
360      fimp(isou,isou,iel) = max(-fimp(isou,isou,iel),zero)
361    enddo
362  enddo
363endif
364
365! --> Physique particulieres
366!     Ordre 2 non pris en compte
367
368if (ippmod(iphpar).ge.1) then
369  call pptsvv(iscal, smbrv)
370endif
371
372! Mass source term
373if (ncesmp.gt.0) then
374
375  ! Entier egal a 1 (pour navsto : nb de sur-iter)
376  iiun = 1
377
378  ! On incremente SMBRV par -Gamma RTPA et ROVSDT par Gamma
379  allocate(gavinj(3,ncelet))
380  call catsmv(ncesmp, 1, icetsm, itypsm(:,ivar),                     &
381              cell_f_vol, cvara_var, smacel(:,ivar), smacel(:,ipr),  &
382              smbrv, fimp, gavinj)
383
384  ! Si on extrapole les TS on met Gamma Pinj dans cproa_vect_st
385  if (st_prv_id .ge. 0) then
386    do iel = 1, ncel
387      do isou = 1, 3
388        cproa_vect_st(isou,iel) = cproa_vect_st(isou,iel) + gavinj(isou,iel)
389      enddo
390    enddo
391  ! Sinon on le met directement dans SMBRV
392  else
393    do iel = 1, ncel
394      do isou = 1, 3
395        smbrv(isou,iel) = smbrv(isou,iel) + gavinj(isou,iel)
396      enddo
397    enddo
398  endif
399
400endif
401
402if (st_prv_id .ge. 0) then
403  thetp1 = 1.d0 + thets
404  do iel = 1, ncel
405    do isou = 1, 3
406      smbrv(isou,iel) = smbrv(isou,iel) + thetp1 * cproa_vect_st(isou,iel)
407    enddo
408  enddo
409endif
410
411! Compressible algorithm
412! or Low Mach compressible algos with mass flux prediction
413if (ippmod(icompf).ge.0.or.(idilat.gt.1.and.ipredfl.eq.1.and.irovar.eq.1)) then
414  pcrom => croma
415
416! Low Mach compressible algos (conservative in time).
417! Same algo for Volume of Fluid method.
418else if ((idilat.gt.1.or.ivofmt.gt.0).and.irovar.eq.1) then
419  if (iterns.eq.1) then
420    call field_get_val_prev2_s(icrom, pcrom)
421  else
422    call field_get_val_prev_s(icrom, pcrom)
423  endif
424
425! Deprecated algo or constant density
426else
427  call field_get_val_s(icrom, pcrom)
428endif
429
430! Get the the order of the diffusivity tensor of the variable
431call field_get_key_int_by_name(iflid, "diffusivity_tensor", idftnp)
432
433! "VITESSE" DE DIFFUSION FACETTE
434
435! On prend le MAX(mu_t,0) car en LES dynamique mu_t peut etre negatif
436! (clipping sur (mu + mu_t)). On aurait pu prendre
437! MAX(K + K_t,0) mais cela autoriserait des K_t negatif, ce qui est
438! considere ici comme non physique.
439if (vcopt%idiff.ge.1) then
440  ! Scalar diffusivity
441  if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then
442
443    allocate(w1(ncelet))
444    idifftp = vcopt%idifft
445    if (ifcvsl.lt.0) then
446      do iel = 1, ncel
447        w1(iel) = visls_0                                         &
448           + idifftp*max(visct(iel),zero)/turb_schmidt
449      enddo
450    else
451      do iel = 1, ncel
452        w1(iel) = cpro_viscls(iel)                                &
453           + idifftp*max(visct(iel),zero)/turb_schmidt
454      enddo
455    endif
456
457    if (vcopt%iwgrec.eq.1) then
458      ! Weighting for gradient
459      do iel = 1, ncel
460        cpro_wgrec_s(iel) = w1(iel)
461      enddo
462      call synsca(cpro_wgrec_s)
463    endif
464
465    imvisp = vcopt%imvisf
466
467    call viscfa &
468   ( imvisp ,                      &
469     w1     ,                      &
470     viscf  , viscb  )
471
472    deallocate(w1)
473
474  ! Symmetric tensor diffusivity (GGDH)
475  elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then
476
477    ! Allocate temporary arrays
478    allocate(viscce(6,ncelet))
479
480    if (iturb.ne.32) then
481      call field_get_val_v(ivsten, visten)
482    else ! EBRSM and (GGDH or AFM)
483      call field_get_val_v(ivstes, visten)
484    endif
485
486    if (ifcvsl.lt.0) then
487      do iel = 1, ncel
488
489        temp = vcopt%idifft*ctheta(iscal)/csrij
490        viscce(1,iel) = temp*visten(1,iel) + visls_0
491        viscce(2,iel) = temp*visten(2,iel) + visls_0
492        viscce(3,iel) = temp*visten(3,iel) + visls_0
493        viscce(4,iel) = temp*visten(4,iel)
494        viscce(5,iel) = temp*visten(5,iel)
495        viscce(6,iel) = temp*visten(6,iel)
496
497      enddo
498    else
499      do iel = 1, ncel
500
501        temp = vcopt%idifft*ctheta(iscal)/csrij
502        viscce(1,iel) = temp*visten(1,iel) + cpro_viscls(iel)
503        viscce(2,iel) = temp*visten(2,iel) + cpro_viscls(iel)
504        viscce(3,iel) = temp*visten(3,iel) + cpro_viscls(iel)
505        viscce(4,iel) = temp*visten(4,iel)
506        viscce(5,iel) = temp*visten(5,iel)
507        viscce(6,iel) = temp*visten(6,iel)
508
509      enddo
510    endif
511
512    if (vcopt%iwgrec.eq.1) then
513      ! Weighting for gradient
514      do iel = 1, ncel
515        do isou = 1, 6
516          cpro_wgrec_v(isou,iel) = viscce(isou,iel)
517        enddo
518      enddo
519      call syntis(cpro_wgrec_v)
520    endif
521
522    iwarnp = vcopt%iwarni
523
524    call vitens &
525    ( viscce , iwarnp ,             &
526      weighf , weighb ,             &
527      viscf  , viscb  )
528
529  endif
530
531else
532
533  do ifac = 1, nfac
534    viscf(ifac) = 0.d0
535  enddo
536  do ifac = 1, nfabor
537    viscb(ifac) = 0.d0
538  enddo
539
540endif
541
542if (ippmod(idarcy).eq.1) then
543  call field_get_key_struct_gwf_soilwater_partition(iflid, sorption_scal)
544  call field_get_val_s(sorption_scal%idel, cpro_delay)
545  call field_get_val_prev_s(sorption_scal%idel, cproa_delay)
546  call field_get_val_prev_s_by_name('saturation', cproa_sat)
547  call field_get_val_s_by_name('saturation', cpro_sat)
548endif
549
550! Not Darcy
551if (ippmod(idarcy).eq.-1) then
552  ! --> Unsteady term and mass aggregation term
553  do iel = 1, ncel
554    do isou = 1, 3
555      fimp(isou,isou,iel) = fimp(isou,isou,iel)                               &
556                          + vcopt%istat*pcrom(iel)*cell_f_vol(iel)/dt(iel)
557    enddo
558  enddo
559! Darcy : we take into account the porosity and delay for underground transport
560else
561  do iel = 1, ncel
562    do isou = 1, 3
563      smbrv(isou,iel) = smbrv(isou,iel)*cpro_delay(iel)*cpro_sat(iel)
564      fimp(isou,isou,iel) = (fimp(isou,isou,iel)                              &
565                          + vcopt%istat*pcrom(iel)*volume(iel)/dt(iel) )      &
566                            * cpro_delay(iel)*cpro_sat(iel)
567    enddo
568  enddo
569endif
570
571! Scalar with a Drift:
572! compute the convective flux
573!----------------------------
574
575call field_get_key_int(iflid, keydri, iscdri)
576
577if (iscdri.ge.1) then
578  allocate(divflu(ncelet))
579
580  call driflu &
581  ( iflid  ,                                                       &
582    dt     ,                                                       &
583    imasfl , bmasfl ,                                              &
584    divflu )
585
586  iconvp = vcopt%iconv
587  thetap = vcopt%thetav
588
589  ! NB: if the porosity module is swiched on, the the porosity is already
590  ! taken into account in divflu
591
592  ! --> mass aggregation term
593  do iel = 1, ncel
594    do isou = 1, 3
595      fimp(isou,isou,iel) = fimp(isou,isou,iel) + iconvp*thetap*divflu(iel)
596      smbrv(isou,iel) = smbrv(isou,iel) - iconvp*divflu(iel)*cvara_var(isou,iel)
597    enddo
598  enddo
599
600  deallocate(divflu)
601endif
602
603! Darcy:
604! This step is necessary because the divergence of the
605! hydraulic head (or pressure), which is taken into account as the
606! 'aggregation term' via
607! cs_equation_iterative_solve_scalar -> cs_balance_scalar,
608! is not exactly equal to the loss of mass of water in the unsteady case
609! (just as far as the precision of the Newton scheme is good),
610! and does not take into account the sorbption of the tracer
611! (represented by the 'delay' coefficient).
612! We choose to 'cancel' the aggregation term here, and to add to smbr
613! (right hand side of the transport equation)
614! the necessary correction to take sorption into account and get the exact
615! conservation of the mass of tracer.
616
617if (ippmod(idarcy).eq.1) then
618  if (darcy_unsteady.eq.1) then
619    call divmas(1, imasfl , bmasfl , diverg)
620    do iel = 1, ncel
621      do isou = 1, 3
622        smbrv(isou,iel) = smbrv(isou,iel) - diverg(iel)*cvar_var(isou,iel)     &
623          + cell_f_vol(iel)/dt(iel)*cvar_var(isou,iel)                         &
624          *( cproa_delay(iel)*cproa_sat(iel)                                   &
625           - cpro_delay(iel)*cpro_sat(iel) )
626      enddo
627    enddo
628    do iel = 1, ncel
629      do isou = 1, 3
630        fimp(isou,isou,iel) = fimp(isou,isou,iel) + thetv*diverg(iel)
631      enddo
632    enddo
633  endif
634endif
635
636!===============================================================================
637! 3. Solving
638!===============================================================================
639
640iescap = 0
641! all boundary convective flux with upwind
642icvflb = 0
643! transposed gradient term only for NS
644ivissv = 0
645
646call field_get_coefa_v(iflid, coefap)
647call field_get_coefb_v(iflid, coefbp)
648call field_get_coefaf_v(iflid, cofafp)
649call field_get_coefbf_v(iflid, cofbfp)
650
651vcopt_loc = vcopt
652
653vcopt_loc%istat  = -1
654vcopt_loc%idifft = -1
655vcopt_loc%iwgrec = 0
656vcopt_loc%thetav = thetv
657vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field
658
659p_k_value => vcopt_loc
660c_k_value = equation_param_from_vcopt(c_loc(p_k_value))
661
662call cs_equation_iterative_solve_vector                     &
663 ( idtvar , iterns ,                                        &
664   iflid  , c_null_char ,                                   &
665   ivissv , iescap , c_k_value       ,                      &
666   cvara_var       , cvara_var       ,                      &
667   coefap , coefbp , cofafp , cofbfp ,                      &
668   imasfl , bmasfl , viscf  ,                               &
669   viscb  , viscf  , viscb  , rvoid  ,                      &
670   rvoid  , viscce , weighf , weighb ,                      &
671   icvflb , ivoid  ,                                        &
672   fimp   , smbrv  , cvar_var        , rvoid )
673
674!===============================================================================
675! 4. Writing
676!===============================================================================
677
678! BILAN EXPLICITE (VOIR CODITS : ON ENLEVE L'INCREMENT)
679! Ceci devrait etre valable avec le theta schema sur les Termes source
680
681if (vcopt%iwarni.ge.2) then
682  if (vcopt%nswrsm.gt.1) then
683    ibcl = 1
684  else
685    ibcl = 0
686  endif
687  do iel = 1, ncel
688    do isou = 1, 3
689      smbrv(isou,iel) = smbrv(isou,iel)                                    &
690                      - vcopt%istat*(pcrom(iel)/dt(iel))*cell_f_vol(iel)   &
691                        *(cvar_var(isou,iel)-cvara_var(isou,iel))*ibcl
692    enddo
693  enddo
694  sclnor = sqrt(cs_gdot(3*ncel,smbrv,smbrv))
695  write(nfecra,1200)chaine(1:16) ,sclnor
696endif
697
698! Free memory
699deallocate(smbrv, fimp)
700deallocate(weighf, weighb)
701if (allocated(viscce)) deallocate(viscce)
702if (allocated(diverg)) deallocate(diverg)
703
704!--------
705! Formats
706!--------
707
708 1000 format(/,                                                   &
709'   ** SOLVING VARIABLE ',A16                                  ,/,&
710'      ----------------'                                       ,/)
711 1200 format(1X,A16,' : EXPLICIT BALANCE = ',E14.5)
712! 9000 format( &
713!'@'                                                            ,/,&
714!'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
715!'@'                                                            ,/,&
716!'@ @@ WARNING: ERROR IN COVOFV'                                ,/,&
717!'@    ========'                                                ,/,&
718!'@    IVARSC MUST BE A STRICTLY POSITIVE INTEGER'              ,/,&
719!'@    ITS VALUE IS ',I10                                       ,/,&
720!'@'                                                            ,/,&
721!'@  The calculation will not be run.'                          ,/,&
722!'@'                                                            ,/,&
723!'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,&
724!'@                                                            ',/)
725
726!----
727! End
728!----
729
730return
731
732end subroutine
733