← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/turb/turbkw.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/turb/turbkw.f90
 
1
!-------------------------------------------------------------------------------
 
2
 
 
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.
 
20
 
 
21
!-------------------------------------------------------------------------------
 
22
 
 
23
subroutine turbkw &
 
24
!================
 
25
 
 
26
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
 
27
   icepdc , icetsm , itypsm ,                                     &
 
28
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
 
29
   tslagr , coefa  , coefb  , ckupdc , smacel )
 
30
 
 
31
!===============================================================================
 
32
! FONCTION :
 
33
! ----------
 
34
 
 
35
! RESOLUTION DES EQUATIONS K-OMEGA SST 1 PHASE INCOMPRESSIBLE OU
 
36
! RHO VARIABLE SUR UN PAS DE TEMPS
 
37
 
 
38
!-------------------------------------------------------------------------------
 
39
!ARGU                             ARGUMENTS
 
40
!__________________.____._____.________________________________________________.
 
41
! name             !type!mode ! role                                           !
 
42
!__________________!____!_____!________________________________________________!
 
43
! nvar             ! i  ! <-- ! total number of variables                      !
 
44
! nscal            ! i  ! <-- ! total number of scalars                        !
 
45
! ncepdp           ! i  ! <-- ! number of cells with head loss                 !
 
46
! ncesmp           ! i  ! <-- ! number of cells with mass source term          !
 
