← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/base/rijech.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/base/rijech.f90
1
 
!-------------------------------------------------------------------------------
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
25
 
 
26
 
!-------------------------------------------------------------------------------
27
 
 
28
 
subroutine rijech &
29
 
!================
30
 
 
31
 
 ( idbia0 , idbra0 ,                                              &
32
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
33
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
34
 
   nvar   , nscal  , nphas  ,                                     &
35
 
   nideve , nrdeve , nituse , nrtuse ,                            &
36
 
   iphas  , ivar   , isou   , ipp    ,                            &
37
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
38
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
39
 
   idevel , ituser , ia     ,                                     &
40
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
41
 
   rtp    , rtpa   , propce , propfa , propfb ,                   &
42
 
   coefa  , coefb  , produc , smbr   ,                            &
43
 
   coefax , coefbx ,                                              &
44
 
   produk , w2     , w3     , w4     , epsk   , w6     ,          &
45
 
   rdevel , rtuser , ra     )
46
 
 
47
 
!===============================================================================
48
 
! FONCTION :
49
 
! ----------
50
 
 
51
 
! TERMES D'ECHO DE PAROI
52
 
!   POUR Rij
53
 
! VAR  = R11 R22 R33 R12 R13 R23
54
 
! ISOU =  1   2   3   4   5   6
55
 
 
56
 
!-------------------------------------------------------------------------------
57
 
!ARGU                             ARGUMENTS
58
 
!__________________.____._____.________________________________________________.
59
 
! name             !type!mode ! role                                           !
60
 
!__________________!____!_____!________________________________________________!
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
 
! nvar             ! i  ! <-- ! total number of variables                      !
75
 
! nscal            ! i  ! <-- ! total number of scalars                        !
76
 
! nphas            ! i  ! <-- ! number of phases                               !
77
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
78
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
79
 
! iphas            ! i  ! <-- ! phase number                                   !
80
 
! ivar             ! i  ! <-- ! variable number                                !
81
 
! isou             ! e  ! <-- ! numero de passage                              !
82
 
! ipp              ! e  ! <-- ! numero de variable pour sorties post           !
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           ! ia ! <-- ! property numbers per family                    !
88
 
!  (nfml, nprfml)  !    !     !                                                !
89
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
90
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
91
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
92
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
93
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
94
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
95
 
! ia(*)            ! ia ! --- ! main integer work array                        !
96
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
97
 
!  (ndim, ncelet)  !    !     !                                                !
98
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
99
 
!  (ndim, nfac)    !    !     !                                                !
100
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
101
 
!  (ndim, nfabor)  !    !     !                                                !
102
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
103
 
!  (ndim, nfac)    !    !     !                                                !
104
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
105
 
!  (ndim, nfabor)  !    !     !                                                !
106
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
107
 
!  (ndim, nnod)    !    !     !                                                !
108
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
109
 
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
110
 
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
111
 
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
112
 
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
113
 
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
114
 
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
115
 
!  (nfabor, *)     !    !     !                                                !
116
 
! produc           ! tr ! <-- ! production                                     !
117
 
!  (6,ncelet)      !    !     !                                                !
118
 
