← Back to branch summary

~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to src/cplv/cptsvc.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-01 17:43:32 UTC
  • mto: (6.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 11.
  • Revision ID: package-import@ubuntu.com-20111101174332-tl4vk45no0x3emc3
Tags: upstream-2.1.0
Import upstream version 2.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/cplv/cptsvc.f90
1
1
!-------------------------------------------------------------------------------
2
2
 
3
 
!     This file is part of the Code_Saturne Kernel, element of the
4
 
!     Code_Saturne CFD tool.
5
 
 
6
 
!     Copyright (C) 1998-2009 EDF S.A., France
7
 
 
8
 
!     contact: saturne-support@edf.fr
9
 
 
10
 
!     The Code_Saturne Kernel is free software; you can redistribute it
11
 
!     and/or modify it under the terms of the GNU General Public License
12
 
!     as published by the Free Software Foundation; either version 2 of
13
 
!     the License, or (at your option) any later version.
14
 
 
15
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
16
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18
 
!     GNU General Public License for more details.
19
 
 
20
 
!     You should have received a copy of the GNU General Public License
21
 
!     along with the Code_Saturne Kernel; if not, write to the
22
 
!     Free Software Foundation, Inc.,
23
 
!     51 Franklin St, Fifth Floor,
24
 
!     Boston, MA  02110-1301  USA
 
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
4
!
 
5
! Copyright (C) 1998-2011 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.
25
20
 
26
21
!-------------------------------------------------------------------------------
27
22
 
28
23
subroutine cptsvc &
29
24
!================
30
25
 
31
 
 ( idbia0 , idbra0 ,                                              &
32
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
33
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
34
 
   nvar   , nscal  , nphas  , ncepdp , ncesmp ,                   &
35
 
   nideve , nrdeve , nituse , nrtuse , iscal  , iscala ,          &
36
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , itypfb ,          &
37
 
   ipnfac , nodfac , ipnfbr , nodfbr , icepdc , icetsm , itypsm , &
38
 
   idevel , ituser , ia     ,                                     &
39
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
 
26
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
 
27
   iscal  , iscala ,                                              &
 
28
   itypfb ,                                                       &
 
29
   icepdc , icetsm , itypsm ,                                     &
40
30
   dt     , rtpa   , rtp    , propce , propfa , propfb ,          &
41
31
   coefa  , coefb  ,                                              &
42
 
   smbrs  , rovsdt ,                                              &
43
 
   wfb    ,                                                       &
44
 
   w1     , w2     , w3     , w4     , w5     ,                   &
45
 
   w6     , w7     , w8     ,                                     &
46
 
   rdevel , rtuser , ra     )
 
32
   smbrs  , rovsdt )
47
33
 
48
34
!===============================================================================
49
35
! FONCTION :
58
44
!__________________.____._____.________________________________________________.
59
45
! name             !type!mode ! role                                           !
60
46
!__________________!____!_____!________________________________________________!
61
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
62
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
63
 
! ndim             ! i  ! <-- ! spatial dimension                              !
64
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
65
 
! ncel             ! i  ! <-- ! number of cells                                !
66
 
! nfac             ! i  ! <-- ! number of interior faces                       !
67
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
68
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
69
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
70
 
! nnod             ! i  ! <-- ! number of vertices                             !
71
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
72
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
73
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
74
47
! nvar             ! i  ! <-- ! total number of variables                      !
75
48
! nscal            ! i  ! <-- ! total number of scalars                        !
76
 
! nphas            ! i  ! <-- ! number of phases                               !
77
49
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
78
50
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
79
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
80
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
81
51
! iscal            ! i  ! <-- ! scalar number                                  !
82
52
! iscala           ! e  ! <-- ! numero du scalaire associe                     !
83
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
84
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
85
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
86
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
87
 
! iprfml           ! te ! <-- ! proprietes d'une famille                       !
88
53
! itypfb           ! ia ! <-- ! boundary face types                            !
89
 
!  (nfabor, nphas) !    !     !                                                !
90
 