47
! icepdc(ncelet    ! te ! <-- ! numero des ncepdp cellules avec pdc            !
 
48
! icetsm(ncesmp    ! te ! <-- ! numero des cellules a source de masse          !
 
49
! itypsm           ! te ! <-- ! type de source de masse pour les               !
 
50
! (ncesmp,nvar)    !    !     !  variables (cf. ustsma)                        !
 
51
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
 
52
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
 
53
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
 
54
! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers            !
 
55
! propfa(nfac, *)  ! ra ! <-- ! physical properties at interior face centers   !
 
56
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
 
57
! tslagr           ! tr ! <-- ! terme de couplage retour du                    !
 
58
!(ncelet,*)        !    !     !     lagrangien                                 !
 
59
! coefa, coefb     ! ra ! <-- ! boundary conditions                            !
 
60
!  (nfabor, *)     !    !     !                                                !
 
61
! ckupdc           ! tr ! <-- ! tableau de travail pour pdc                    !
 
62
!  (ncepdp,6)      !    !     !                                                !
 
63
! smacel           ! tr ! <-- ! valeur des variables associee a la             !
 
64
! (ncesmp,*   )    !    !     !  source de masse                               !
 
65
!                  !    !     !  pour ivar=ipr, smacel=flux de masse           !
 
66
!__________________!____!_____!________________________________________________!
 
67
 
 
68
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
 
69
!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
 
70
!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
 
71
!            --- tableau de travail
 
72
!===============================================================================
 
73
 
 
74
!===============================================================================
 
75
! Module files
 
76
!===============================================================================
 
77
 
 
78
use paramx
 
79
use dimens, only: ndimfb
 
80
use numvar
 
81
use entsor
 
82
use cstnum
 
83
use cstphy
 
84
use optcal
 
85
use lagran
 
86
use pointe, only: s2kw, divukw, ifapat, dispar
 
87
use parall
 
88
use mesh
 
89
 
 
90
!===============================================================================
 
91
 
 
92
implicit none
 
93
 
 
94
! Arguments
 
95
 
 
96
integer          nvar   , nscal
 
97
integer          ncepdp , ncesmp
 
98
 
 
99
integer          icepdc(ncepdp)
 
100
integer          icetsm(ncesmp), itypsm(ncesmp,nvar)
 
101
 
 
102
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
 
103
double precision propce(ncelet,*)
 
104
double precision propfa(nfac,*), propfb(ndimfb,*)
 
105
double precision tslagr(ncelet,*)
 
106
double precision coefa(ndimfb,*), coefb(ndimfb,*)
 
107
double precision ckupdc(ncepdp,6), smacel(ncesmp,nvar)
 
108
 
 
109
! Local variables
 
110
 
 
111
character*80     chaine
 
112
integer          iel   , ifac  , init  , inc   , iccocg, ivar
 
113
integer          ii, iivar , iiun  , ifacpt
 
114
integer          iclipk, iclipw, isqrt
 
115
integer          nswrgp, imligp
 
116
integer          iclikp, iclomg
 
117
integer          iclvar, iclvaf
 
118
integer          iconvp, idiffp, ndircp, ireslp
 
119
integer          nitmap, nswrsp, ircflp, ischcp, isstpp, iescap
 
120
integer          imgrp , ncymxp, nitmfp
 
121
integer          ipcrom, ipbrom, ipcvst, ipcvis, iflmas, iflmab
 
122
integer          iwarnp, ipp
 
123
integer          iptsta
 
124
integer          ipcroo, ipbroo, ipcvto, ipcvlo
 
125
 
 
126
double precision rnorm , d2s3, divp23, epz2
 
127
double precision deltk , deltw, a11, a12, a22, a21
 
128
double precision unsdet, romvsd
 
129
double precision prdtur, xk, xw, xeps, xnu
 
130
double precision visct , rom
 
131
double precision blencp, epsilp, epsrgp, climgp, extrap, relaxp
 
132
double precision thetp1, thetak, thetaw, thets, thetap, epsrsp
 
133
double precision tuexpk, tuexpw
 
134
double precision cdkw, xarg1, xxf1, xgamma, xbeta, sigma, produc
 
135
double precision var, vrmin, vrmax
 
136
 
 
137
double precision rvoid(1)
 
138
 
 
139
double precision, allocatable, dimension(:) :: viscf, viscb
 
140
double precision, allocatable, dimension(:) :: dam
 
141
double precision, allocatable, dimension(:) :: smbrk, smbrw, rovsdt
 
142
double precision, allocatable, dimension(:) :: tinstk, tinstw, xf1
 
143
double precision, allocatable, dimension(:,:) :: gradk, grado, grad
 
144
double precision, allocatable, dimension(:) :: w1, w2, w3
 
145
double precision, allocatable, dimension(:) :: w4, w5, w6
 
146
double precision, allocatable, dimension(:) :: w7, w8
 
147
 
 
148
!===============================================================================
 
149
 
 
150
!===============================================================================
 
151
! 1. INITIALISATION
 
152
!===============================================================================
 
153
 
 
154
! Allocate temporary arrays for the turbulence resolution
 
155
allocate(viscf(nfac), viscb(nfabor))
 
156
allocate(dam(ncelet))
 
157
allocate(smbrk(ncelet), smbrw(ncelet), rovsdt(ncelet))
 
158
allocate(tinstk(ncelet), tinstw(ncelet), xf1(ncelet))
 
159
 
 
160
! Allocate work arrays
 
161
allocate(w1(ncelet), w2(ncelet), w3(ncelet))
 
162
allocate(w4(ncelet), w5(ncelet), w6(ncelet))
 
163
allocate(w7(ncelet), w8(ncelet))
 
164
 
 
165
 
 
166
epz2 = epzero**2
 
167
 
 
168
iclikp = iclrtp(ik ,icoef)
 
169
iclomg = iclrtp(iomg,icoef)
 
170
 
 
171
ipcrom = ipproc(irom  )
 
172
ipcvst = ipproc(ivisct)
 
173
ipcvis = ipproc(iviscl)
 
174
iflmas = ipprof(ifluma(ik))
 
175
iflmab = ipprob(ifluma(ik))
 
176
ipbrom = ipprob(irom  )
 
177
 
 
178
thets  = thetst
 
179
 
 
180
ipcroo = ipcrom
 
181
ipbroo = ipbrom
 
182
ipcvto = ipcvst
 
183
ipcvlo = ipcvis
 
184
if(isto2t.gt.0) then
 
185
  if (iroext.gt.0) then
 
186
    ipcroo = ipproc(iroma)
 
187
    ipbroo = ipprob(iroma)
 
188
  endif
 
189
  if(iviext.gt.0) then
 
190
    ipcvto = ipproc(ivista)
 
191
    ipcvlo = ipproc(ivisla)
 
192
  endif
 
193
endif
 
194
 
 
195
if(isto2t.gt.0) then
 
196
  iptsta = ipproc(itstua)
 
197
else
 
198
  iptsta = 0
 
199
endif
 
200
 
 
201
if(iwarni(ik).ge.1) then
 
202
  write(nfecra,1000)
 
203
endif
 
204
 
 
205
 
 
206
!===============================================================================
 
207
! 2. CALCUL DE dk/dxj.dw/dxj
 
208
!      Le terme est stocke dans         W1
 
209
!      En sortie de l'etape on conserve W1
 
210
!===============================================================================
 
211
 
 
212
! Allocate temporary arrays for gradients calculation
 
213
allocate(gradk(ncelet,3), grado(ncelet,3))
 
214
 
 
215
iccocg = 1
 
216
inc = 1
 
217
 
 
218
nswrgp = nswrgr(ik)
 
219
imligp = imligr(ik)
 
220
iwarnp = iwarni(ik)
 
221
epsrgp = epsrgr(ik)
 
222
climgp = climgr(ik)
 
223
extrap = extrag(ik)
 
224
 
 
225
call grdcel &
 
226
!==========
 
227
 ( ik  , imrgra , inc    , iccocg , nswrgp , imligp ,             &
 
228
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
 
229
   rtpa(1,ik)   , coefa(1,iclikp) , coefb(1,iclikp) ,             &
 
230
   gradk  )
 
231
 
 
232
nswrgp = nswrgr(iomg)
 
233
imligp = imligr(iomg)
 
234
iwarnp = iwarni(iomg)
 
235
epsrgp = epsrgr(iomg)
 
236
climgp = climgr(iomg)
 
237
extrap = extrag(iomg)
 
238
 
 
239
call grdcel &
 
240
!==========
 
241
 ( iomg , imrgra , inc    , iccocg , nswrgp , imligp ,            &
 
242
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
 
243
   rtpa(1,iomg)  , coefa(1,iclomg) , coefb(1,iclomg) ,            &
 
244
   grado  )
 
245
 
 
246
do iel = 1, ncel
 
247
  w1(iel) = gradk(iel,1)*grado(iel,1) &
 
248
          + gradk(iel,2)*grado(iel,2) &
 
249
          + gradk(iel,3)*grado(iel,3)
 
250
enddo
 
251
 
 
252
! Free memory
 
253
deallocate(gradk, grado)
 
254
 
 
255
!====================================================
 
256
! 3. CALCUL DU COEFFICIENT DE PONDERATION F1
 
257
!      Le terme est stocke dans         XF1
 
258
!      En sortie de l'etape on conserve W1,XF1
 
259
!====================================================
 
260
 
 
261
 
 
262
if(abs(icdpar).eq.2) then
 
263
  do iel = 1, ncel
 
264
    ifacpt = ifapat(iel)
 
265
    w2(iel) = (cdgfbo(1,ifacpt)-xyzcen(1,iel))**2 &
 
266
            + (cdgfbo(2,ifacpt)-xyzcen(2,iel))**2 &
 
267
            + (cdgfbo(3,ifacpt)-xyzcen(3,iel))**2
 
268
    w2(iel) = sqrt(w2(iel))
 
269
  enddo
 
270
else
 
271
  do iel = 1, ncel
 
272
    w2(iel) =  max(dispar(iel),epzero)
 
273
  enddo
 
274
endif
 
275
 
 
276
!     En cas d'ordre 2 on utilise les valeurs en n car le terme en (1-F1)*W1
 
277
!     sera dans PROPCE. Du coup, on aura quand meme certaines "constantes"
 
278
!     intervenant dans des termes en n+1/2 (ex sigma_k pour la diffusion) calcules
 
279
!     a partir de F1 en n -> mais l'effet sur les "constantes" est faible
 
280
!     -> a garder en tete si on fait vraiment de l'ordre 2 en temps en k-omega
 
281
do iel = 1, ncel
 
282
  rom = propce(iel,ipcroo)
 
283
  xnu = propce(iel,ipcvlo)/rom
 
284
  xk = rtpa(iel,ik)
 
285
  xw  = rtpa(iel,iomg)
 
286
  cdkw = 2*rom/ckwsw2/xw*w1(iel)
 
287
  cdkw = max(cdkw,1.d-20)
 
288
  xarg1 = max(sqrt(xk)/cmu/xw/w2(iel), 500.d0*xnu/xw/w2(iel)**2)
 
289
  xarg1 = min(xarg1, 4.d0*rom*xk/ckwsw2/cdkw/w2(iel)**2)
 
290
  xf1(iel) = tanh(xarg1**4)
 
291
enddo
 
292
 
 
293
!===============================================================================
 
294
! 4. CALCUL DU TERME DE PRODUCTION
 
295
!      Les termes sont stockes dans     TINSTK,TINSTW
 
296
!      En sortie de l'etape on conserve W1,XF1,TINSTK,TINSTW
 
297
!===============================================================================
 
298
 
 
299
d2s3 = 2.d0/3.d0
 
300
do iel = 1, ncel
 
301
  xk   = rtpa(iel,ik)
 
302
  xeps = cmu*rtpa(iel,iomg)*xk
 
303
  tinstk(iel) = s2kw(iel) - d2s3*divukw(iel)*divukw(iel)
 
304
  tinstw(iel) = propce(iel,ipcvto)*tinstk(iel)                    &
 
305
       -d2s3*propce(iel,ipcroo)*xk*divukw(iel)
 
306
  tinstk(iel) = min(tinstw(iel),ckwc1*propce(iel,ipcroo)*xeps)
 
307
enddo
 
308
 
 
309
!===============================================================================
 
310
! 4. CALCUL DU TERME DE GRAVITE
 
311
!      Les termes sont stockes dans     TINSTK,TINSTW,W2
 
312
!      En sortie de l'etape on conserve W1,W2,XF1,TINSTK,TINSTW
 
313
!===============================================================================
 
314
 
 
315
if(igrake.eq.1) then
 
316
 
 
317
  ! Allocate a temporary array for the gradient calculation
 
318
  allocate(grad(ncelet,3))
 
319
 
 
320
! --- Terme de gravite G = BETA*G*GRAD(SCA)/PRDTUR/RHO
 
321
!     Ici on calcule   G =-G*GRAD(RHO)/PRDTUR/RHO
 
322
 
 
323
  iccocg = 1
 
324
  inc = 1
 
325
 
 
326
!     Le choix ci dessous a l'avantage d'etre simple
 
327
 
 
328
  nswrgp = nswrgr(ik)
 
329
  epsrgp = epsrgr(ik)
 
330
  imligp = imligr(ik)
 
331
  iwarnp = iwarni(ik)
 
332
  climgp = climgr(ik)
 
333
  extrap = extrag(ik)
 
334
 
 
335
!     Conditions aux limites sur ROM : Dirichlet ROMB
 
336
!       On utilise VISCB pour stocker le COEFB relatif a ROM
 
337
!       On impose en Dirichlet (COEFA) la valeur ROMB
 
338
 
 
339
  do ifac = 1, nfabor
 
340
    viscb(ifac) = 0.d0
 
341
  enddo
 
342
 
 
343
  iivar = 0
 
344
 
 
345
  call grdcel                                                     &
 
346
  !==========
 
347
 ( iivar  , imrgra , inc    , iccocg , nswrgp , imligp ,          &
 
348
   iwarnp , nfecra , epsrgp , climgp , extrap ,                   &
 
349
   propce(1,ipcroo), propfb(1,ipbroo), viscb  ,                   &
 
350
   grad   )
 
351
 
 
352
 
 
353
!      Production et terme de gravite
 
354
!        TINSTK=MIN(P,C1*EPS)+G et TINSTW=P+(1-CE3)*G
 
355
!        On conserve G dans W2 pour la phase de couplage des termes sources
 
356
 
 
357
  if(iscalt.gt.0.and.nscal.ge.iscalt) then
 
358
    prdtur = sigmas(iscalt)
 
359
  else
 
360
    prdtur = 1.d0
 
361
  endif
 
362
 
 
363
  do iel = 1, ncel
 
364
    w2(iel) = -(grad(iel,1)*gx + grad(iel,2)*gy + grad(iel,3)*gz) / &
 
365
               (propce(iel,ipcroo)*prdtur)
 
366
    tinstw(iel)=tinstw(iel)+propce(iel,ipcvto)*max(w2(iel),zero)
 
367
    tinstk(iel)=tinstk(iel)+propce(iel,ipcvto)*w2(iel)
 
368
  enddo
 
369
 
 
370
  ! Free memory
 
371
  deallocate(grad)
 
372
 
 
373
endif
 
374
 
 
375
 
 
376
!===============================================================================
 
377
! 5. PRISE EN COMPTE DES TERMES SOURCES UTILISATEURS
 
378
!      Les termes sont stockes dans     SMBRK,SMBRW,DAM,W3
 
379
!      En sortie de l'etape on conserve W1-3,XF1,TINSTK,TINSTW,SMBRK,SMBRW,DAM
 
380
!===============================================================================
 
381
 
 
382
do iel = 1, ncel
 
383
  smbrk(iel) = 0.d0
 
384
  smbrw(iel) = 0.d0
 
385
  dam  (iel) = 0.d0
 
386
  w3   (iel) = 0.d0
 
387
enddo
 
388
 
 
389
call ustskw                                                       &
 
390
!==========
 
391
 ( nvar   , nscal  , ncepdp , ncesmp ,                            &
 
392
   icepdc , icetsm , itypsm ,                                     &
 
393
   dt     , rtpa   , propce , propfa , propfb ,                   &
 
394
   coefa  , coefb  , ckupdc , smacel , s2kw   , divukw ,          &
 
395
   w1     , w2     , xf1    ,                                     &
 
396
   smbrk  , smbrw  , dam    , w3     )
 
397
!        ------   ------   ------   ------
 
398
 
 
399
!===============================================================================
 
400
! 6. TERME D'ACCUMULATION DE MASSE -(dRO/dt)*Volume
 
401
 
 
402
!      Le terme est stocke dans         W4
 
403
!      En sortie de l'etape on conserve W1-4,XF1,TINSTK,TINSTW,SMBRK,SMBRW,DAM
 
404
!===============================================================================
 
405
 
 
406
init = 1
 
407
call divmas(ncelet,ncel,nfac,nfabor,init,nfecra,                  &
 
408
               ifacel,ifabor,propfa(1,iflmas),propfb(1,iflmab),w4)
 
409
 
 
410
!===============================================================================
 
411
! 7. ON FINALISE LE CALCUL DES TERMES SOURCES
 
412
 
 
413
!      Les termes sont stockes dans     SMBRK, SMBRW
 
414
!      En sortie de l'etape on conserve SMBRK,SMBRW,DAM,W1-4
 
415
!===============================================================================
 
416
 
 
417
do iel = 1, ncel
 
418
 
 
419
  visct  = propce(iel,ipcvto)
 
420
  rom    = propce(iel,ipcroo)
 
421
  xk     = rtpa(iel,ik)
 
422
  xw     = rtpa(iel,iomg)
 
423
  xxf1   = xf1(iel)
 
424
  xgamma = xxf1*ckwgm1 + (1.d0-xxf1)*ckwgm2
 
425
  xbeta  = xxf1*ckwbt1 + (1.d0-xxf1)*ckwbt2
 
426
 
 
427
  smbrk(iel) = smbrk(iel) + volume(iel)*(                         &
 
428
       tinstk(iel)                                                &
 
429
       -cmu*rom*xw*xk )
 
430
 
 
431
  smbrw(iel) = smbrw(iel) + volume(iel)*(                         &
 
432
       rom*xgamma/visct*tinstw(iel)                               &
 
433
       -xbeta*rom*xw**2                                           &
 
434
       +2.d0*rom/xw*(1.d0-xxf1)/ckwsw2*w1(iel) )
 
435
 
 
436
enddo
 
437
 
 
438
!===============================================================================
 
439
! 8. PRISE EN COMPTE DES TERMES D'ACCUMULATION DE MASSE ET
 
440
!         DE LA DEUXIEME PARTIE DES TS UTILISATEURS (PARTIE EXPLICITE)
 
441
!         STOCKAGE POUR EXTRAPOLATION EN TEMPS
 
442
!      On utilise                       SMBRK,SMBRW
 
443
!      En sortie de l'etape on conserve SMBRK,SMBRW,DAM,W2-4
 
444
 
 
445
!    Remarque : l'extrapolation telle qu'elle est ecrite n'a pas grand
 
446
!               sens si IKECOU=1
 
447
!===============================================================================
 
448
 
 
449
!     Si on extrapole les T.S.
 
450
if(isto2t.gt.0) then
 
451
 
 
452
  do iel = 1, ncel
 
453
 
 
454
!       Sauvegarde pour echange
 
455
    tuexpk = propce(iel,iptsta)
 
456
!       Pour la suite et le pas de temps suivant
 
457
    propce(iel,iptsta) = smbrk(iel)
 
458
!       Termes dependant de la variable resolue et theta PROPCE
 
459
    smbrk(iel) = iconv(ik)*w4(iel)*rtpa(iel,ik) - thets*tuexpk
 
460
!       On suppose -DAM > 0 : on implicite
 
461
!         le terme utilisateur dependant de la variable resolue
 
462
    smbrk(iel) = dam(iel)*rtpa(iel,ik) + smbrk(iel)
 
463
 
 
464
!       Sauvegarde pour echange
 
465
    tuexpw = propce(iel,iptsta+1)
 
466
!       Pour la suite et le pas de temps suivant
 
467
    propce(iel,iptsta+1) = smbrw(iel)
 
468
!       Termes dependant de la variable resolue et theta PROPCE
 
469
    smbrw(iel) = iconv(iomg)*w4(iel)*rtpa(iel,iomg) - thets*tuexpw
 
470
!       On suppose -W3 > 0 : on implicite
 
471
!         le terme utilisateur dependant de la variable resolue
 
472
    smbrw(iel) =  w3(iel)*rtpa(iel,iomg) + smbrw(iel)
 
473
 
 
474
  enddo
 
475
 
 
476
!     Si on n'extrapole pas les T.S.
 
477
else
 
478
  do iel = 1, ncel
 
479
    smbrk(iel) = smbrk(iel) + dam(iel)*rtpa(iel,ik)            &
 
480
         +iconv(ik)*w4(iel)*rtpa(iel,ik)
 
481
    smbrw(iel) = smbrw(iel) + w3 (iel)*rtpa(iel,iomg)           &
 
482
         +iconv(iomg)*w4(iel)*rtpa(iel,iomg)
 
483
  enddo
 
484
endif
 
485
 
 
486
!===============================================================================
 
487
! 8.1 PRISE EN COMPTE DES TERMES SOURCES LAGRANGIEN : PARTIE EXPLICITE
 
488
!     COUPLAGE RETOUR
 
489
!===============================================================================
 
490
 
 
491
!     Ordre 2 non pris en compte
 
492
if (iilagr.eq.2 .and. ltsdyn.eq.1) then
 
493
 
 
494
  do iel = 1,ncel
 
495
 
 
496
! Termes sources explicte et implicte sur k
 
497
 
 
498
    smbrk(iel)  = smbrk(iel) + tslagr(iel,itske)
 
499
 
 
500
! Termes sources explicte sur omega : on reprend la constante CE4 directement
 
501
!    du k-eps sans justification ... a creuser si necessaire           !
 
502
 
 
503
    smbrw(iel)  = smbrw(iel)                                      &
 
504
                + ce4 *tslagr(iel,itske) * propce(iel,ipcroo)     &
 
505
                /propce(iel,ipcvto)
 
506
  enddo
 
507
 
 
508
endif
 
509
 
 
510
 
 
511
!===============================================================================
 
512
! 9. PRISE EN COMPTE DES TERMES DE CONV/DIFF DANS LE SECOND MEMBRE
 
513
 
 
514
!      Tableaux de travail              W7, W8, W1, TINSTK, TINSTW
 
515
!      Les termes sont stockes dans     W5 et W6, puis ajoutes a SMBRK, SMBRW
 
516
!      En sortie de l'etape on conserve W2-6,SMBRK,SMBRW,DAM
 
517
!===============================================================================
 
518
 
 
519
!     Ceci ne sert a rien si IKECOU n'est pas egal a 1
 
520
 
 
521
if (ikecou.eq.1) then
 
522
 
 
523
  do iel = 1, ncel
 
524
    w5 (iel) = 0.d0
 
525
    w6 (iel) = 0.d0
 
526
  enddo
 
527
 
 
528
! ---> Traitement de k
 
529
 
 
530
  ivar   = ik
 
531
 
 
532
  ipp    = ipprtp(ivar)
 
533
 
 
534
  iclvar = iclrtp(ivar,icoef )
 
535
  iclvaf = iclrtp(ivar,icoeff)
 
536
  chaine = nomvar(ipp)
 
537
 
 
538
  if( idiff(ivar).ge. 1 ) then
 
539
 
 
540
    do iel = 1, ncel
 
541
      xxf1 = xf1(iel)
 
542
      sigma = xxf1*ckwsk1 + (1.d0-xxf1)*ckwsk2
 
543
      w7(iel) = propce(iel,ipcvis)                                &
 
544
           + idifft(ivar)*propce(iel,ipcvst)/sigma
 
545
    enddo
 
546
    call viscfa                                                   &
 
547
    !==========
 
548
 ( imvisf ,                                                       &
 
549
   w7     ,                                                       &
 
550
   viscf  , viscb  )
 
551
 
 
552
  else
 
553
 
 
554
    do ifac = 1, nfac
 
555
      viscf(ifac) = 0.d0
 
556
    enddo
 
557
    do ifac = 1, nfabor
 
558
      viscb(ifac) = 0.d0
 
559
    enddo
 
560
 
 
561
  endif
 
562
 
 
563
  iccocg = 1
 
564
  inc    = 1
 
565
  iconvp = iconv (ivar)
 
566
  idiffp = idiff (ivar)
 
567
  nswrgp = nswrgr(ivar)
 
568
  imligp = imligr(ivar)
 
569
  ircflp = ircflu(ivar)
 
570
  ischcp = ischcv(ivar)
 
571
  isstpp = isstpc(ivar)
 
572
  iwarnp = iwarni(ivar)
 
573
  blencp = blencv(ivar)
 
574
  epsrgp = epsrgr(ivar)
 
575
  climgp = climgr(ivar)
 
576
  extrap = extrag(ivar)
 
577
  relaxp = relaxv(ivar)
 
578
  thetap = thetav(ivar)
 
579
 
 
580
  call bilsc2                                                     &
 
581
  !==========
 
582
 ( nvar   , nscal  ,                                              &
 
583
   idtvar , ivar   , iconvp , idiffp , nswrgp , imligp , ircflp , &
 
584
   ischcp , isstpp , inc    , imrgra , iccocg ,                   &
 
585
   ipp    , iwarnp ,                                              &
 
586
   blencp , epsrgp , climgp , extrap , relaxp , thetap ,          &
 
587
   rtpa(1,ivar)    , rtpa(1,ivar)    ,                            &
 
588
   coefa(1,iclvar) , coefb(1,iclvar) ,                            &
 
589
   coefa(1,iclvaf) , coefb(1,iclvaf) ,                            &
 
590
   propfa(1,iflmas), propfb(1,iflmab), viscf  , viscb  ,          &
 
591
   w5     )
 
592
 
 
593
 
 
594
  if (iwarni(ivar).ge.2) then
 
595
    isqrt = 1
 
596
    call prodsc(ncelet,ncel,isqrt,smbrk,smbrk,rnorm)
 
597
    write(nfecra,1100) chaine(1:8) ,rnorm
 
598
  endif
 
599
 
 
600
 
 
601
! ---> Traitement de omega
 
602
 
 
603
  ivar   = iomg
 
604
 
 
605
  ipp    = ipprtp(ivar)
 
606
 
 
607
  iclvar = iclrtp(ivar,icoef )
 
608
  iclvaf = iclrtp(ivar,icoeff)
 
609
  chaine = nomvar(ipp)
 
610
 
 
611
  if( idiff(ivar).ge. 1 ) then
 
612
    do iel = 1, ncel
 
613
      xxf1 = xf1(iel)
 
614
      sigma = xxf1*ckwsw1 + (1.d0-xxf1)*ckwsw2
 
615
      w7(iel) = propce(iel,ipcvis)                                &
 
616
           + idifft(ivar)*propce(iel,ipcvst)/sigma
 
617
    enddo
 
618
    call viscfa                                                   &
 
619
    !==========
 
620
 ( imvisf ,                                                       &
 
621
   w7     ,                                                       &
 
622
   viscf  , viscb  )
 
623
 
 
624
  else
 
625
 
 
626
    do ifac = 1, nfac
 
627
      viscf(ifac) = 0.d0
 
628
    enddo
 
629
    do ifac = 1, nfabor
 
630
      viscb(ifac) = 0.d0
 
631
    enddo
 
632
 
 
633
  endif
 
634
 
 
635
  iccocg = 1
 
636
  inc    = 1
 
637
  iconvp = iconv (ivar)
 
638
  idiffp = idiff (ivar)
 
639
  nswrgp = nswrgr(ivar)
 
640
  imligp = imligr(ivar)
 
641
  ircflp = ircflu(ivar)
 
642
  ischcp = ischcv(ivar)
 
643
  isstpp = isstpc(ivar)
 
644
  iwarnp = iwarni(ivar)
 
645
  blencp = blencv(ivar)
 
646
  epsrgp = epsrgr(ivar)
 
647
  climgp = climgr(ivar)
 
648
  extrap = extrag(ivar)
 
649
  relaxp = relaxv(ivar)
 
650
  thetap = thetav(ivar)
 
651
 
 
652
  call bilsc2                                                     &
 
653
  !==========
 
654
 ( nvar   , nscal  ,                                              &
 
655
   idtvar , ivar   , iconvp , idiffp , nswrgp , imligp , ircflp , &
 
656
   ischcp , isstpp , inc    , imrgra , iccocg ,                   &
 
657
   ipp    , iwarnp ,                                              &
 
658
   blencp , epsrgp , climgp , extrap , relaxp , thetap ,          &
 
659
   rtpa(1,ivar)    , rtpa(1,ivar)    ,                            &
 
660
   coefa(1,iclvar) , coefb(1,iclvar) ,                            &
 
661
   coefa(1,iclvaf) , coefb(1,iclvaf) ,                            &
 
662
   propfa(1,iflmas), propfb(1,iflmab), viscf  , viscb  ,          &
 
663
   w6  )
 
664
!        --
 
665
 
 
666
  if (iwarni(ivar).ge.2) then
 
667
    isqrt = 1
 
668
    call prodsc(ncelet,ncel,isqrt,smbrw,smbrw,rnorm)
 
669
    write(nfecra,1100) chaine(1:8) ,rnorm
 
670
  endif
 
671
 
 
672
  do iel = 1,ncel
 
673
    smbrk(iel) = smbrk(iel) + w5(iel)
 
674
    smbrw(iel) = smbrw(iel) + w6(iel)
 
675
  enddo
 
676
 
 
677
endif
 
678
 
 
679
!===============================================================================
 
680
! 10. AJOUT DES TERMES SOURCES DE MASSE EXPLICITES
 
681
 
 
682
!       Les parties implicites eventuelles sont conservees dans W7 et W8
 
683
!         et utilisees dans la phase d'implicitation cv/diff
 
684
 
 
685
!       Les termes sont stockes dans     SMBRK, SMBRW, W7, W8
 
686
!       En sortie de l'etape on conserve W2-8,SMBRK,SMBRW,DAM
 
687
!===============================================================================
 
688
 
 
689
if (ncesmp.gt.0) then
 
690
 
 
691
  do iel = 1, ncel
 
692
    w7(iel) = 0.d0
 
693
    w8(iel) = 0.d0
 
694
  enddo
 
695
 
 
696
!       Entier egal a 1 (pour navsto : nb de sur-iter)
 
697
  iiun = 1
 
698
 
 
699
!       On incremente SMBRS par -Gamma RTPA et ROVSDT par Gamma (*theta)
 
700
  ivar = ik
 
701
  call catsma &
 
702
  !==========
 
703
 ( ncelet , ncel   , ncesmp , iiun   ,                            &
 
704
                                 isto2t , thetav(ivar) ,   &
 
705
   icetsm , itypsm(1,ivar) ,                                      &
 
706
   volume , rtpa(1,ivar) , smacel(1,ivar) , smacel(1,ipr) ,       &
 
707
   smbrk  , w7     , tinstk )
 
708
  ivar = iomg
 
709
  call catsma &
 
710
  !==========
 
711
 ( ncelet , ncel   , ncesmp , iiun   ,                            &
 
712
                                 isto2t , thetav(ivar) ,   &
 
713
   icetsm , itypsm(1,ivar) ,                                      &
 
714
   volume , rtpa(1,ivar) , smacel(1,ivar) , smacel(1,ipr) ,       &
 
715
   smbrw  , w8     , tinstw )
 
716
 
 
717
!       Si on extrapole les TS on met Gamma Pinj dans PROPCE
 
718
  if(isto2t.gt.0) then
 
719
    do iel = 1, ncel
 
720
      propce(iel,iptsta  ) = propce(iel,iptsta  ) + tinstk(iel)
 
721
      propce(iel,iptsta+1) = propce(iel,iptsta+1) + tinstw(iel)
 
722
    enddo
 
723
!       Sinon on le met directement dans SMBR
 
724
  else
 
725
    do iel = 1, ncel
 
726
      smbrk(iel) = smbrk(iel) + tinstk(iel)
 
727
      smbrw(iel) = smbrw(iel) + tinstw(iel)
 
728
    enddo
 
729
  endif
 
730
 
 
731
endif
 
732
 
 
733
!     Finalisation des termes sources
 
734
if(isto2t.gt.0) then
 
735
  thetp1 = 1.d0 + thets
 
736
  do iel = 1, ncel
 
737
    smbrk(iel) = smbrk(iel) + thetp1 * propce(iel,iptsta)
 
738
    smbrw(iel) = smbrw(iel) + thetp1 * propce(iel,iptsta+1)
 
739
  enddo
 
740
endif
 
741
 
 
742
!===============================================================================
 
743
! 11. INCREMENTS DES TERMES SOURCES DANS LE SECOND MEMBRE
 
744
 
 
745
!       Les termes sont stockes dans     SMBRK, SMBRW
 
746
!       En sortie de l'etape on conserve W3-8,SMBRK,SMBRW
 
747
!===============================================================================
 
748
 
 
749
!     Ordre 2 non pris en compte
 
750
if(ikecou.eq.1) then
 
751
 
 
752
  do iel = 1, ncel
 
753
 
 
754
    rom = propce(iel,ipcrom)
 
755
 
 
756
!   RESOLUTION COUPLEE
 
757
 
 
758
    romvsd     = 1.d0/(rom*volume(iel))
 
759
    smbrk(iel) = smbrk(iel)*romvsd
 
760
    smbrw(iel) = smbrw(iel)*romvsd
 
761
    divp23     = d2s3*max(divukw(iel),zero)
 
762
    produc     = s2kw(iel)-d2s3*divukw(iel)**2+w2(iel)
 
763
    xk         = rtpa(iel,ik)
 
764
    xw         = rtpa(iel,iomg)
 
765
    xxf1       = xf1(iel)
 
766
    xgamma     = xxf1*ckwgm1 + (1.d0-xxf1)*ckwgm2
 
767
    xbeta      = xxf1*ckwbt1 + (1.d0-xxf1)*ckwbt2
 
768
 
 
769
    a11 = 1.d0/dt(iel)                                            &
 
770
         - 1.d0/xw*min(produc,zero)+divp23+cmu*xw
 
771
    a12 = cmu*xk
 
772
    a21 = 0.d0
 
773
    a22 = 1.d0/dt(iel)+xgamma*divp23+2.d0*xbeta*xw
 
774
 
 
775
    unsdet = 1.d0/(a11*a22 -a12*a21)
 
776
 
 
777
    deltk = ( a22*smbrk(iel) -a12*smbrw(iel) )*unsdet
 
778
    deltw = (-a21*smbrk(iel) +a11*smbrw(iel) )*unsdet
 
779
 
 
780
!     NOUVEAU TERME SOURCE POUR CODITS
 
781
 
 
782
    romvsd = rom*volume(iel)/dt(iel)
 
783
 
 
784
    smbrk(iel) = romvsd*deltk
 
785
    smbrw(iel) = romvsd*deltw
 
786
 
 
787
  enddo
 
788
 
 
789
endif
 
790
 
 
791
!===============================================================================
 
792
! 12. TERMES INSTATIONNAIRES
 
793
 
 
794
!     Les termes sont stockes dans     TINSTK, TINSTW
 
795
!     En sortie de l'etape on conserve SMBRK, SMBRW,  TINSTK, TINSTW
 
796
!===============================================================================
 
797
 
 
798
! --- PARTIE EXPLICITE
 
799
 
 
800
!     on enleve la convection/diffusion au temps n a SMBRK et SMBRW
 
801
!     s'ils ont ete calcules
 
802
if (ikecou.eq.1) then
 
803
  do iel = 1, ncel
 
804
    smbrk(iel) = smbrk(iel) - w5(iel)
 
805
    smbrw(iel) = smbrw(iel) - w6(iel)
 
806
  enddo
 
807
endif
 
808
 
 
809
! --- RHO/DT et DIV
 
810
!     Extrapolation ou non, meme forme par coherence avec bilsc2
 
811
 
 
812
do iel = 1, ncel
 
813
  rom = propce(iel,ipcrom)
 
814
  romvsd = rom*volume(iel)/dt(iel)
 
815
  tinstk(iel) = istat(ik)*romvsd   -iconv(ik)*w4(iel)*thetav(ik)
 
816
  tinstw(iel) = istat(iomg)*romvsd -iconv(iomg)*w4(iel)*thetav(iomg)
 
817
enddo
 
818
 
 
819
! --- Source de masse (le theta est deja inclus par catsma)
 
820
if (ncesmp.gt.0) then
 
821
  do iel = 1, ncel
 
822
    tinstk(iel) = tinstk(iel) + w7(iel)
 
823
    tinstw(iel) = tinstw(iel) + w8(iel)
 
824
  enddo
 
825
endif
 
826
 
 
827
! --- Termes sources utilisateurs
 
828
if(isto2t.gt.0) then
 
829
  thetak = thetav(ik)
 
830
  thetaw = thetav(iomg)
 
831
  do iel = 1, ncel
 
832
    tinstk(iel) = tinstk(iel) -dam(iel)*thetak
 
833
    tinstw(iel) = tinstw(iel) -w3 (iel)*thetaw
 
834
  enddo
 
835
else
 
836
  do iel = 1, ncel
 
837
    tinstk(iel) = tinstk(iel) + max(-dam(iel),zero)
 
838
    tinstw(iel) = tinstw(iel) + max(-w3 (iel),zero)
 
839
  enddo
 
840
endif
 
841
 
 
842
! --- PRISE EN COMPTE DES TERMES LAGRANGIEN : COUPLAGE RETOUR
 
843
 
 
844
!     Ordre 2 non pris en compte
 
845
if (iilagr.eq.2 .and. ltsdyn.eq.1) then
 
846
 
 
847
  do iel = 1,ncel
 
848
 
 
849
! Termes sources implicite sur k
 
850
 
 
851
    tinstk(iel) = tinstk(iel) + max(-tslagr(iel,itsli),zero)
 
852
 
 
853
! Termes sources implicte sur omega
 
854
 
 
855
    tinstw(iel) = tinstw(iel)                                     &
 
856
          + max( (-ce4*tslagr(iel,itske)/rtpa(iel,ik)) , zero)
 
857
 
 
858
  enddo
 
859
 
 
860
endif
 
861
 
 
862
! Si IKECOU=0, on implicite plus fortement k et omega
 
863
 
 
864
if(ikecou.eq.0)then
 
865
  do iel=1,ncel
 
866
    xw    = rtpa(iel,iomg)
 
867
    xxf1  = xf1(iel)
 
868
    xbeta = xxf1*ckwbt1 + (1.d0-xxf1)*ckwbt2
 
869
    rom = propce(iel,ipcrom)
 
870
      tinstk(iel) = tinstk(iel) +                                 &
 
871
         volume(iel)*cmu*rom*xw
 
872
      tinstw(iel) = tinstw(iel) +                                 &
 
873
         volume(iel)*xbeta*rom*xw
 
874
  enddo
 
875
endif
 
876
 
 
877
 
 
878
!===============================================================================
 
879
! 13. RESOLUTION
 
880
 
 
881
!       On utilise                      SMBRK, SMBRW,  TINSTK, TINSTW
 
882
!===============================================================================
 
883
 
 
884
! ---> Traitement de k
 
885
 
 
886
ivar = ik
 
887
iclvar = iclrtp(ivar,icoef )
 
888
iclvaf = iclrtp(ivar,icoeff)
 
889
 
 
890
ipp    = ipprtp(ivar)
 
891
 
 
892
!     "VITESSE" DE DIFFUSION FACETTE
 
893
 
 
894
if( idiff(ivar).ge. 1 ) then
 
895
 
 
896
  do iel = 1, ncel
 
897
    xxf1 = xf1(iel)
 
898
    sigma = xxf1*ckwsk1 + (1.d0-xxf1)*ckwsk2
 
899
    w1(iel) = propce(iel,ipcvis)                                  &
 
900
                        + idifft(ivar)*propce(iel,ipcvst)/sigma
 
901
  enddo
 
902
  call viscfa                                                     &
 
903
  !==========
 
904
 ( imvisf ,                                                       &
 
905
   w1     ,                                                       &
 
906
   viscf  , viscb  )
 
907
 
 
908
else
 
909
 
 
910
  do ifac = 1, nfac
 
911
    viscf(ifac) = 0.d0
 
912
  enddo
 
913
  do ifac = 1, nfabor
 
914
    viscb(ifac) = 0.d0
 
915
  enddo
 
916
 
 
917
endif
 
918
 
 
919
!     RESOLUTION POUR K
 
920
 
 
921
iconvp = iconv (ivar)
 
922
idiffp = idiff (ivar)
 
923
ireslp = iresol(ivar)
 
924
ndircp = ndircl(ivar)
 
925
nitmap = nitmax(ivar)
 
926
nswrsp = nswrsm(ivar)
 
927
nswrgp = nswrgr(ivar)
 
928
imligp = imligr(ivar)
 
929
ircflp = ircflu(ivar)
 
930
ischcp = ischcv(ivar)
 
931
isstpp = isstpc(ivar)
 
932
iescap = 0
 
933
imgrp  = imgr  (ivar)
 
934
ncymxp = ncymax(ivar)
 
935
nitmfp = nitmgf(ivar)
 
936
!MO      IPP    =
 
937
iwarnp = iwarni(ivar)
 
938
blencp = blencv(ivar)
 
939
epsilp = epsilo(ivar)
 
940
epsrsp = epsrsm(ivar)
 
941
epsrgp = epsrgr(ivar)
 
942
climgp = climgr(ivar)
 
943
extrap = extrag(ivar)
 
944
relaxp = relaxv(ivar)
 
945
thetap = thetav(ivar)
 
946
 
 
947
call codits                                                       &
 
948
!==========
 
949
 ( nvar   , nscal  ,                                              &
 
950
   idtvar , ivar   , iconvp , idiffp , ireslp , ndircp , nitmap , &
 
951
   imrgra , nswrsp , nswrgp , imligp , ircflp ,                   &
 
952
   ischcp , isstpp , iescap ,                                     &
 
953
   imgrp  , ncymxp , nitmfp , ipp    , iwarnp ,                   &
 
954
   blencp , epsilp , epsrsp , epsrgp , climgp , extrap ,          &
 
955
   relaxp , thetap ,                                              &
 
956
   rtpa(1,ivar)    , rtpa(1,ivar)    ,                            &
 
957
                     coefa(1,iclvar) , coefb(1,iclvar) ,          &
 
958
                     coefa(1,iclvaf) , coefb(1,iclvaf) ,          &
 
959
                     propfa(1,iflmas), propfb(1,iflmab),          &
 
960
   viscf  , viscb  , viscf  , viscb  ,                            &
 
961
   tinstk , smbrk  , rtp(1,ivar)     ,                            &
 
962
   rvoid  )
 
963
 
 
964
 
 
965
! ---> Traitement de omega
 
966
 
 
967
ivar = iomg
 
968
iclvar = iclrtp(ivar,icoef )
 
969
iclvaf = iclrtp(ivar,icoeff)
 
970
 
 
971
ipp    = ipprtp(ivar)
 
972
 
 
973
 
 
974
!     "VITESSE" DE DIFFUSION FACETTE
 
975
 
 
976
if( idiff(ivar).ge. 1 ) then
 
977
  do iel = 1, ncel
 
978
    xxf1 = xf1(iel)
 
979
    sigma = xxf1*ckwsw1 + (1.d0-xxf1)*ckwsw2
 
980
    w1(iel) = propce(iel,ipcvis)                                  &
 
981
                        + idifft(ivar)*propce(iel,ipcvst)/sigma
 
982
  enddo
 
983
  call viscfa                                                     &
 
984
  !==========
 
985
 ( imvisf ,                                                       &
 
986
   w1     ,                                                       &
 
987
   viscf  , viscb  )
 
988
 
 
989
else
 
990
 
 
991
  do ifac = 1, nfac
 
992
    viscf(ifac) = 0.d0
 
993
  enddo
 
994
  do ifac = 1, nfabor
 
995
    viscb(ifac) = 0.d0
 
996
  enddo
 
997
 
 
998
endif
 
999
 
 
1000
! RESOLUTION POUR OMEGA
 
1001
 
 
1002
iconvp = iconv (ivar)
 
1003
idiffp = idiff (ivar)
 
1004
ireslp = iresol(ivar)
 
1005
ndircp = ndircl(ivar)
 
1006
nitmap = nitmax(ivar)
 
1007
nswrsp = nswrsm(ivar)
 
1008
nswrgp = nswrgr(ivar)
 
1009
imligp = imligr(ivar)
 
1010
ircflp = ircflu(ivar)
 
1011
ischcp = ischcv(ivar)
 
1012
isstpp = isstpc(ivar)
 
1013
iescap = 0
 
1014
imgrp  = imgr  (ivar)
 
1015
ncymxp = ncymax(ivar)
 
1016
nitmfp = nitmgf(ivar)
 
1017
!MO      IPP    =
 
1018
iwarnp = iwarni(ivar)
 
1019
blencp = blencv(ivar)
 
1020
epsilp = epsilo(ivar)
 
1021
epsrsp = epsrsm(ivar)
 
1022
epsrgp = epsrgr(ivar)
 
1023
climgp = climgr(ivar)
 
1024
extrap = extrag(ivar)
 
1025
relaxp = relaxv(ivar)
 
1026
thetap = thetav(ivar)
 
1027
 
 
1028
call codits                                                       &
 
1029
!==========
 
1030
 ( nvar   , nscal  ,                                              &
 
1031
   idtvar , ivar   , iconvp , idiffp , ireslp , ndircp , nitmap , &
 
1032
   imrgra , nswrsp , nswrgp , imligp , ircflp ,                   &
 
1033
   ischcp , isstpp , iescap ,                                     &
 
1034
   imgrp  , ncymxp , nitmfp , ipp    , iwarnp ,                   &
 
1035
   blencp , epsilp , epsrsp , epsrgp , climgp , extrap ,          &
 
1036
   relaxp , thetap ,                                              &
 
1037
   rtpa(1,ivar)    , rtpa(1,ivar)    ,                            &
 
1038
                     coefa(1,iclvar) , coefb(1,iclvar) ,          &
 
1039
                     coefa(1,iclvaf) , coefb(1,iclvaf) ,          &
 
1040
                     propfa(1,iflmas), propfb(1,iflmab),          &
 
1041
   viscf  , viscb  , viscf  , viscb  ,                            &
 
1042
   tinstw , smbrw  , rtp(1,ivar)     ,                            &
 
1043
   rvoid  )
 
1044
 
 
1045
 
 
1046
!===============================================================================
 
1047
! 14. CLIPPING
 
1048
!===============================================================================
 
1049
 
 
1050
!     Calcul des Min/Max avant clipping, pour affichage
 
1051
do ii = 1, 2
 
1052
  if(ii.eq.1) then
 
1053
    ivar = ik
 
1054
  elseif(ii.eq.2) then
 
1055
    ivar = iomg
 
1056
  endif
 
1057
  ipp  = ipprtp(ivar)
 
1058
 
 
1059
  vrmin =  grand
 
1060
  vrmax = -grand
 
1061
  do iel = 1, ncel
 
1062
    var = rtp(iel,ivar)
 
1063
    vrmin = min(vrmin,var)
 
1064
    vrmax = max(vrmax,var)
 
1065
  enddo
 
1066
  if (irangp.ge.0) then
 
1067
    call parmax (vrmax)
 
1068
    !==========
 
1069
    call parmin (vrmin)
 
1070
    !==========
 
1071
  endif
 
1072
  varmna(ipp) = vrmin
 
1073
  varmxa(ipp) = vrmax
 
1074
 
 
1075
enddo
 
1076
 
 
1077
!     On clippe simplement k et omega par valeur absolue
 
1078
iclipk = 0
 
1079
iclipw = 0
 
1080
do iel = 1, ncel
 
1081
  xk = rtp(iel,ik )
 
1082
  xw = rtp(iel,iomg)
 
1083
  if (abs(xk).le.epz2) then
 
1084
    iclipk = iclipk + 1
 
1085
    rtp(iel,ik) = max(rtp(iel,ik),epz2)
 
1086
  elseif(xk.le.0.d0) then
 
1087
    iclipk = iclipk + 1
 
1088
    rtp(iel,ik) = -xk
 
1089
  endif
 
1090
  if (abs(xw).le.epz2) then
 
1091
    iclipw = iclipw + 1
 
1092
    rtp(iel,iomg) = max(rtp(iel,iomg),epz2)
 
1093
  elseif(xw.le.0.d0) then
 
1094
    iclipw = iclipw + 1
 
1095
    rtp(iel,iomg) = -xw
 
1096
  endif
 
1097
enddo
 
1098
 
 
1099
if (irangp.ge.0) then
 
1100
  call parcpt (iclipk)
 
1101
  !==========
 
1102
  call parcpt (iclipw)
 
1103
  !==========
 
1104
endif
 
1105
 
 
1106
! ---  Stockage nb de clippings pour listing
 
1107
 
 
1108
iclpmn(ipprtp(ik )) = iclipk
 
1109
iclpmn(ipprtp(iomg)) = iclipw
 
1110
 
 
1111
 
 
1112
! Free memory
 
1113
deallocate(viscf, viscb)
 
1114
deallocate(dam)
 
1115
deallocate(smbrk, smbrw, rovsdt)
 
1116
deallocate(tinstk, tinstw, xf1)
 
1117
deallocate(w1, w2, w3)
 
1118
deallocate(w4, w5, w6)
 
1119
deallocate(w7, w8)
 
1120
 
 
1121
!--------
 
1122
! FORMATS
 
1123
!--------
 
1124
 
 
1125
#if defined(_CS_LANG_FR)
 
1126
 
 
1127
 1000 format(/,                                                   &
 
1128
'   ** RESOLUTION DU K-OMEGA                     ',/,&
 
1129
'      ---------------------                     ',/)
 
1130
 1100 format(1X,A8,' : BILAN EXPLICITE = ',E14.5)
 
1131
 
 
1132
#else
 
1133
 
 
1134
 1000 format(/,                                                   &
 
1135
'   ** SOLVING K-OMEGA'                           ,/,&
 
1136
'      ---------------'                           ,/)
 
1137
 1100 format(1X,A8,' : EXPLICIT BALANCE = ',E14.5)
 
1138
 
 
1139
#endif
 
1140
 
 
1141
!----
 
1142
! FIN
 
1143
!----
 
1144
 
 
1145
return
 
1146
 
 
1147
end subroutine