← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/base/reseps.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/reseps.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 reseps &
29
 
!================
30
 
 
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 ,                            &
36
 
   iphas  , ivar   , isou   , ipp    ,                            &
37
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
38
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
39
 
   icepdc , icetsm , itpsmp ,                                     &
40
 
   idevel , ituser , ia     ,                                     &
41
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
42
 
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
43
 
   coefa  , coefb  , grdvit , produc , grarox , graroy , graroz , &
44
 
   ckupdc , smcelp , gamma  ,                                     &
45
 
   viscf  , viscb  ,                                              &
46
 
   tslagr ,                                                       &
47
 
   dam    , xam    , drtp   , smbr   , rovsdt ,                   &
48
 
   w1     , w2     , w3     , w4     ,                            &
49
 
   w5     , w6     , w7     , w8     , w9     ,                   &
50
 
   rdevel , rtuser , ra     )
51
 
 
52
 
!===============================================================================
53
 
! FONCTION :
54
 
! ----------
55
 
 
56
 
! RESOLUTION DES EQUATIONS CONVECTION DIFFUSION TERME SOURCE
57
 
!   POUR EPSILON DANS LE CAS DU MODELE Rij-epsilon
58
 
 
59
 
!   On a ici ISOU   = 7 (7 ieme variable du Rij-epsilon)
60
 
 
61
 
!-------------------------------------------------------------------------------
62
 
!ARGU                             ARGUMENTS
63
 
!__________________.____._____.________________________________________________.
64
 
! name             !type!mode ! role                                           !
65
 
!__________________!____!_____!________________________________________________!
66
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
67
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
68
 
! ndim             ! i  ! <-- ! spatial dimension                              !
69
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
70
 
! ncel             ! i  ! <-- ! number of cells                                !
71
 
! nfac             ! i  ! <-- ! number of interior faces                       !
72
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
73
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
74
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
75
 
! nnod             ! i  ! <-- ! number of vertices                             !
76
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
77
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
78
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
79
 
! nvar             ! i  ! <-- ! total number of variables                      !
80
 
! nscal            ! i  ! <-- ! total number of scalars                        !
81
 
! nphas            ! i  ! <-- ! number of phases                               !
82
 
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
83
 
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
84
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
85
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
86
 
! iphas            ! i  ! <-- ! phase number                                   !
87
 
! ivar             ! i  ! <-- ! variable number                                !
88
 
! isou             ! e  ! <-- ! numero de passage                              !
89
 
! ipp              ! e  ! <-- ! numero de variable pour sorties post           !
90
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
91
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
92
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
93
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
94
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
95
 
!  (nfml, nprfml)  !    !     !                                                !
96
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
97
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
98
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
99
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
100
 
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
101
 
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
102
 
! itpsmp           ! te ! <-- ! type de source de masse pour la                !
103
 
! (ncesmp)         !    !     !  variables (cf. ustsma)                        !
104
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
105
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
106
 
! ia(*)            ! ia ! --- ! main integer work array                        !
107
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
108
 
!  (ndim, ncelet)  !    !     !                                                !
109
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
110
 
!  (ndim, nfac)    !    !     !                                                !
111
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
112
 
!  (ndim, nfabor)  !    !     !                                                !
113
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
114
 
!  (ndim, nfac)    !    !     !                                                !
115
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
116
 
!  (ndim, nfabor)  !    !     !                                                !
117
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
118
 
!  (ndim, nnod)    !    !     !                                                !
119
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
120
 
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
121
 
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
122
 
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
123
 
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
124
 
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
125
 
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
126
 
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
127
 
!  (nfabor, *)     !    !     !                                                !
128
 
! grdvit           ! tr ! --- ! tableau de travail pour terme grad             !
129
 
!  (ncelet,3,3)    !    !     !    de vitesse     uniqt pour iturb=31          !
130
 
! produc           ! tr ! <-- ! tableau de travail pour production             !
131
 
!  (6,ncelet)      !    !     ! (sans rho volume) uniqt pour iturb=30          !
132
 
! grarox,y,z       ! tr ! <-- ! tableau de travail pour grad rom               !
133
 
!  (ncelet)        !    !     !                                                !
134
 
! ckupdc           ! tr ! <-- ! tableau de travail pour pdc                    !
135
 
!  (ncepdp,6)      !    !     !                                                !
136
 