! nfml  ,nprfml    !    !     !                                                !
91
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
92
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
93
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
94
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
95
54
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
96
55
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
97
56
! itypsm           ! te ! <-- ! type de source de masse pour les               !
98
57
! (ncesmp,nvar)    !    !     !  variables (cf. ustsma)                        !
99
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
100
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
101
 
! ia(*)            ! ia ! --- ! main integer work array                        !
102
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
103
 
!  (ndim, ncelet)  !    !     !                                                !
104
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
105
 
!  (ndim, nfac)    !    !     !                                                !
106
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
107
 
!  (ndim, nfabor)  !    !     !                                                !
108
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
109
 
!  (ndim, nfac)    !    !     !                                                !
110
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
111
 
!  (ndim, nfabor)  !    !     !                                                !
112
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
113
 
!  (ndim, nnod)    !    !     !                                                !
114
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
115
58
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
116
59
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
117
60
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
122
65
!  (nfabor, *)     !    !     !                                                !
123
66
! smbrs(ncelet)    ! tr ! --> ! second membre explicite                        !
124
67
! rovsdt(ncelet    ! tr ! --> ! partie diagonale implicite                     !
125
 
! wfb(nfabor)      ! tr ! --- ! tableau de travail    faces de bord            !
126
 
! w1..8(ncelet)    ! tr ! --- ! tableau de travail    cellules                 !
127
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
128
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
129
 
! ra(*)            ! ra ! --- ! main real work array                           !
130
68
!__________________!____!_____!________________________________________________!
131
69
 
132
70
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
133
71
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
134
72
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
135
73
!            --- tableau de travail
136
 
!-------------------------------------------------------------------------------
 
74
!===============================================================================
 
75
 
 
76
!===============================================================================
 
77
! Module files
 
78
!===============================================================================
 
79
 
 
80
use paramx
 
81
use numvar
 
82
use entsor
 
83
use optcal
 
84
use cstphy
 
85
use cstnum
 
86
use parall
 
87
use period
 
88
use ppppar
 
89
use ppthch
 
90
use coincl
 
91
use cpincl
 
92
use ppincl
 
93
use mesh
 
94
 
137
95
!===============================================================================
138
96
 
139
97
implicit none
140
98
 
141
 
!===============================================================================
142
 
! Common blocks
143
 
!===============================================================================
144
 
 
145
 
include "paramx.h"
146
 
include "numvar.h"
147
 
include "entsor.h"
148
 
include "optcal.h"
149
 
include "cstphy.h"
150
 
include "cstnum.h"
151
 
include "parall.h"
152
 
include "period.h"
153
 
include "ppppar.h"
154
 
include "ppthch.h"
155
 
include "coincl.h"
156
 
include "cpincl.h"
157
 
include "ppincl.h"
158
 
 
159
 
!===============================================================================
160
 
 
161
99
! Arguments
162
100
 
163
 
integer          idbia0 , idbra0
164
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
165
 
integer          nfml   , nprfml
166
 
integer          nnod   , lndfac , lndfbr , ncelbr
167
 
integer          nvar   , nscal  , nphas
 
101
integer          nvar   , nscal
168
102
integer          ncepdp , ncesmp
169
 
integer          nideve , nrdeve , nituse , nrtuse
170
103
integer          iscal  , iscala
171
104
 
172
 
integer          ifacel(2,nfac) , ifabor(nfabor)
173
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
174
 
integer          iprfml(nfml,nprfml) , itypfb(nfabor,nphas)
175
 
integer          ipnfac(nfac+1), nodfac(lndfac)
176
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
 
105
integer          itypfb(nfabor)
177
106
integer          icepdc(ncepdp)
178
107
integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
179
 
integer          idevel(nideve)
180
 
integer          ituser(nituse), ia(*)
181
108
 
182
 
double precision xyzcen(ndim,ncelet)
183
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
184
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
185
 
double precision xyznod(ndim,nnod), volume(ncelet)
186
109
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
187
110
double precision propce(ncelet,*)
188
111
double precision propfa(nfac,*), propfb(nfabor,*)
189
112
double precision coefa(nfabor,*), coefb(nfabor,*)
190
113
double precision smbrs(ncelet), rovsdt(ncelet)
191
 
double precision wfb(nfabor)
192
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
193
 
double precision w4(ncelet), w5(ncelet), w6(ncelet)
194
 
double precision w7(ncelet), w8(ncelet)
195
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
196
114
 
197
115
! Local variables
198
116
 
199
 
integer          idebia , idebra
200
117
integer          ivar   , ivarsc , ivarut , ivar0
201
 
integer          iel    , iphas  , ifac
 
118
integer          iel    , ifac
202
119
integer          ipcrom , ipcvst
203
 
integer          ikiph  , ieiph  , iomgip , iphydp
204
 
integer          ir11ip , ir22ip , ir33ip
205
120
integer          ixchcl , ixckcl , ixnpcl , icla   , icha
206
121
integer          inc    , iccocg , nswrgp , imligp , iwarnp
207
 
integer          ifinra , icoefa , icoefb
208
 
integer          idimte , itenso
209
122
 
210
123
double precision x2     , xk     , xe     , rhovst
211
124
double precision epsrgp , climgp , extrap
212
125
 
 
126
double precision, allocatable, dimension(:) :: coefap, coefbp
 
127
double precision, allocatable, dimension(:,:) :: grad
 
128
double precision, allocatable, dimension(:) :: w1, w2
 
129
double precision, allocatable, dimension(:) :: w7, w8
213
130
 
214
131
!===============================================================================
215
132
 
217
134
! 1. INITIALISATION
218
135
!===============================================================================
219
136
 
220
 
idebia = idbia0
221
 
idebra = idbra0
 
137
! Initialize variables to avoid compiler warnings
 
138
 
 
139
xe = 0.d0
 
140
xk = 0.d0
 
141
 
 
142
! Memoire
 
143
 
222
144
 
223
145
! --- Numero du scalaire a traiter : ISCAL
224
146
 
233
155
  ivarsc = 0
234
156
endif
235
157
 
236
 
! --- Numero de phase associee au scalaire ISCAL
237
 
iphas = iphsca(iscal)
238
 
 
239
 
! --- Numero des variables de calcul
240
 
if ( itytur(iphas).eq.2 .or. iturb(iphas).eq.50 ) then
241
 
  ikiph  = ik  (iphas)
242
 
  ieiph  = iep (iphas)
243
 
elseif ( itytur(iphas).eq.3 ) then
244
 
  ir11ip = ir11(iphas)
245
 
  ir22ip = ir22(iphas)
246
 
  ir33ip = ir33(iphas)
247
 
  ieiph  = iep (iphas)
248
 
elseif ( iturb(iphas).eq.60 ) then
249
 
  ikiph  = ik  (iphas)
250
 
  iomgip = iomg(iphas)
251
 
endif
252
 
 
253
158
! --- Numero des grandeurs physiques
254
 
ipcrom = ipproc(irom(iphas))
255
 
ipcvst = ipproc(ivisct(iphas))
 
159
ipcrom = ipproc(irom)
 
160
ipcvst = ipproc(ivisct)
256
161
 
257
162
 
258
163
!===============================================================================
260
165
!    ET DE DISSIPATION
261
166
!===============================================================================
262
167
 
263
 
if ( itytur(iphas).eq.2 .or. itytur(iphas).eq.3                   &
264
 
     .or. iturb(iphas).eq.50 .or. iturb(iphas).eq.60) then
 
168
if ( itytur.eq.2 .or. itytur.eq.3                   &
 
169
     .or. iturb.eq.50 .or. iturb.eq.60) then
265
170
 
266
171
  inc = 1
267
172
  iccocg = 1
278
183
  climgp = climgr(ivarut)
279
184
  extrap = extrag(ivarut)
280
185
 
 
186
  ! Allocate work arrays
 
187
  allocate(w7(ncelet), w8(ncelet))
 
188
 
281
189
! --> Calcul de FIM et X2 dans W7 et W8
282
190
 
283
191
  do iel = 1, ncel
284
 
    w1(iel) = zero
285
 
    w2(iel) = zero
286
 
    w3(iel) = zero
287
 
    w4(iel) = zero
288
 
    w5(iel) = zero
289
 
    w6(iel) = zero
290
192
    w7(iel) = zero
291
193
    w8(iel) = 1.d0
292
194
  enddo
307
209
! ---- W7 = FJM (kg/kg du melange gazeux)
308
210
 
309
211
  if (ivarsc.eq.0) then
 
212
    ! Allocate work arrays
 
213
    allocate(w1(ncelet), w2(ncelet))
 
214
    do iel = 1, ncel
 
215
      w1(iel) = zero
 
216
      w2(iel) = zero
 
217
    enddo
310
218
    do icha = 1, ncharb
311
219
      do iel = 1, ncel
312
220
        w1(iel) =  w1(iel) + rtp(iel,isca(if1m(icha)))
318
226
              ( (w1(iel) + w2(iel) + rtp(iel,isca(if3m)))         &
319
227
                / w8(iel) )
320
228
    enddo
 
229
    ! Free some work arrays
 
230
    deallocate(w1, w2)
321
231
  else
322
232
    do iel = 1, ncel
323
233
      w7(iel) = rtp(iel,ivarsc) / w8(iel)
324
234
    enddo
325
235
  endif
326
236
 
327
 
! --> Calcul des COEFA et COEFB de FIM afin d'en calculer son gradient
328
 
!     On alloue localement 2 tableaux de NFABOR pour le calcul
329
 
!       de COEFA et COEFB de FIM
330
 
 
331
 
  icoefa = idebra
332
 
  icoefb = icoefa + nfabor
333
 
  ifinra = icoefb + nfabor
334
 
  CALL RASIZE ('CPTSVC',IFINRA)
335
 
  !==========
 
237
  ! Allocate temporary arrays
 
238
  allocate(coefap(nfabor), coefbp(nfabor))
336
239
 
337
240
  do ifac = 1, nfabor
338
 
    ra(icoefa+ifac-1) = zero
339
 
    ra(icoefb+ifac-1) = 1.d0
340
 
    if ( itypfb(ifac,iphas).eq.ientre ) then
341
 
      ra(icoefa+ifac-1) = zero
342
 
      ra(icoefb+ifac-1) = zero
343
 
      if (ivarsc.eq.0) ra(icoefa+ifac-1) = 1.d0
 
241
    coefap(ifac) = zero
 
242
    coefbp(ifac) = 1.d0
 
243
    if ( itypfb(ifac).eq.ientre ) then
 
244
      coefap(ifac) = zero
 
245
      coefbp(ifac) = zero
 
246
      if (ivarsc.eq.0) coefap(ifac) = 1.d0
344
247
    endif
345
248
  enddo
346
249
 
347
 
! En periodique et parallele, echange avant calcul du gradient
348
 
 
349
 
!    Parallele
350
 
  if(irangp.ge.0) then
351
 
    call parcom(w7)
352
 
    !==========
353
 
  endif
354
 
 
355
 
!    Periodique
356
 
  if(iperio.eq.1) then
357
 
    idimte = 0
358
 
    itenso = 0
359
 
    call percom                                                   &
360
 
    !==========
361
 
  ( idimte , itenso ,                                             &
362
 
    w7     , w7     , w7     ,                                    &
363
 
    w7     , w7     , w7     ,                                    &
364
 
    w7     , w7     , w7     )
365
 
  endif
 
250
  ! En periodique et parallele, echange avant calcul du gradient
 
251
  if (irangp.ge.0.or.iperio.eq.1) then
 
252
    call synsca(w7)
 
253
    !==========
 
254
  endif
 
255
 
 
256
  ! Allocate a temporary array for gradient computation
 
257
  allocate(grad(ncelet,3))
366
258
 
367
259
!  IVAR0 = 0 (indique pour la periodicite de rotation que la variable
368
260
!     n'est pas la vitesse ni Rij)
369
261
  ivar0  = 0
370
 
  iphydp = 0
371
262
  call grdcel                                                     &
372
263
  !==========
373
 
 ( idebia , ifinra ,                                              &
374
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
375
 
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
376
 
   nideve , nrdeve , nituse , nrtuse ,                            &
377
 
   ivar0  , imrgra , inc    , iccocg , nswrgp , imligp , iphydp , &
 
264
 ( ivar0  , imrgra , inc    , iccocg , nswrgp , imligp ,          &
378
265
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
379
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
380
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
381
 
   idevel , ituser , ia     ,                                     &
382
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
383
 
   w7     , w7     , w7     ,                                     &
384
 
   w7     , ra(icoefa) , ra(icoefb)  ,                            &
 
266
   w7     , coefap , coefbp ,                                     &
385
267
!          FIM      COEFA        COEFB
386
 
   w1     , w2     , w3     ,                                     &
387
 
!        ------   ------   ------
388
 
   w4     , w5     , w6     ,                                     &
389
 
   rdevel , rtuser , ra     )
 
268
   grad   )
 
269
 
 
270
  ! Free memory
 
271
  deallocate(coefap, coefbp)
390
272
 
391
273
  do iel = 1, ncel
392
 
    if ( itytur(iphas).eq.2 .or. iturb(iphas).eq.50 ) then
393
 
      xk = rtpa(iel,ikiph)
394
 
      xe = rtpa(iel,ieiph)
395
 
    elseif ( itytur(iphas).eq.3 ) then
396
 
      xk =                                                        &
397
 
       0.5d0*(rtpa(iel,ir11ip)+rtpa(iel,ir22ip)+rtpa(iel,ir33ip))
398
 
      xe = rtpa(iel,ieiph)
399
 
    elseif ( iturb(iphas).eq.60 ) then
400
 
      xk = rtpa(iel,ikiph)
401
 
      xe = cmu*xk*rtpa(iel,iomgip)
 
274
    if ( itytur.eq.2 .or. iturb.eq.50 ) then
 
275
      xk = rtpa(iel,ik)
 
276
      xe = rtpa(iel,iep)
 
277
    elseif ( itytur.eq.3 ) then
 
278
      xk = 0.5d0*(rtpa(iel,ir11)+rtpa(iel,ir22)+rtpa(iel,ir33))
 
279
      xe = rtpa(iel,iep)
 
280
    elseif ( iturb.eq.60 ) then
 
281
      xk = rtpa(iel,ik)
 
282
      xe = cmu*xk*rtpa(iel,iomg)
402
283
    endif
403
284
 
404
285
    rhovst = propce(iel,ipcrom)*xe/                               &
405
286
             (xk * rvarfl(iscal))*volume(iel)
406
287
    rovsdt(iel) = rovsdt(iel) + max(zero,rhovst)
407
 
    smbrs(iel) = smbrs(iel) +                                     &
408
 
                2.d0*propce(iel,ipcvst)*volume(iel)/sigmas(iscal) &
409
 
                * (w1(iel)**2 + w2(iel)**2 + w3(iel)**2) * w8(iel)&
 
288
    smbrs(iel) = smbrs(iel) +                                                 &
 
289
                2.d0*propce(iel,ipcvst)*volume(iel)/sigmas(iscal)             &
 
290
               * (grad(iel,1)**2 + grad(iel,2)**2 + grad(iel,3)**2) * w8(iel) &
410
291
                - rhovst*rtpa(iel,ivar)
411
292
  enddo
412
293
 
413
 
!     On libere COEFA COEFB
414
 
  ifinra = idebra
 
294
  ! Free memory
 
295
  deallocate(grad)
 
296
  deallocate(w7, w8)
415
297
 
416
298
endif
417
299
 
418
 
 
419
 
 
420
300
!--------
421
301
! FORMATS
422
302
!--------