← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/cfbl/cfdivs.f90

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-24 00:00:08 UTC
  • mfrom: (6.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20111124000008-2vo99e38267942q5
Tags: 2.1.0-3
Install a missing file

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/cfbl/cfdivs.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 cfdivs &
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 , iphas  ,                   &
36
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
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
   icepdc , icetsm , itypsm ,                                     &
40
28
   rtp    , propce , propfa , propfb ,                            &
41
29
   coefa  , coefb  , ckupdc , smacel ,                            &
42
 
   diverg , ux     , uy     , uz     ,                            &
43
 
   vistot ,                                                       &
44
 
   w1     , w2     , w3     , w4     , w5     , w6     ,          &
45
 
   rdevel , rtuser , ra     )
 
30
   diverg , ux     , uy     , uz     )
46
31
 
47
32
!===============================================================================
48
33
! FONCTION :
73
58
!__________________.____._____.________________________________________________.
74
59
! name             !type!mode ! role                                           !
75
60
!__________________!____!_____!________________________________________________!
76
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
77
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
78
 
! ndim             ! i  ! <-- ! spatial dimension                              !
79
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
80
 
! ncel             ! i  ! <-- ! number of cells                                !
81
 
! nfac             ! i  ! <-- ! number of interior faces                       !
82
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
83
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
84
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
85
 
! nnod             ! i  ! <-- ! number of vertices                             !
86
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
87
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
88
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
89
61
! nvar             ! i  ! <-- ! total number of variables                      !
90
62
! nscal            ! i  ! <-- ! total number of scalars                        !
91
 
! nphas            ! i  ! <-- ! number of phases                               !
92
63
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
93
64
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
94
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
95
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
96
 
! iphas            ! i  ! <-- ! phase number                                   !
97
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
98
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
99
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
100
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
101
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
102
 
!  (nfml, nprfml)  !    !     !                                                !
103
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
104
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
105
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
106
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
107
65
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
108
66
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
109
67
! itypsm           ! te ! <-- ! type de source de masse pour les               !
110
68
! (ncesmp,nvar)    !    !     !  variables (cf. ustsma)                        !
111
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
112
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
113
 
! ia(*)            ! ia ! --- ! main integer work array                        !
114
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
115
 
!  (ndim, ncelet)  !    !     !                                                !
116
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
117
 
!  (ndim, nfac)    !    !     !                                                !
118
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
119
 
!  (ndim, nfabor)  !    !     !                                                !
120
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
121
 
!  (ndim, nfac)    !    !     !                                                !
122
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
123
 
!  (ndim, nfabor)  !    !     !                                                !
124
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
125
 