! smbr(ncelet      ! tr ! <-- ! tableau de travail pour sec mem                !
119
 
! coefax,coefbx    ! tr ! --- ! tableau de travail pour les cond.              !
120
 
!  (nfabor)        !    !     !    aux limites de la dist. paroi               !
121
 
! produk(ncelet    ! tr ! --- ! tableau de travail production                  !
122
 
! epsk  (ncelet    ! tr ! --- ! tableau de travail epsilon/k                   !
123
 
! w2...6(ncelet    ! tr ! --- ! tableau de travail production                  !
124
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
125
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
126
 
! ra(*)            ! ra ! --- ! main real work array                           !
127
 
!__________________!____!_____!________________________________________________!
128
 
 
129
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
130
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
131
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
132
 
!            --- tableau de travail
133
 
!-------------------------------------------------------------------------------
134
 
!===============================================================================
135
 
 
136
 
implicit none
137
 
 
138
 
!===============================================================================
139
 
! Common blocks
140
 
!===============================================================================
141
 
 
142
 
include "paramx.h"
143
 
include "numvar.h"
144
 
include "entsor.h"
145
 
include "optcal.h"
146
 
include "cstphy.h"
147
 
include "cstnum.h"
148
 
include "pointe.h"
149
 
include "period.h"
150
 
include "parall.h"
151
 
 
152
 
!===============================================================================
153
 
 
154
 
! Arguments
155
 
 
156
 
integer          idbia0 , idbra0
157
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
158
 
integer          nfml   , nprfml
159
 
integer          nnod   , lndfac , lndfbr , ncelbr
160
 
integer          nvar   , nscal  , nphas
161
 
integer          nideve , nrdeve , nituse , nrtuse
162
 
integer          iphas  , ivar   , isou   , ipp
163
 
 
164
 
integer          ifacel(2,nfac) , ifabor(nfabor)
165
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
166
 
integer          iprfml(nfml,nprfml)
167
 
integer          ipnfac(nfac+1), nodfac(lndfac)
168
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
169
 
integer          idevel(nideve), ituser(nituse), ia(*)
170
 
 
171
 
double precision xyzcen(ndim,ncelet)
172
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
173
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
174
 
double precision xyznod(ndim,nnod), volume(ncelet)
175
 
double precision rtp(ncelet,*), rtpa(ncelet,*)
176
 
double precision propce(ncelet,*)
177
 
double precision propfa(nfac,*), propfb(nfabor,*)
178
 
double precision coefa(nfabor,*), coefb(nfabor,*)
179
 
double precision produc(6,ncelet)
180
 
double precision smbr(ncelet)
181
 
double precision coefax(nfabor), coefbx(nfabor)
182
 
double precision produk(ncelet),w2(ncelet),w3(ncelet)
183
 
double precision w4(ncelet),epsk(ncelet),w6(ncelet)
184
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
185
 
 
186
 
! Local variables
187
 
 
188
 
integer          idebia, idebra
189
 
integer          ifacpt, iel   , ii    , jj    , kk    , mm
190
 
integer          irkm  , irki  , irkj  , iskm  , iski  , iskj
191
 
integer          ir11ip, ir22ip, ir33ip, ir12ip, ir13ip, ir23ip
192
 
integer          ieiph , ipcrom, ipcroo
193
 
integer          ifac  , idimte, itenso
194
 
integer          inc   , iccocg, iphydp, ivar0 , iityph
195
 
 
196
 
double precision cmu075, distxn, d2s3  , trrij , xk
197
 
double precision unssur, vnk   , vnm   , vni   , vnj
198
 
double precision deltki, deltkj, deltkm, deltij
199
 
double precision aa    , bb    , surfbn, xnorme
200
 
 
201
 
 
202
 
!===============================================================================
203
 
 
204
 
!===============================================================================
205
 
! 1. INITIALISATION
206
 
!===============================================================================
207
 
 
208
 
idebia = idbia0
209
 
idebra = idbra0
210
 
 
211
 
ir11ip = ir11(iphas)
212
 
ir22ip = ir22(iphas)
213
 
ir33ip = ir33(iphas)
214
 
ir12ip = ir12(iphas)
215
 
ir13ip = ir13(iphas)
216
 
ir23ip = ir23(iphas)
217
 
ieiph  = iep (iphas)
218
 
 
219
 
ipcrom = ipproc(irom  (iphas))
220
 
ipcroo = ipcrom
221
 
if(isto2t(iphas).gt.0.and.iroext(iphas).gt.0) then
222
 
  ipcroo = ipproc(iroma(iphas))
223
 
endif
224
 
 
225
 
deltij = 1.0d0
226
 
if(isou.gt.3) then
227
 
  deltij = 0.0d0
228
 
endif
229
 
 
230
 
cmu075 = cmu**0.75d0
231
 
d2s3   = 2.d0/3.d0
232
 
 
233
 
!===============================================================================
234
 
! 2. CALCUL AUX CELLULES DES NORMALES EN PAROI CORRESPONDANTES
235
 
!===============================================================================
236
 
 
237
 
!     On stocke les composantes dans les tableaux de travail W2, W3 et W4
238
 
 
239
 
if(abs(icdpar).eq.2) then
240
 
 
241
 
!     On connait la face de paroi correspondante
242
 
 
243
 
    do iel = 1, ncel
244
 
      ifacpt = ia(iifapa(iphas)-1+iel)
245
 
      surfbn = ra(isrfbn-1+ifacpt)
246
 
      unssur = 1.d0/surfbn
247
 
      w2(iel)= surfbo(1,ifacpt)*unssur
248
 
      w3(iel)= surfbo(2,ifacpt)*unssur
249
 
      w4(iel)= surfbo(3,ifacpt)*unssur
250
 
    enddo
251
 
 
252
 
elseif(abs(icdpar).eq.1) then
253
 
 
254
 
!     La normale est definie comme - gradient de la distance
255
 
!       a la paroi
256
 
 
257
 
!       La distance a la paroi vaut 0 en paroi
258
 
!         par definition et obeit a un flux nul ailleurs
259
 
 
260
 
  iityph = iitypf+(iphas-1)*nfabor
261
 
  do ifac = 1, nfabor
262
 
    if(ia(iityph+ifac-1).eq.iparoi .or.                           &
263
 
       ia(iityph+ifac-1).eq.iparug) then
264
 
      coefax(ifac) = 0.0d0
265
 
      coefbx(ifac) = 0.0d0
266
 
    else
267
 
      coefax(ifac) = 0.0d0
268
 
      coefbx(ifac) = 1.0d0
269
 
    endif
270
 
  enddo
271
 
 
272
 
!       Calcul du gradient
273
 
 
274
 
  if(irangp.ge.0) call parcom (ra(idipar))
275
 
                  !==========
276
 
  if(iperio.eq.1) then
277
 
    idimte = 0
278
 
    itenso = 0
279
 
    call percom                                                   &
280
 
    !==========
281
 
   ( idimte , itenso ,                                            &
282
 
     ra(idipar), ra(idipar), ra(idipar),                          &
283
 
     ra(idipar), ra(idipar), ra(idipar),                          &
284
 
     ra(idipar), ra(idipar), ra(idipar))
285
 
  endif
286
 
 
287
 
  inc    = 1
288
 
  iccocg = 1
289
 
  iphydp = 0
290
 
  ivar0  = 0
291
 
 
292
 
  call grdcel                                                     &
293
 
  !==========
294
 
 ( idebia , idebra ,                                              &
295
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
296
 
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
297
 
   nideve , nrdeve , nituse , nrtuse ,                            &
298
 
   ivar0  , imrgra , inc    , iccocg , nswrgy , imligy , iphydp , &
299
 
   iwarny , nfecra ,                                              &
300
 
   epsrgy , climgy , extray ,                                     &
301
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
302
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
303
 
   idevel , ituser , ia     ,                                     &
304
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
305
 
   ra     , ra     , ra     ,                                     &
306
 
   ra(idipar) , coefax , coefbx ,                                 &
307
 
   w2     , w3     , w4     ,                                     &
308
 
!        ------   ------   ------
309
 
   produk , epsk   , w6     ,                                     &
310
 
   rdevel , rtuser , ra     )
311
 
 
312
 
 
313
 
!     Normalisation (attention, le gradient peut etre nul, parfois)
314
 
 
315
 
  do iel = 1 ,ncel
316
 
    xnorme = max(sqrt(w2(iel)**2+w3(iel)**2+w4(iel)**2),epzero)
317
 
    w2(iel) = -w2(iel)/xnorme
318
 
    w3(iel) = -w3(iel)/xnorme
319
 
    w4(iel) = -w4(iel)/xnorme
320
 
  enddo
321
 
 
322
 
endif
323
 
 
324
 
!===============================================================================
325
 
! 2. CALCUL DE VARIABLES DE TRAVAIL
326
 
!===============================================================================
327
 
 
328
 
 
329
 
 
330
 
! ---> Production et k
331
 
 
332
 
do iel = 1 , ncel
333
 
  produk(iel) = 0.5d0 * (produc(1,iel)  + produc(2,iel)  +        &
334
 
                         produc(3,iel))
335
 
  xk          = 0.5d0 * (rtpa(iel,ir11ip) + rtpa(iel,ir22ip) +    &
336
 
                         rtpa(iel,ir33ip))
337
 
  epsk(iel)   = rtpa(iel,ieiph)/xk
338
 
enddo
339
 
 
340
 
 
341
 
 
342
 
! ---> Indices des tensions
343
 
 
344
 
if     ((isou.eq.1).or.(isou.eq.4).or.(isou.eq.5)) then
345
 
  ii = 1
346
 
elseif ((isou.eq.2).or.(isou.eq.6)) then
347
 
  ii = 2
348
 
elseif  (isou.eq.3) then
349
 
  ii = 3
350
 
endif
351
 
 
352
 
if     ((isou.eq.3).or.(isou.eq.5).or.(isou.eq.6)) then
353
 
  jj = 3
354
 
elseif ((isou.eq.2).or.(isou.eq.4)) then
355
 
  jj = 2
356
 
elseif  (isou.eq.1) then
357
 
  jj = 1
358
 
endif
359
 
 
360
 
! ---> Boucle pour construction des termes sources
361
 
 
362
 
do iel = 1, ncel
363
 
  w6(iel) = 0.d0
364
 
enddo
365
 
 
366
 
do kk = 1, 3
367
 
 
368
 
! ---> Sommes sur m
369
 
 
370
 
  do mm = 1, 3
371
 
 
372
 
!   --> Delta km
373
 
 
374
 
    if(kk.eq.mm) then
375
 
      deltkm = 1.0d0
376
 
    else
377
 
      deltkm = 0.0d0
378
 
    endif
379
 
 
380
 
!  --> R km
381
 
 
382
 
    if     ((kk*mm).eq.1) then
383
 
      irkm = ir11ip
384
 
      iskm = 1
385
 
    elseif ((kk*mm).eq.4) then
386
 
      irkm = ir22ip
387
 
      iskm = 2
388
 
    elseif ((kk*mm).eq.9) then
389
 
      irkm = ir33ip
390
 
      iskm = 3
391
 
    elseif ((kk*mm).eq.2) then
392
 
      irkm = ir12ip
393
 
      iskm = 4
394
 
    elseif ((kk*mm).eq.3) then
395
 
      irkm = ir13ip
396
 
      iskm = 5
397
 
    elseif ((kk*mm).eq.6) then
398
 
      irkm = ir23ip
399
 
      iskm = 6
400
 
    endif
401
 
 
402
 
!  --> Termes en R km et Phi km
403
 
 
404
 
    do iel = 1, ncel
405
 
 
406
 
      if    (kk.eq.1) then
407
 
        vnk    = w2(iel)
408
 
      elseif(kk.eq.2) then
409
 
        vnk    = w3(iel)
410
 
      elseif(kk.eq.3) then
411
 
        vnk    = w4(iel)
412
 
      endif
413
 
 
414
 
      if    (mm.eq.1) then
415
 
        vnm    = w2(iel)
416
 
      elseif(mm.eq.2) then
417
 
        vnm    = w3(iel)
418
 
      elseif(mm.eq.3) then
419
 
        vnm    = w4(iel)
420
 
      endif
421
 
 
422
 
      w6(iel) = w6(iel) + vnk*vnm*deltij*(                        &
423
 
             crijp1*rtpa(iel,irkm)*epsk(iel)                      &
424
 
            -crijp2                                               &
425
 
             *crij2*(produc(iskm,iel)-d2s3*produk(iel)*deltkm) )
426
 
    enddo
427
 
 
428
 
  enddo
429
 
 
430
 
! ---> Fin des sommes sur m
431
 
 
432
 
 
433
 
!  --> R ki
434
 
 
435
 
  if     ((kk*ii).eq.1) then
436
 
    irki = ir11ip
437
 
    iski = 1
438
 
  elseif ((kk*ii).eq.4) then
439
 
    irki = ir22ip
440
 
    iski = 2
441
 
  elseif ((kk*ii).eq.9) then
442
 
    irki = ir33ip
443
 
    iski = 3
444
 
  elseif ((kk*ii).eq.2) then
445
 
    irki = ir12ip
446
 
    iski = 4
447
 
  elseif ((kk*ii).eq.3) then
448
 
    irki = ir13ip
449
 
    iski = 5
450
 
  elseif ((kk*ii).eq.6) then
451
 
    irki = ir23ip
452
 
    iski = 6
453
 
  endif
454
 
 
455
 
!  --> R kj
456
 
 
457
 
  if     ((kk*jj).eq.1) then
458
 
    irkj = ir11ip
459
 
    iskj = 1
460
 
  elseif ((kk*jj).eq.4) then
461
 
    irkj = ir22ip
462
 
    iskj = 2
463
 
  elseif ((kk*jj).eq.9) then
464
 
    irkj = ir33ip
465
 
    iskj = 3
466
 
  elseif ((kk*jj).eq.2) then
467
 
    irkj = ir12ip
468
 
    iskj = 4
469
 
  elseif ((kk*jj).eq.3) then
470
 
    irkj = ir13ip
471
 
    iskj = 5
472
 
  elseif ((kk*jj).eq.6) then
473
 
    irkj = ir23ip
474
 
    iskj = 6
475
 
  endif
476
 
 
477
 
!   --> Delta ki
478
 
 
479
 
  if (kk.eq.ii) then
480
 
    deltki = 1.d0
481
 
  else
482
 
    deltki = 0.d0
483
 
  endif
484
 
 
485
 
!   --> Delta kj
486
 
 
487
 
  if (kk.eq.jj) then
488
 
    deltkj = 1.d0
489
 
  else
490
 
    deltkj = 0.d0
491
 
  endif
492
 
 
493
 
  do iel = 1, ncel
494
 
 
495
 
      if    (kk.eq.1) then
496
 
        vnk    = w2(iel)
497
 
      elseif(kk.eq.2) then
498
 
        vnk    = w3(iel)
499
 
      elseif(kk.eq.3) then
500
 
        vnk    = w4(iel)
501
 
      endif
502
 
 
503
 
      if    (ii.eq.1) then
504
 
        vni    = w2(iel)
505
 
      elseif(ii.eq.2) then
506
 
        vni    = w3(iel)
507
 
      elseif(ii.eq.3) then
508
 
        vni    = w4(iel)
509
 
      endif
510
 
 
511
 
      if    (jj.eq.1) then
512
 
        vnj    = w2(iel)
513
 
      elseif(jj.eq.2) then
514
 
        vnj    = w3(iel)
515
 
      elseif(jj.eq.3) then
516
 
        vnj    = w4(iel)
517
 
      endif
518
 
 
519
 
    w6(iel) = w6(iel) + 1.5d0*vnk*(                               &
520
 
    -crijp1*(rtpa(iel,irki)*vnj+rtpa(iel,irkj)*vni)*epsk(iel)     &
521
 
    +crijp2                                                       &
522
 
     *crij2*((produc(iski,iel)-d2s3*produk(iel)*deltki)*vnj       &
523
 
            +(produc(iskj,iel)-d2s3*produk(iel)*deltkj)*vni) )
524
 
 
525
 
  enddo
526
 
 
527
 
enddo
528
 
 
529
 
 
530
 
! ---> Distance a la paroi et fonction d'amortissement : W3
531
 
!     Pour chaque mode de calcul : meme code, test
532
 
!       en dehors de la boucle
533
 
 
534
 
if(abs(icdpar).eq.2) then
535
 
  do iel = 1 , ncel
536
 
    ifacpt = ia(iifapa(iphas)-1+iel)
537
 
    distxn =                                                      &
538
 
          (cdgfbo(1,ifacpt)-xyzcen(1,iel))**2                     &
539
 
         +(cdgfbo(2,ifacpt)-xyzcen(2,iel))**2                     &
540
 
         +(cdgfbo(3,ifacpt)-xyzcen(3,iel))**2
541
 
    distxn = sqrt(distxn)
542
 
    trrij  = 0.5d0 * (rtpa(iel,ir11ip) + rtpa(iel,ir22ip) +       &
543
 
                         rtpa(iel,ir33ip))
544
 
    aa = 1.d0
545
 
    bb = cmu075*trrij**1.5d0/(xkappa*rtpa(iel,ieiph)*distxn)
546
 
    w3(iel) = min(aa, bb)
547
 
  enddo
548
 
else
549
 
  do iel = 1 , ncel
550
 
    distxn =  max(ra(idipar+iel-1),epzero)
551
 
    trrij  = 0.5d0 * (rtpa(iel,ir11ip) + rtpa(iel,ir22ip) +       &
552
 
                         rtpa(iel,ir33ip))
553
 
    aa = 1.d0
554
 
    bb = cmu075*trrij**1.5d0/(xkappa*rtpa(iel,ieiph)*distxn)
555
 
    w3(iel) = min(aa, bb)
556
 
  enddo
557
 
endif
558
 
 
559
 
 
560
 
! ---> Increment du terme source
561
 
 
562
 
do iel = 1, ncel
563
 
  smbr(iel) = smbr(iel)                                           &
564
 
            + propce(iel,ipcroo)*volume(iel)*w6(iel)*w3(iel)
565
 
enddo
566
 
 
567
 
 
568
 
return
569
 
 
570
 
end subroutine