! smcelp(ncesmp    ! tr ! <-- ! valeur de la variable associee a la            !
137
 
!                  !    !     !  source de masse                               !
138
 
! gamma(ncesmp)    ! tr ! <-- ! valeur du flux de masse                        !
139
 
! viscf(nfac)      ! tr ! --- ! visc*surface/dist aux faces internes           !
140
 
! viscb(nfabor     ! tr ! --- ! visc*surface/dist aux faces de bord            !
141
 
! tslagr           ! tr ! <-- ! terme de couplage retour du                    !
142
 
!  (ncelet,*)      !    !     !   lagrangien                                   !
143
 
! dam(ncelet       ! tr ! --- ! tableau de travail pour matrice                !
144
 
! xam(nfac,*)      ! tr ! --- ! tableau de travail pour matrice                !
145
 
! drtp(ncelet      ! tr ! --- ! tableau de travail pour increment              !
146
 
! smbr(ncelet      ! tr ! --- ! tableau de travail pour sec mem                !
147
 
! rovsdt(ncelet    ! tr ! --- ! tableau de travail pour terme instat           !
148
 
! w1...9(ncelet    ! tr ! --- ! tableau de travail                             !
149
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
150
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
151
 
! ra(*)            ! ra ! --- ! main real work array                           !
152
 
!__________________!____!_____!________________________________________________!
153
 
 
154
 
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
155
 
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
156
 
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
157
 
!            --- tableau de travail
158
 
!-------------------------------------------------------------------------------
159
 
!===============================================================================
160
 
 
161
 
implicit none
162
 
 
163
 
!===============================================================================
164
 
! Common blocks
165
 
!===============================================================================
166
 
 
167
 
include "dimfbr.h"
168
 
include "paramx.h"
169
 
include "numvar.h"
170
 
include "entsor.h"
171
 
include "optcal.h"
172
 
include "cstnum.h"
173
 
include "cstphy.h"
174
 
include "pointe.h"
175
 
include "period.h"
176
 
include "parall.h"
177
 
include "lagpar.h"
178
 
include "lagran.h"
179
 
 
180
 
!===============================================================================
181
 
 
182
 
! Arguments
183
 
 
184
 
integer          idbia0 , idbra0
185
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
186
 
integer          nfml   , nprfml
187
 
integer          nnod   , lndfac , lndfbr , ncelbr
188
 
integer          nvar   , nscal  , nphas
189
 
integer          ncepdp , ncesmp
190
 
integer          nideve , nrdeve , nituse , nrtuse
191
 
integer          iphas  , ivar   , isou   , ipp
192
 
 
193
 
integer          ifacel(2,nfac) , ifabor(nfabor)
194
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
195
 
integer          iprfml(nfml,nprfml)
196
 
integer          ipnfac(nfac+1), nodfac(lndfac)
197
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
198
 
integer          icepdc(ncepdp)
199
 
integer          icetsm(ncesmp), itpsmp(ncesmp)
200
 
integer          idevel(nideve), ituser(nituse)
201
 
integer          ia(*)
202
 
 
203
 
double precision xyzcen(ndim,ncelet)
204
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
205
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
206
 
double precision xyznod(ndim,nnod), volume(ncelet)
207
 
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
208
 
double precision propce(ncelet,*)
209
 
double precision propfa(nfac,*), propfb(ndimfb,*)
210
 
double precision coefa(ndimfb,*), coefb(ndimfb,*)
211
 
double precision produc(6,ncelet), grdvit(ncelet,3,3)
212
 
double precision grarox(ncelet), graroy(ncelet), graroz(ncelet)
213
 
double precision ckupdc(ncepdp,6)
214
 
double precision smcelp(ncesmp), gamma(ncesmp)
215
 
double precision viscf(nfac), viscb(nfabor)
216
 
double precision tslagr(ncelet,*)
217
 
double precision dam(ncelet), xam(nfac,2)
218
 
double precision drtp(ncelet), smbr(ncelet), rovsdt(ncelet)
219
 
double precision w1(ncelet), w2(ncelet), w3(ncelet)
220
 
double precision w4(ncelet), w5(ncelet), w6(ncelet)
221
 
double precision w7(ncelet), w8(ncelet), w9(ncelet)
222
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
223
 
 
224
 
! Local variables
225
 
 
226
 
integer          idebia, idebra, ifinia
227
 
integer          init  , ifac  , iel   , inc   , iccocg
228
 
integer          ii    ,jj     , iiun
229
 
integer          ir11ip, ir22ip, ir33ip, ir12ip, ir13ip, ir23ip
230
 
integer          ieiph , iuiph
231
 
integer          iclvar, iclvaf
232
 
integer          ipcrom, ipcroo, ipcvis, ipcvst, iflmas, iflmab
233
 
integer          nswrgp, imligp, iwarnp, iphydp
234
 
integer          iconvp, idiffp, ndircp, ireslp
235
 
integer          nitmap, nswrsp, ircflp, ischcp, isstpp, iescap
236
 
integer          imgrp , ncymxp, nitmfp
237
 
integer          idimte, itenso
238
 
integer          iptsta
239
 
integer          maxelt, ils
240
 
double precision blencp, epsilp, epsrgp, climgp, extrap, relaxp
241
 
double precision epsrsp
242
 
double precision trprod , trrij ,csteps, rctse
243
 
double precision grdpx,grdpy,grdpz,grdsn
244
 
double precision surfn2
245
 
double precision tseps , kseps , ceps2
246
 
double precision tuexpe, thets , thetv , thetap, thetp1
247
 
 
248
 
!===============================================================================
249
 
 
250
 
!===============================================================================
251
 
! 1. INITIALISATION
252
 
!===============================================================================
253
 
 
254
 
idebia = idbia0
255
 
idebra = idbra0
256
 
 
257
 
if(iwarni(ivar).ge.1) then
258
 
  write(nfecra,1000) nomvar(ipp)
259
 
endif
260
 
 
261
 
iuiph  = iu(iphas)
262
 
ir11ip = ir11(iphas)
263
 
ir22ip = ir22(iphas)
264
 
ir33ip = ir33(iphas)
265
 
ir12ip = ir12(iphas)
266
 
ir13ip = ir13(iphas)
267
 
ir23ip = ir23(iphas)
268
 
ieiph  = iep (iphas)
269
 
 
270
 
ipcrom = ipproc(irom  (iphas))
271
 
ipcvis = ipproc(iviscl(iphas))
272
 
ipcvst = ipproc(ivisct(iphas))
273
 
iflmas = ipprof(ifluma(iuiph))
274
 
iflmab = ipprob(ifluma(iuiph))
275
 
 
276
 
iclvar = iclrtp(ivar,icoef)
277
 
iclvaf = iclrtp(ivar,icoeff)
278
 
 
279
 
!     Constante Ce2, qui vaut CE2 pour ITURB=30 et CSSGE2 pour ITRUB=31
280
 
if (iturb(iphas).eq.30) then
281
 
  ceps2 = ce2
282
 
else
283
 
  ceps2 = cssge2
284
 
endif
285
 
 
286
 
!     S pour Source, V pour Variable
287
 
thets  = thetst(iphas)
288
 
thetv  = thetav(ivar )
289
 
 
290
 
ipcroo = ipcrom
291
 
if(isto2t(iphas).gt.0.and.iroext(iphas).gt.0) then
292
 
  ipcroo = ipproc(iroma(iphas))
293
 
endif
294
 
if(isto2t(iphas).gt.0) then
295
 
  iptsta = ipproc(itstua(iphas))
296
 
else
297
 
  iptsta = 0
298
 
endif
299
 
 
300
 
do iel = 1, ncel
301
 
  smbr(iel) = 0.d0
302
 
enddo
303
 
do iel = 1, ncel
304
 
  rovsdt(iel) = 0.d0
305
 
enddo
306
 
 
307
 
!===============================================================================
308
 
! 2. TERMES SOURCES  UTILISATEURS
309
 
!===============================================================================
310
 
 
311
 
maxelt = max(ncelet, nfac, nfabor)
312
 
ils    = idebia
313
 
ifinia = ils + maxelt
314
 
CALL IASIZE('RESEPS',IFINIA)
315
 
 
316
 
call ustsri                                                       &
317
 
!==========
318
 
 ( ifinia , idebra ,                                              &
319
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
320
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
321
 
   nvar   , nscal  , nphas  , ncepdp , ncesmp ,                   &
322
 
   nideve , nrdeve , nituse , nrtuse ,                            &
323
 
   iphas  , ivar   ,                                              &
324
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , ia(ils), &
325
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
326
 
   icepdc , icetsm , itpsmp ,                                     &
327
 
   idevel , ituser , ia     ,                                     &
328
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
329
 
   dt     , rtpa   , propce , propfa , propfb ,                   &
330
 
   coefa  , coefb  , ckupdc , smcelp , gamma  , grdvit , produc , &
331
 
   smbr   , rovsdt ,                                              &
332
 
!        ------   ------
333
 
   viscf  , viscb  , xam    ,                                     &
334
 
   w1     , w2     , w3     , w4     , w5     , w6     ,          &
335
 
   w7     , w8     , w9     , dam    , drtp   ,                   &
336
 
   rdevel , rtuser , ra     )
337
 
 
338
 
!     Si on extrapole les T.S.
339
 
if(isto2t(iphas).gt.0) then
340
 
  do iel = 1, ncel
341
 
!       Sauvegarde pour echange
342
 
    tuexpe = propce(iel,iptsta+isou-1)
343
 
!       Pour la suite et le pas de temps suivant
344
 
    propce(iel,iptsta+isou-1) = smbr(iel)
345
 
!       Second membre du pas de temps precedent
346
 
!       On suppose -ROVSDT > 0 : on implicite
347
 
!          le terme source utilisateur (le reste)
348
 
    smbr(iel) = rovsdt(iel)*rtpa(iel,ivar) - thets*tuexpe
349
 
!       Diagonale
350
 
    rovsdt(iel) = - thetv*rovsdt(iel)
351
 
  enddo
352
 
else
353
 
  do iel = 1, ncel
354
 
    smbr(iel)   = rovsdt(iel)*rtpa(iel,ivar) + smbr(iel)
355
 
    rovsdt(iel) = max(-rovsdt(iel),zero)
356
 
  enddo
357
 
endif
358
 
 
359
 
!===============================================================================
360
 
! 3.  TERMES SOURCES LAGRANGIEN : COUPLAGE RETOUR
361
 
!===============================================================================
362
 
 
363
 
!     Ordre 2 non pris en compte
364
 
 if (iilagr.eq.2 .and. ltsdyn.eq.1 .and. iphas.eq.ilphas) then
365
 
 
366
 
   do iel = 1,ncel
367
 
! Ts sur eps
368
 
     tseps = -0.5d0 * ( tslagr(iel,itsr11)                        &
369
 
                      + tslagr(iel,itsr22)                        &
370
 
                      + tslagr(iel,itsr33) )
371
 
! rapport k/eps
372
 
     kseps = 0.5d0 * ( rtpa(iel,ir11ip)                           &
373
 
                     + rtpa(iel,ir22ip)                           &
374
 
                     + rtpa(iel,ir33ip) )                         &
375
 
                     / rtpa(iel,ieiph)
376
 
 
377
 
     smbr(iel)   = smbr(iel) + ce4 *tseps *rtpa(iel,ieiph) /kseps
378
 
     rovsdt(iel) = rovsdt(iel) + max( (-ce4*tseps/kseps) , zero)
379
 
   enddo
380
 
 
381
 
 endif
382
 
 
383
 
!===============================================================================
384
 
! 4. TERME SOURCE DE MASSE
385
 
!===============================================================================
386
 
 
387
 
 
388
 
if (ncesmp.gt.0) then
389
 
 
390
 
!       Entier egal a 1 (pour navsto : nb de sur-iter)
391
 
  iiun = 1
392
 
 
393
 
!       On incremente SMBR par -Gamma RTPA et ROVSDT par Gamma (*theta)
394
 
  call catsma                                                     &
395
 
  !==========
396
 
 ( ncelet , ncel   , ncesmp , iiun   , isto2t(iphas) , thetv ,    &
397
 
   icetsm , itpsmp ,                                              &
398
 
   volume , rtpa(1,ivar) , smcelp , gamma  ,                      &
399
 
   smbr   ,  rovsdt , w1 )
400
 
 
401
 
!       Si on extrapole les TS on met Gamma Pinj dans PROPCE
402
 
  if(isto2t(iphas).gt.0) then
403
 
    do iel = 1, ncel
404
 
      propce(iel,iptsta+isou-1) =                                 &
405
 
      propce(iel,iptsta+isou-1) + w1(iel)
406
 
    enddo
407
 
!       Sinon on le met directement dans SMBR
408
 
  else
409
 
    do iel = 1, ncel
410
 
      smbr(iel) = smbr(iel) + w1(iel)
411
 
    enddo
412
 
  endif
413
 
 
414
 
endif
415
 
 
416
 
!===============================================================================
417
 
! 5. TERME D'ACCUMULATION DE MASSE -(dRO/dt)*VOLUME
418
 
!    ET TERME INSTATIONNAIRE
419
 
!===============================================================================
420
 
 
421
 
! ---> Calcul de mij
422
 
 
423
 
init = 1
424
 
call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,                  &
425
 
               ifacel,ifabor,propfa(1,iflmas),propfb(1,iflmab),w1)
426
 
 
427
 
! ---> Ajout au second membre
428
 
 
429
 
do iel = 1, ncel
430
 
  smbr(iel) = smbr(iel)                                           &
431
 
              + iconv(ivar)*w1(iel)*rtpa(iel,ivar)
432
 
enddo
433
 
 
434
 
! ---> Ajout dans la diagonale de la matrice
435
 
!     Extrapolation ou non, meme forme par coherence avec bilsc2
436
 
 
437
 
do iel = 1, ncel
438
 
  rovsdt(iel) = rovsdt(iel)                                       &
439
 
           + istat(ivar)*(propce(iel,ipcrom)/dt(iel))*volume(iel) &
440
 
           - iconv(ivar)*w1(iel)*thetv
441
 
enddo
442
 
 
443
 
 
444
 
!===============================================================================
445
 
! 6. PRODUCTION RHO * Ce1 * epsilon / k * P
446
 
!    DISSIPATION RHO*Ce2.epsilon/k*epsilon
447
 
!===============================================================================
448
 
 
449
 
 
450
 
! ---> Calcul de k pour la suite du sous-programme
451
 
!       on utilise un tableau de travail puisqu'il y en a...
452
 
do iel = 1, ncel
453
 
  w8(iel) = 0.5d0 * (rtpa(iel,ir11ip) + rtpa(iel,ir22ip) +        &
454
 
                     rtpa(iel,ir33ip))
455
 
enddo
456
 
! ---> Calcul de la trace de la production, suivant qu'on est en
457
 
!     Rij standard ou en SSG (utilisation de PRODUC ou GRDVIT)
458
 
if (iturb(iphas).eq.30) then
459
 
  do iel = 1, ncel
460
 
    w9(iel) = 0.5d0*(produc(1,iel)+produc(2,iel)+produc(3,iel))
461
 
  enddo
462
 
else
463
 
  do iel = 1, ncel
464
 
    w9(iel) = -( rtpa(iel,ir11ip)*grdvit(iel,1,1) +               &
465
 
                 rtpa(iel,ir12ip)*grdvit(iel,1,2) +               &
466
 
                 rtpa(iel,ir13ip)*grdvit(iel,1,3) +               &
467
 
                 rtpa(iel,ir12ip)*grdvit(iel,2,1) +               &
468
 
                 rtpa(iel,ir22ip)*grdvit(iel,2,2) +               &
469
 
                 rtpa(iel,ir23ip)*grdvit(iel,2,3) +               &
470
 
                 rtpa(iel,ir13ip)*grdvit(iel,3,1) +               &
471
 
                 rtpa(iel,ir23ip)*grdvit(iel,3,2) +               &
472
 
                 rtpa(iel,ir33ip)*grdvit(iel,3,3) )
473
 
  enddo
474
 
endif
475
 
 
476
 
 
477
 
!     Terme explicite (Production)
478
 
 
479
 
do iel = 1, ncel
480
 
!     Demi-traces
481
 
  trprod = w9(iel)
482
 
  trrij  = w8(iel)
483
 
  w1(iel)   =             propce(iel,ipcroo)*volume(iel)*         &
484
 
       ce1*rtpa(iel,ieiph)/trrij*trprod
485
 
enddo
486
 
 
487
 
!     Si on extrapole les T.S : PROPCE
488
 
if(isto2t(iphas).gt.0) then
489
 
  do iel = 1, ncel
490
 
    propce(iel,iptsta+isou-1) =                                   &
491
 
    propce(iel,iptsta+isou-1) + w1(iel)
492
 
  enddo
493
 
!     Sinon : SMBR
494
 
else
495
 
  do iel = 1, ncel
496
 
    smbr(iel) = smbr(iel) + w1(iel)
497
 
  enddo
498
 
endif
499
 
 
500
 
!     Terme implicite (Dissipation)
501
 
do iel = 1, ncel
502
 
  trrij  = w8(iel)
503
 
  smbr(iel) = smbr(iel) - propce(iel,ipcrom)*volume(iel)*         &
504
 
                         ceps2*rtpa(iel,ieiph)**2/trrij
505
 
enddo
506
 
 
507
 
! ---> Matrice
508
 
 
509
 
if(isto2t(iphas).gt.0) then
510
 
  thetap = thetv
511
 
else
512
 
  thetap = 1.d0
513
 
endif
514
 
do iel = 1, ncel
515
 
  trrij  = w8(iel)
516
 
  rovsdt(iel) = rovsdt(iel) + ceps2*propce(iel,ipcrom)*volume(iel)&
517
 
                     *rtpa(iel,ieiph)/trrij*thetap
518
 
enddo
519
 
 
520
 
!===============================================================================
521
 
! 7. TERMES DE GRAVITE
522
 
!===============================================================================
523
 
 
524
 
if(igrari(iphas).eq.1) then
525
 
 
526
 
  do iel = 1, ncel
527
 
    w7(iel) = 0.d0
528
 
  enddo
529
 
 
530
 
  call rijthe                                                     &
531
 
  !==========
532
 
 ( idebia , idebra ,                                              &
533
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
534
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
535
 
   nvar   , nscal  , nphas  ,                                     &
536
 
   nideve , nrdeve , nituse , nrtuse ,                            &
537
 
   iphas  , ivar   , isou   , ipp    ,                            &
538
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
539
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
540
 
   idevel , ituser , ia     ,                                     &
541
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
542
 
   rtp    , rtpa   , propce , propfa , propfb ,                   &
543
 
   coefa  , coefb  , grarox , graroy , graroz , w7     ,          &
544
 
   w1     , w2     , w3     , w4     , w5     , w6     ,          &
545
 
   rdevel , rtuser , ra     )
546
 
 
547
 
!     Si on extrapole les T.S. : PROPCE
548
 
if(isto2t(iphas).gt.0) then
549
 
  do iel = 1, ncel
550
 
     propce(iel,iptsta+isou-1) =                                  &
551
 
     propce(iel,iptsta+isou-1) + w7(iel)
552
 
   enddo
553
 
!     Sinon SMBR
554
 
 else
555
 
   do iel = 1, ncel
556
 
     smbr(iel) = smbr(iel) + w7(iel)
557
 
   enddo
558
 
 endif
559
 
 
560
 
endif
561
 
 
562
 
 
563
 
!===============================================================================
564
 
! 8.a TERMES DE DIFFUSION  A.grad(Eps) : PARTIE EXTRADIAGONALE EXPLICITE
565
 
!     RIJ STANDARD
566
 
!===============================================================================
567
 
 
568
 
if (iturb(iphas).eq.30) then
569
 
 
570
 
! ---> Calcul du grad(Eps)
571
 
 
572
 
 
573
 
  iccocg = 1
574
 
  inc = 1
575
 
 
576
 
  nswrgp = nswrgr(ivar )
577
 
  imligp = imligr(ivar )
578
 
  iwarnp = iwarni(ivar )
579
 
  epsrgp = epsrgr(ivar )
580
 
  climgp = climgr(ivar )
581
 
  extrap = extrag(ivar )
582
 
  iphydp = 0
583
 
 
584
 
  call grdcel                                                     &
585
 
  !==========
586
 
 ( idebia , idebra ,                                              &
587
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
588
 
   nnod   , lndfac , lndfbr , ncelbr , nphas  ,                   &
589
 
   nideve , nrdeve , nituse , nrtuse ,                            &
590
 
   ivar   , imrgra , inc    , iccocg , nswrgp , imligp , iphydp , &
591
 
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
592
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
593
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
594
 
   idevel , ituser , ia     ,                                     &
595
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
596
 
   w1     , w1     , w1     ,                                     &
597
 
   rtpa(1,ivar )   , coefa(1,iclvar) , coefb(1,iclvar) ,          &
598
 
   w1     , w2     , w3     ,                                     &
599
 
!        ------   ------   ------
600
 
   w4     , w5     , w6     ,                                     &
601
 
   rdevel , rtuser , ra     )
602
 
 
603
 
! ---> Calcul des termes extradiagonaux de A.grad(Eps)
604
 
 
605
 
  do iel = 1 , ncel
606
 
    trrij  = w8(iel)
607
 
    csteps  = propce(iel,ipcroo) * crijep *trrij / rtpa(iel,ieiph)
608
 
    w4(iel) = csteps * ( rtpa(iel,ir12ip) * w2(iel) +             &
609
 
         rtpa(iel,ir13ip) * w3(iel) )
610
 
    w5(iel) = csteps * ( rtpa(iel,ir12ip) * w1(iel) +             &
611
 
         rtpa(iel,ir23ip) * w3(iel) )
612
 
    w6(iel) = csteps * ( rtpa(iel,ir13ip) * w1(iel) +             &
613
 
         rtpa(iel,ir23ip) * w2(iel) )
614
 
  enddo
615
 
 
616
 
! ---> Assemblage de { A.grad(Eps) } .S aux faces
617
 
 
618
 
  call vectds                                                     &
619
 
  !==========
620
 
 (ndim   , ncelet , ncel   , nfac   , nfabor   ,                  &
621
 
  ifacel , ifabor , ia     ,                                      &
622
 
  surfac , surfbo , ra(ipond) ,                                   &
623
 
  w4     , w5     , w6     ,                                      &
624
 
  viscf  , viscb  , ra     )
625
 
 
626
 
  init = 1
627
 
  call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,                &
628
 
                                   ifacel,ifabor,viscf,viscb,w4)
629
 
 
630
 
!     Si on extrapole les termes sources
631
 
  if(isto2t(iphas).gt.0) then
632
 
    do iel = 1, ncel
633
 
      propce(iel,iptsta+isou-1) =                                 &
634
 
           propce(iel,iptsta+isou-1) + w4(iel)
635
 
    enddo
636
 
!     Sinon
637
 
  else
638
 
    do iel = 1, ncel
639
 
      smbr(iel) = smbr(iel) + w4(iel)
640
 
    enddo
641
 
  endif
642
 
 
643
 
 
644
 
!===============================================================================
645
 
! 8.b TERMES DE DIFFUSION  A.grad(Eps) : PARTIE DIAGONALE
646
 
!     RIJ STANDARD
647
 
!===============================================================================
648
 
!     Implicitation de (grad(eps).n)n en gradient facette
649
 
!     Si IDIFRE=1, terme correctif explicite
650
 
!        grad(eps)-(grad(eps).n)n calcule en gradient cellule
651
 
!     Les termes de bord sont uniquement pris en compte dans la partie
652
 
!        en (grad(eps).n)n
653
 
!     (W1,W2,W3) contient toujours le gradient de la variable traitee
654
 
 
655
 
!     La parcom/percom-isation du gradient de epsilon a ete faite dans
656
 
!       grdcel. Pas utile de recommencer.
657
 
 
658
 
  if (idifre(iphas).eq.1) then
659
 
 
660
 
    do iel = 1, ncel
661
 
      trrij  = w8(iel)
662
 
      csteps = propce(iel,ipcroo) * crijep *trrij/rtpa(iel,ieiph)
663
 
      w4(iel)=csteps*rtpa(iel,ir11ip)
664
 
      w5(iel)=csteps*rtpa(iel,ir22ip)
665
 
      w6(iel)=csteps*rtpa(iel,ir33ip)
666
 
    enddo
667
 
 
668
 
! --->  TRAITEMENT DU PARALLELISME (MEMES DOUTES QUE PERIODICITE ?)
669
 
 
670
 
    if(irangp.ge.0) then
671
 
      call parcom (w4)
672
 
      !==========
673
 
      call parcom (w5)
674
 
      !==========
675
 
      call parcom (w6)
676
 
      !==========
677
 
    endif
678
 
 
679
 
! -->   TRAITEMENT DE LA PERIODICITE
680
 
    if(iperio.eq.1) then
681
 
 
682
 
      idimte = 21
683
 
      itenso = 0
684
 
 
685
 
      call percom                                                 &
686
 
      !==========
687
 
    ( idimte , itenso ,                                           &
688
 
      w4     , w4     , w4    ,                                   &
689
 
      w5     , w5     , w5    ,                                   &
690
 
      w6     , w6     , w6    )
691
 
 
692
 
    endif
693
 
 
694
 
    do ifac = 1, nfac
695
 
 
696
 
      ii=ifacel(1,ifac)
697
 
      jj=ifacel(2,ifac)
698
 
 
699
 
      surfn2 = ra(isrfan-1+ifac)**2
700
 
 
701
 
      grdpx=0.5d0*(w1(ii)+w1(jj))
702
 
      grdpy=0.5d0*(w2(ii)+w2(jj))
703
 
      grdpz=0.5d0*(w3(ii)+w3(jj))
704
 
      grdsn=grdpx*surfac(1,ifac)+grdpy*surfac(2,ifac)             &
705
 
           +grdpz*surfac(3,ifac)
706
 
      grdpx=grdpx-grdsn*surfac(1,ifac)/surfn2
707
 
      grdpy=grdpy-grdsn*surfac(2,ifac)/surfn2
708
 
      grdpz=grdpz-grdsn*surfac(3,ifac)/surfn2
709
 
 
710
 
      viscf(ifac)= 0.5d0*(                                        &
711
 
           (w4(ii)+w4(jj))*grdpx*surfac(1,ifac)                   &
712
 
           +(w5(ii)+w5(jj))*grdpy*surfac(2,ifac)                  &
713
 
           +(w6(ii)+w6(jj))*grdpz*surfac(3,ifac))
714
 
 
715
 
    enddo
716
 
 
717
 
    do ifac = 1, nfabor
718
 
      viscb(ifac) = 0.d0
719
 
    enddo
720
 
 
721
 
    init = 1
722
 
    call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,              &
723
 
       ifacel,ifabor,viscf,viscb,w1)
724
 
 
725
 
!     Si on extrapole les termes sources
726
 
    if(isto2t(iphas).gt.0) then
727
 
      do iel = 1, ncel
728
 
        propce(iel,iptsta+isou-1) =                               &
729
 
             propce(iel,iptsta+isou-1) + w1(iel)
730
 
      enddo
731
 
!     Sinon
732
 
    else
733
 
      do iel = 1, ncel
734
 
        smbr(iel) = smbr(iel) + w1(iel)
735
 
      enddo
736
 
    endif
737
 
 
738
 
  endif
739
 
 
740
 
! ---> Viscosite orthotrope pour partie implicite
741
 
 
742
 
  if( idiff(ivar).ge. 1 ) then
743
 
    do iel = 1, ncel
744
 
      trrij  = w8(iel)
745
 
      rctse = propce(iel,ipcrom) * crijep * trrij/rtpa(iel,ieiph)
746
 
      w1(iel) = propce(iel,ipcvis)                                &
747
 
           + idifft(ivar)*rctse*rtpa(iel,ir11ip)
748
 
      w2(iel) = propce(iel,ipcvis)                                &
749
 
           + idifft(ivar)*rctse*rtpa(iel,ir22ip)
750
 
      w3(iel) = propce(iel,ipcvis)                                &
751
 
           + idifft(ivar)*rctse*rtpa(iel,ir33ip)
752
 
    enddo
753
 
 
754
 
    call visort                                                   &
755
 
    !==========
756
 
 ( idebia , idebra ,                                              &
757
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
758
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
759
 
   nideve , nrdeve , nituse , nrtuse , imvisf ,                   &
760
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
761
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
762
 
   idevel , ituser , ia     ,                                     &
763
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
764
 
   w1     , w2     , w3     ,                                     &
765
 
   viscf  , viscb  ,                                              &
766
 
   rdevel , rtuser , ra     )
767
 
 
768
 
  else
769
 
 
770
 
    do ifac = 1, nfac
771
 
      viscf(ifac) = 0.d0
772
 
    enddo
773
 
    do ifac = 1, nfabor
774
 
      viscb(ifac) = 0.d0
775
 
    enddo
776
 
 
777
 
  endif
778
 
 
779
 
else
780
 
!===============================================================================
781
 
! 8.c TERMES DE DIFFUSION
782
 
!     RIJ SSG
783
 
!===============================================================================
784
 
! ---> Viscosite
785
 
 
786
 
  if( idiff(ivar).ge. 1 ) then
787
 
    do iel = 1, ncel
788
 
      w1(iel) = propce(iel,ipcvis)                                &
789
 
           + idifft(ivar)*propce(iel,ipcvst)/sigmae
790
 
    enddo
791
 
 
792
 
    call viscfa                                                   &
793
 
   !==========
794
 
 ( idebia , idebra ,                                              &
795
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
796
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
797
 
   nideve , nrdeve , nituse , nrtuse , imvisf ,                   &
798
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
799
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
800
 
   idevel , ituser , ia     ,                                     &
801
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
802
 
   w1     ,                                                       &
803
 
   viscf  , viscb  ,                                              &
804
 
   rdevel , rtuser , ra     )
805
 
 
806
 
  else
807
 
 
808
 
    do ifac = 1, nfac
809
 
      viscf(ifac) = 0.d0
810
 
    enddo
811
 
    do ifac = 1, nfabor
812
 
      viscb(ifac) = 0.d0
813
 
    enddo
814
 
 
815
 
  endif
816
 
 
817
 
endif
818
 
 
819
 
 
820
 
!===============================================================================
821
 
! 9. RESOLUTION
822
 
!===============================================================================
823
 
 
824
 
if(isto2t(iphas).gt.0) then
825
 
  thetp1 = 1.d0 + thets
826
 
  do iel = 1, ncel
827
 
    smbr(iel) = smbr(iel) + thetp1*propce(iel,iptsta+isou-1)
828
 
  enddo
829
 
endif
830
 
 
831
 
iconvp = iconv (ivar)
832
 
idiffp = idiff (ivar)
833
 
ndircp = ndircl(ivar)
834
 
ireslp = iresol(ivar)
835
 
nitmap = nitmax(ivar)
836
 
nswrsp = nswrsm(ivar)
837
 
nswrgp = nswrgr(ivar)
838
 
imligp = imligr(ivar)
839
 
ircflp = ircflu(ivar)
840
 
ischcp = ischcv(ivar)
841
 
isstpp = isstpc(ivar)
842
 
iescap = 0
843
 
imgrp  = imgr  (ivar)
844
 
ncymxp = ncymax(ivar)
845
 
nitmfp = nitmgf(ivar)
846
 
!MO      IPP    =
847
 
iwarnp = iwarni(ivar)
848
 
blencp = blencv(ivar)
849
 
epsilp = epsilo(ivar)
850
 
epsrsp = epsrsm(ivar)
851
 
epsrgp = epsrgr(ivar)
852
 
climgp = climgr(ivar)
853
 
extrap = extrag(ivar)
854
 
relaxp = relaxv(ivar)
855
 
 
856
 
call codits                                                       &
857
 
!==========
858
 
 ( idebia , idebra ,                                              &
859
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
860
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
861
 
   nvar   , nscal  , nphas  ,                                     &
862
 
   nideve , nrdeve , nituse , nrtuse ,                            &
863
 
   idtvar , ivar   , iconvp , idiffp , ireslp , ndircp , nitmap , &
864
 
   imrgra , nswrsp , nswrgp , imligp , ircflp ,                   &
865
 
   ischcp , isstpp , iescap ,                                     &
866
 
   imgrp  , ncymxp , nitmfp , ipp    , iwarnp ,                   &
867
 
   blencp , epsilp , epsrsp , epsrgp , climgp , extrap ,          &
868
 
   relaxp , thetv  ,                                              &
869
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml ,                   &
870
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
871
 
   idevel , ituser , ia     ,                                     &
872
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
873
 
   rtpa(1,ivar)    , rtpa(1,ivar)    ,                            &
874
 
                     coefa(1,iclvar) , coefb(1,iclvar) ,          &
875
 
                     coefa(1,iclvaf) , coefb(1,iclvaf) ,          &
876
 
                     propfa(1,iflmas), propfb(1,iflmab),          &
877
 
   viscf  , viscb  , viscf  , viscb  ,                            &
878
 
   rovsdt , smbr   , rtp(1,ivar)     ,                            &
879
 
   dam    , xam    , drtp   ,                                     &
880
 
   w1     , w2     , w3     , w4     , w5     ,                   &
881
 
   w6     , w7     , w8     , w9     ,                            &
882
 
   rdevel , rtuser , ra     )
883
 
 
884
 
!===============================================================================
885
 
! 10. IMPRESSIONS
886
 
!===============================================================================
887
 
 
888
 
 
889
 
!--------
890
 
! FORMATS
891
 
!--------
892
 
 
893
 
#if defined(_CS_LANG_FR)
894
 
 
895
 
 1000 format(/,'           RESOLUTION POUR LA VARIABLE ',A8,/)
896
 
 
897
 
#else
898
 
 
899
 
 1000 format(/,'           SOLVING VARIABLE ',A8           ,/)
900
 
 
901
 
#endif
902
 
 
903
 
!12345678 : MAX: 12345678901234 MIN: 12345678901234 NORM: 12345678901234
904
 
!----
905
 
! FIN
906
 
!----
907
 
 
908
 
return
909
 
 
910
 
end subroutine