!  (ndim, nnod)    !    !     !                                                !
126
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
127
69
! rtp              ! tr ! <-- ! variables de calcul au centre des              !
128
70
! (ncelet,*)       !    !     !    cellules (instant courant)                  !
129
71
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
138
80
!                  !    !     !  pour ivar=ipr, smacel=flux de masse           !
139
81
! diverg(ncelet    ! tr ! --> ! div(sigma.u)                                   !
140
82
! ux,y,z(ncelet    ! tr ! <-- ! composantes du vecteur u                       !
141
 
! vistot(ncelet    ! tr ! --- ! tableau de travail pour mu                     !
142
 
! w1...6(ncelet    ! tr ! --- ! tableau de travail                             !
143
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
144
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
145
 
! ra(*)            ! ra ! --- ! main real work array                           !
146
83
!__________________!____!_____!________________________________________________!
147
84
 
148
85
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
152
89
 
153
90
!===============================================================================
154
91
 
 
92
!===============================================================================
 
93
! Module files
 
94
!===============================================================================
 
95
 
 
96
use paramx
 
97
use cstphy
 
98
use entsor
 
99
use numvar
 
100
use optcal
 
101
use parall
 
102
use period
 
103
use ppppar
 
104
use ppthch
 
105
use ppincl
 
106
use mesh
 
107
 
 
108
!===============================================================================
 
109
 
155
110
implicit none
156
111
 
157
 
!===============================================================================
158
 
! Common blocks
159
 
!===============================================================================
160
 
 
161
 
include "paramx.h"
162
 
include "cstphy.h"
163
 
include "entsor.h"
164
 
include "numvar.h"
165
 
include "optcal.h"
166
 
include "vector.h"
167
 
include "period.h"
168
 
include "parall.h"
169
 
include "ppppar.h"
170
 
include "ppthch.h"
171
 
include "ppincl.h"
172
 
 
173
 
!===============================================================================
174
 
 
175
112
! Arguments
176
113
 
177
 
integer          idbia0 , idbra0
178
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
179
 
integer          nfml   , nprfml
180
 
integer          nnod   , lndfac , lndfbr , ncelbr
181
 
integer          nvar   , nscal  , nphas
 
114
integer          nvar   , nscal
182
115
integer          ncepdp , ncesmp
183
 
integer          nideve , nrdeve , nituse , nrtuse , iphas
184
116
 
185
 
integer          ifacel(2,nfac) , ifabor(nfabor)
186
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
187
 
integer          iprfml(nfml,nprfml)
188
 
integer          ipnfac(nfac+1), nodfac(lndfac)
189
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
190
117
integer          icepdc(ncepdp)
191
118
integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
192
 
integer          idevel(nideve), ituser(nituse)
193
 
integer          ia(*)
194
119
 
195
 
double precision xyzcen(ndim,ncelet)
196
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
197
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
198
 
double precision xyznod(ndim,nnod), volume(ncelet)
199
120
double precision rtp(ncelet,*)
200
121
double precision propce(ncelet,*)
201
122
double precision propfa(nfac,*), propfb(nfabor,*)
203
124
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
204
125
double precision diverg(ncelet)
205
126
double precision ux(ncelet), uy(ncelet), uz(ncelet)
206
 
double precision vistot(ncelet)
207
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
208
 
double precision w4(ncelet), w5(ncelet), w6(ncelet)
209
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
210
127
 
211
128
! Local variables
212
129
 
213
 
integer          idebia, idebra
214
130
integer          iccocg, inc, iel, ifac, ivar, isou, ii, jj
215
 
integer          iuiph, iviph, iwiph
216
131
integer          iclvar
217
132
integer          nswrgp, imligp, iwarnp
218
133
integer          ipcvis, ipcvst, ipcvsv
219
 
integer          idimte, itenso, iphydp
 
134
 
220
135
double precision epsrgp, climgp, extrap
221
136
double precision vecfac, visttt
222
137
 
 
138
double precision, allocatable, dimension(:) :: vistot
 
139
double precision, allocatable, dimension(:,:) :: grad
 
140
double precision, allocatable, dimension(:) :: w4, w5, w6
 
141
 
223
142
!===============================================================================
224
143
 
225
144
!===============================================================================
226
145
! 1.  INITIALISATION
227
146
!===============================================================================
228
147
 
229
 
idebia = idbia0
230
 
idebra = idbra0
231
 
 
232
 
iuiph  = iu(iphas)
233
 
iviph  = iv(iphas)
234
 
iwiph  = iw(iphas)
235
 
 
236
 
ipcvis = ipproc(iviscl(iphas))
237
 
ipcvst = ipproc(ivisct(iphas))
238
 
if(iviscv(iphas).gt.0) then
239
 
  ipcvsv = ipproc(iviscv(iphas))
 
148
! Allocate temporary arrays
 
149
allocate(vistot(ncelet))
 
150
allocate(grad(ncelet,3))
 
151
 
 
152
! Allocate work arrays
 
153
allocate(w4(ncelet), w5(ncelet), w6(ncelet))
 
154
 
 
155
 
 
156
ipcvis = ipproc(iviscl)
 
157
ipcvst = ipproc(ivisct)
 
158
if(iviscv.gt.0) then
 
159
  ipcvsv = ipproc(iviscv)
240
160
else
241
161
  ipcvsv = 0
242
162
endif
244
164
 
245
165
! --- Calcul de la viscosite totale
246
166
 
247
 
if (itytur(iphas).eq.3 ) then
 
167
if (itytur.eq.3 ) then
248
168
  do iel = 1, ncel
249
169
    vistot(iel) = propce(iel,ipcvis)
250
170
  enddo
260
180
!      cellules halo (calcul sur le halo, exceptionnellement).
261
181
!    Pour le parallelisme, on s'aligne sur la sequence ainsi definie.
262
182
 
263
 
! ---> TRAITEMENT DU PARALLELISME
264
 
 
265
 
if(irangp.ge.0) then
266
 
  call parcom (vistot)
267
 
  !==========
268
 
  if(ipcvsv.gt.0) then
269
 
    call parcom (propce(1,ipcvsv))
270
 
    !==========
271
 
  endif
272
 
endif
273
 
 
274
 
! ---> TRAITEMENT DE LA PERIODICITE
275
 
 
276
 
if(iperio.eq.1) then
277
 
  idimte = 0
278
 
  itenso = 0
279
 
  call percom                                                     &
280
 
  !==========
281
 
( idimte , itenso ,                                               &
282
 
  vistot , vistot , vistot ,                                      &
283
 
  vistot , vistot , vistot ,                                      &
284
 
  vistot , vistot , vistot )
285
 
  if(ipcvsv.gt.0) then
286
 
    call percom                                                   &
287
 
    !==========
288
 
( idimte , itenso ,                                               &
289
 
  propce(1,ipcvsv) , propce(1,ipcvsv) , propce(1,ipcvsv) ,        &
290
 
  propce(1,ipcvsv) , propce(1,ipcvsv) , propce(1,ipcvsv) ,        &
291
 
  propce(1,ipcvsv) , propce(1,ipcvsv) , propce(1,ipcvsv) )
 
183
! ---> TRAITEMENT DU PARALLELISME ET DE LA PERIODICITE
 
184
 
 
185
if (irangp.ge.0.or.iperio.eq.1) then
 
186
  call synsca(vistot)
 
187
  !==========
 
188
  if (ipcvsv.gt.0) then
 
189
    call synsca(propce(1,ipcvsv))
 
190
    !==========
292
191
  endif
293
192
endif
294
193
 
301
200
 
302
201
do isou = 1, 3
303
202
 
304
 
  if (isou.eq.1) ivar = iuiph
305
 
  if (isou.eq.2) ivar = iviph
306
 
  if (isou.eq.3) ivar = iwiph
 
203
  if (isou.eq.1) ivar = iu
 
204
  if (isou.eq.2) ivar = iv
 
205
  if (isou.eq.3) ivar = iw
307
206
 
308
207
! Ceci pointe eventuellement sur ICLRTP(IVAR,ICOEF)
309
208
  iclvar = iclrtp(ivar,icoeff)
318
217
  epsrgp = epsrgr(ivar)
319
218
  climgp = climgr(ivar)
320
219
  extrap = extrag(ivar)
321
 
  iphydp = 0
322
220
 
323
221
  call grdcel                                                     &
324
222
  !==========
325
 
 ( idebia , idebra ,                                              &
326
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
327
 
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
328
 
   nideve , nrdeve , nituse , nrtuse ,                            &
329
 
   ivar   , imrgra , inc    , iccocg , nswrgp , imligp , iphydp , &
 
223
 ( ivar   , imrgra , inc    , iccocg , nswrgp , imligp ,          &
330
224
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
331
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
332
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
333
 
   idevel , ituser , ia     ,                                     &
334
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
335
 
   w6     , w6     , w6     ,                                     &
336
225
   rtp(1,ivar)     , coefa(1,iclvar) , coefb(1,iclvar) ,          &
337
 
   w1     , w2     , w3     ,                                     &
338
 
!        ------   ------   ------
339
 
   w4     , w5     , w6     ,                                     &
340
 
   rdevel , rtuser , ra     )
 
226
   grad   )
341
227
 
342
228
 
343
229
 
349
235
    if    (isou.eq.1) then
350
236
      do iel = 1, ncelet
351
237
        visttt = propce(iel,ipcvsv) - 2.d0/3.d0*vistot(iel)
352
 
        w4(iel) = vistot(iel)*( 2.d0*w1(iel)*ux(iel)              &
353
 
                                   + w2(iel)*uy(iel)              &
354
 
                                   + w3(iel)*uz(iel) )            &
355
 
                     + visttt*w1(iel)*ux(iel)
356
 
        w5(iel) = vistot(iel)*w2(iel)*ux(iel)                     &
357
 
                     + visttt*w1(iel)*uy(iel)
358
 
        w6(iel) = vistot(iel)*w3(iel)*ux(iel)                     &
359
 
                     + visttt*w1(iel)*uz(iel)
 
238
        w4(iel) = vistot(iel)*( 2.d0*grad(iel,1)*ux(iel)          &
 
239
                                   + grad(iel,2)*uy(iel)          &
 
240
                                   + grad(iel,3)*uz(iel) )        &
 
241
                     + visttt*grad(iel,1)*ux(iel)
 
242
        w5(iel) = vistot(iel)*grad(iel,2)*ux(iel)                 &
 
243
                     + visttt*grad(iel,1)*uy(iel)
 
244
        w6(iel) = vistot(iel)*grad(iel,3)*ux(iel)                 &
 
245
                     + visttt*grad(iel,1)*uz(iel)
360
246
      enddo
361
247
 
362
248
    elseif(isou.eq.2) then
363
249
      do iel = 1, ncelet
364
250
        visttt = propce(iel,ipcvsv) - 2.d0/3.d0*vistot(iel)
365
 
        w4(iel) = vistot(iel)*w1(iel)*uy(iel)                     &
366
 
                     + visttt*w2(iel)*ux(iel)
367
 
        w5(iel) = vistot(iel)*(      w1(iel)*ux(iel)              &
368
 
                              + 2.d0*w2(iel)*uy(iel)              &
369
 
                                   + w3(iel)*uz(iel) )            &
370
 
                     + visttt*w2(iel)*uy(iel)
371
 
        w6(iel) = vistot(iel)*w3(iel)*uy(iel)                     &
372
 
                     + visttt*w2(iel)*uz(iel)
 
251
        w4(iel) = vistot(iel)*grad(iel,1)*uy(iel)                 &
 
252
                     + visttt*grad(iel,2)*ux(iel)
 
253
        w5(iel) = vistot(iel)*(      grad(iel,1)*ux(iel)          &
 
254
                              + 2.d0*grad(iel,2)*uy(iel)          &
 
255
                                   + grad(iel,3)*uz(iel) )        &
 
256
                     + visttt*grad(iel,2)*uy(iel)
 
257
        w6(iel) = vistot(iel)*grad(iel,3)*uy(iel)                 &
 
258
                     + visttt*grad(iel,2)*uz(iel)
373
259
      enddo
374
260
 
375
261
    elseif(isou.eq.3) then
376
262
      do iel = 1, ncelet
377
263
        visttt = propce(iel,ipcvsv) - 2.d0/3.d0*vistot(iel)
378
 
        w4(iel) = vistot(iel)*w1(iel)*uz(iel)                     &
379
 
                     + visttt*w3(iel)*ux(iel)
380
 
        w5(iel) = vistot(iel)*w2(iel)*uz(iel)                     &
381
 
                     + visttt*w3(iel)*uy(iel)
382
 
        w6(iel) = vistot(iel)*(      w1(iel)*ux(iel)              &
383
 
                                   + w2(iel)*uy(iel)              &
384
 
                              + 2.d0*w3(iel)*uz(iel) )            &
385
 
                     + visttt*w3(iel)*uz(iel)
 
264
        w4(iel) = vistot(iel)*grad(iel,1)*uz(iel)                 &
 
265
                     + visttt*grad(iel,3)*ux(iel)
 
266
        w5(iel) = vistot(iel)*grad(iel,2)*uz(iel)                 &
 
267
                     + visttt*grad(iel,3)*uy(iel)
 
268
        w6(iel) = vistot(iel)*(      grad(iel,1)*ux(iel)          &
 
269
                                   + grad(iel,2)*uy(iel)          &
 
270
                              + 2.d0*grad(iel,3)*uz(iel) )        &
 
271
                     + visttt*grad(iel,3)*uz(iel)
386
272
      enddo
387
273
 
388
274
    endif
391
277
 
392
278
    if    (isou.eq.1) then
393
279
      do iel = 1, ncelet
394
 
        visttt = viscv0(iphas) - 2.d0/3.d0*vistot(iel)
395
 
        w4(iel) = vistot(iel)*( 2.d0*w1(iel)*ux(iel)              &
396
 
                                   + w2(iel)*uy(iel)              &
397
 
                                   + w3(iel)*uz(iel) )            &
398
 
                     + visttt*w1(iel)*ux(iel)
399
 
        w5(iel) = vistot(iel)*w2(iel)*ux(iel)                     &
400
 
                     + visttt*w1(iel)*uy(iel)
401
 
        w6(iel) = vistot(iel)*w3(iel)*ux(iel)                     &
402
 
                     + visttt*w1(iel)*uz(iel)
 
280
        visttt = viscv0 - 2.d0/3.d0*vistot(iel)
 
281
        w4(iel) = vistot(iel)*( 2.d0*grad(iel,1)*ux(iel)          &
 
282
                                   + grad(iel,2)*uy(iel)          &
 
283
                                   + grad(iel,3)*uz(iel) )        &
 
284
                     + visttt*grad(iel,1)*ux(iel)
 
285
        w5(iel) = vistot(iel)*grad(iel,2)*ux(iel)                 &
 
286
                     + visttt*grad(iel,1)*uy(iel)
 
287
        w6(iel) = vistot(iel)*grad(iel,3)*ux(iel)                 &
 
288
                     + visttt*grad(iel,1)*uz(iel)
403
289
      enddo
404
290
 
405
291
    elseif(isou.eq.2) then
406
292
      do iel = 1, ncelet
407
 
        visttt = viscv0(iphas) - 2.d0/3.d0*vistot(iel)
408
 
        w4(iel) = vistot(iel)*w1(iel)*uy(iel)                     &
409
 
                     + visttt*w2(iel)*ux(iel)
410
 
        w5(iel) = vistot(iel)*(      w1(iel)*ux(iel)              &
411
 
                              + 2.d0*w2(iel)*uy(iel)              &
412
 
                                   + w3(iel)*uz(iel) )            &
413
 
                     + visttt*w2(iel)*uy(iel)
414
 
        w6(iel) = vistot(iel)*w3(iel)*uy(iel)                     &
415
 
                     + visttt*w2(iel)*uz(iel)
 
293
        visttt = viscv0 - 2.d0/3.d0*vistot(iel)
 
294
        w4(iel) = vistot(iel)*grad(iel,1)*uy(iel)                 &
 
295
                     + visttt*grad(iel,2)*ux(iel)
 
296
        w5(iel) = vistot(iel)*(      grad(iel,1)*ux(iel)          &
 
297
                              + 2.d0*grad(iel,2)*uy(iel)          &
 
298
                                   + grad(iel,3)*uz(iel) )        &
 
299
                     + visttt*grad(iel,2)*uy(iel)
 
300
        w6(iel) = vistot(iel)*grad(iel,3)*uy(iel)                 &
 
301
                     + visttt*grad(iel,2)*uz(iel)
416
302
      enddo
417
303
 
418
304
    elseif(isou.eq.3) then
419
305
      do iel = 1, ncelet
420
 
        visttt = viscv0(iphas) - 2.d0/3.d0*vistot(iel)
421
 
        w4(iel) = vistot(iel)*w1(iel)*uz(iel)                     &
422
 
                     + visttt*w3(iel)*ux(iel)
423
 
        w5(iel) = vistot(iel)*w2(iel)*uz(iel)                     &
424
 
                     + visttt*w3(iel)*uy(iel)
425
 
        w6(iel) = vistot(iel)*(      w1(iel)*ux(iel)              &
426
 
                                   + w2(iel)*uy(iel)              &
427
 
                              + 2.d0*w3(iel)*uz(iel) )            &
428
 
                     + visttt*w3(iel)*uz(iel)
 
306
        visttt = viscv0 - 2.d0/3.d0*vistot(iel)
 
307
        w4(iel) = vistot(iel)*grad(iel,1)*uz(iel)                 &
 
308
                     + visttt*grad(iel,3)*ux(iel)
 
309
        w5(iel) = vistot(iel)*grad(iel,2)*uz(iel)                 &
 
310
                     + visttt*grad(iel,3)*uy(iel)
 
311
        w6(iel) = vistot(iel)*(      grad(iel,1)*ux(iel)          &
 
312
                                   + grad(iel,2)*uy(iel)          &
 
313
                              + 2.d0*grad(iel,3)*uz(iel) )        &
 
314
                     + visttt*grad(iel,3)*uz(iel)
429
315
      enddo
430
316
 
431
317
    endif
445
331
 
446
332
 
447
333
 
448
 
  if(ivecti.eq.1) then
449
 
 
450
 
!CDIR NODEP
451
 
    do ifac = 1, nfac
452
 
      ii = ifacel(1,ifac)
453
 
      jj = ifacel(2,ifac)
454
 
!MO             VECFAC = SURFAC(ISOU,IFAC)
455
 
!MO     &                  *(POND(IFAC)*W4(II)+(1.D0-POND(IFAC))*W4(JJ))
456
 
      vecfac = surfac(1,ifac)*(w4(ii)+w4(jj))*0.5d0               &
457
 
             + surfac(2,ifac)*(w5(ii)+w5(jj))*0.5d0               &
458
 
             + surfac(3,ifac)*(w6(ii)+w6(jj))*0.5d0
459
 
      diverg(ii) = diverg(ii) + vecfac
460
 
      diverg(jj) = diverg(jj) - vecfac
461
 
    enddo
462
 
 
463
 
  else
464
 
 
465
 
! VECTORISATION NON FORCEE
466
 
    do ifac = 1, nfac
467
 
      ii = ifacel(1,ifac)
468
 
      jj = ifacel(2,ifac)
469
 
!MO             VECFAC = SURFAC(ISOU,IFAC)
470
 
!MO     &                  *(POND(IFAC)*W4(II)+(1.D0-POND(IFAC))*W4(JJ))
471
 
      vecfac = surfac(1,ifac)*(w4(ii)+w4(jj))*0.5d0               &
472
 
             + surfac(2,ifac)*(w5(ii)+w5(jj))*0.5d0               &
473
 
             + surfac(3,ifac)*(w6(ii)+w6(jj))*0.5d0
474
 
      diverg(ii) = diverg(ii) + vecfac
475
 
      diverg(jj) = diverg(jj) - vecfac
476
 
    enddo
477
 
 
478
 
  endif
 
334
  do ifac = 1, nfac
 
335
    ii = ifacel(1,ifac)
 
336
    jj = ifacel(2,ifac)
 
337
!MO             VECFAC = SURFAC(ISOU,IFAC)
 
338
!MO     &                  *(POND(IFAC)*W4(II)+(1.D0-POND(IFAC))*W4(JJ))
 
339
    vecfac = surfac(1,ifac)*(w4(ii)+w4(jj))*0.5d0               &
 
340
         + surfac(2,ifac)*(w5(ii)+w5(jj))*0.5d0               &
 
341
         + surfac(3,ifac)*(w6(ii)+w6(jj))*0.5d0
 
342
    diverg(ii) = diverg(ii) + vecfac
 
343
    diverg(jj) = diverg(jj) - vecfac
 
344
  enddo
479
345
 
480
346
 
481
347
! --- Assemblage sur les faces de bord
482
348
 
483
 
  if(ivectb.eq.1) then
484
 
 
485
 
!CDIR NODEP
486
 
    do ifac = 1, nfabor
487
 
      ii = ifabor(ifac)
488
 
      vecfac = surfbo(1,ifac)*w4(ii)                              &
489
 
             + surfbo(2,ifac)*w5(ii)                              &
490
 
             + surfbo(3,ifac)*w6(ii)
491
 
      diverg(ii) = diverg(ii) + vecfac
492
 
    enddo
493
 
 
494
 
  else
495
 
 
496
 
! VECTORISATION NON FORCEE
497
 
    do ifac = 1, nfabor
498
 
      ii = ifabor(ifac)
499
 
      vecfac = surfbo(1,ifac)*w4(ii)                              &
500
 
             + surfbo(2,ifac)*w5(ii)                              &
501
 
             + surfbo(3,ifac)*w6(ii)
502
 
      diverg(ii) = diverg(ii) + vecfac
503
 
    enddo
504
 
 
505
 
  endif
 
349
  do ifac = 1, nfabor
 
350
    ii = ifabor(ifac)
 
351
    vecfac = surfbo(1,ifac)*w4(ii)                              &
 
352
         + surfbo(2,ifac)*w5(ii)                              &
 
353
         + surfbo(3,ifac)*w6(ii)
 
354
    diverg(ii) = diverg(ii) + vecfac
 
355
  enddo
506
356
 
507
357
enddo
508
358
 
 
359
! Free memory
 
360
deallocate(w4, w5, w6)
509
361
 
510
362
return
511
363