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

« back to all changes in this revision

Viewing changes to src/cogz/ebuini.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:
1
1
!-------------------------------------------------------------------------------
2
2
 
3
 
!     This file is part of the Code_Saturne Kernel, element of the
4
 
!     Code_Saturne CFD tool.
5
 
 
6
 
!     Copyright (C) 1998-2009 EDF S.A., France
7
 
 
8
 
!     contact: saturne-support@edf.fr
9
 
 
10
 
!     The Code_Saturne Kernel is free software; you can redistribute it
11
 
!     and/or modify it under the terms of the GNU General Public License
12
 
!     as published by the Free Software Foundation; either version 2 of
13
 
!     the License, or (at your option) any later version.
14
 
 
15
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
16
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18
 
!     GNU General Public License for more details.
19
 
 
20
 
!     You should have received a copy of the GNU General Public License
21
 
!     along with the Code_Saturne Kernel; if not, write to the
22
 
!     Free Software Foundation, Inc.,
23
 
!     51 Franklin St, Fifth Floor,
24
 
!     Boston, MA  02110-1301  USA
 
3
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
4
!
 
5
! Copyright (C) 1998-2011 EDF S.A.
 
6
!
 
7
! This program is free software; you can redistribute it and/or modify it under
 
8
! the terms of the GNU General Public License as published by the Free Software
 
9
! Foundation; either version 2 of the License, or (at your option) any later
 
10
! version.
 
11
!
 
12
! This program is distributed in the hope that it will be useful, but WITHOUT
 
13
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 
14
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 
15
! details.
 
16
!
 
17
! You should have received a copy of the GNU General Public License along with
 
18
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
 
19
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
25
20
 
26
21
!-------------------------------------------------------------------------------
27
22
 
28
23
subroutine ebuini &
29
24
!================
30
25
 
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
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , &
37
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
38
 
   idevel , ituser , ia     ,                                     &
39
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
40
 
   dt     , rtp    , propce , propfa , propfb , coefa  , coefb  , &
41
 
   rdevel , rtuser , ra     )
 
26
 ( nvar   , nscal  ,                                              &
 
27
   dt     , rtp    , propce , propfa , propfb , coefa  , coefb  )
42
28
 
43
29
!===============================================================================
44
30
! FONCTION :
66
52
!     PROPCE (prop au centre), PROPFA (aux faces internes),
67
53
!     PROPFB (prop aux faces de bord)
68
54
!     Ainsi,
69
 
!      PROPCE(IEL,IPPROC(IROM  (IPHAS))) designe ROM   (IEL ,IPHAS)
70
 
!      PROPCE(IEL,IPPROC(IVISCL(IPHAS))) designe VISCL (IEL ,IPHAS)
71
 
!      PROPCE(IEL,IPPROC(ICP   (IPHAS))) designe CP    (IEL ,IPHAS)
 
55
!      PROPCE(IEL,IPPROC(IROM  )) designe ROM   (IEL)
 
56
!      PROPCE(IEL,IPPROC(IVISCL)) designe VISCL (IEL)
 
57
!      PROPCE(IEL,IPPROC(ICP   )) designe CP    (IEL)
72
58
!      PROPCE(IEL,IPPROC(IVISLS(ISCAL))) designe VISLS (IEL ,ISCAL)
73
59
 
74
60
!      PROPFA(IFAC,IPPROF(IFLUMA(IVAR ))) designe FLUMAS(IFAC,IVAR)
75
61
 
76
 
!      PROPFB(IFAC,IPPROB(IROM  (IPHAS))) designe ROMB  (IFAC,IPHAS)
 
62
!      PROPFB(IFAC,IPPROB(IROM  )) designe ROMB  (IFAC)
77
63
!      PROPFB(IFAC,IPPROB(IFLUMA(IVAR ))) designe FLUMAB(IFAC,IVAR)
78
64
 
79
65
! LA MODIFICATION DES PROPRIETES PHYSIQUES (ROM, VISCL, VISCLS, CP)
84
70
!__________________.____._____.________________________________________________.
85
71
! name             !type!mode ! role                                           !
86
72
!__________________!____!_____!________________________________________________!
87
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
88
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
89
 
! ndim             ! i  ! <-- ! spatial dimension                              !
90
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
91
 
! ncel             ! i  ! <-- ! number of cells                                !
92
 
! nfac             ! i  ! <-- ! number of interior faces                       !
93
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
94
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
95
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
96
 
! nnod             ! i  ! <-- ! number of vertices                             !
97
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
98
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
99
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
100
73
! nvar             ! i  ! <-- ! total number of variables                      !
101
74
! nscal            ! i  ! <-- ! total number of scalars                        !
102
 
! nphas            ! i  ! <-- ! number of phases                               !
103
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
104
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
105
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
106
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
107
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
108
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
109
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
110
 
!  (nfml, nprfml)  !    !     !                                                !
111
 
! maxelt           ! i  ! <-- ! max number of cells and faces (int/boundary)   !
112
 
! lstelt(maxelt)   ! ia ! --- ! work array                                     !
113
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
114
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
115
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
116
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
117
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
118
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
119
 
! ia(*)            ! ia ! --- ! main integer work array                        !
120
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
121
 
!  (ndim, ncelet)  !    !     !                                                !
122
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
123
 
!  (ndim, nfac)    !    !     !                                                !
124
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
125
 
!  (ndim, nfabor)  !    !     !                                                !
126
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
127
 
!  (ndim, nfac)    !    !     !                                                !
128
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
129
 
!  (ndim, nfabor)  !    !     !                                                !
130
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
131
 
!  (ndim, nnod)    !    !     !                                                !
132
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
133
75
! dt(ncelet)       ! tr ! <-- ! valeur du pas de temps                         !
134
76
! rtp              ! tr ! <-- ! variables de calcul au centre des              !
135
77
! (ncelet,*)       !    !     !    cellules                                    !
138
80
! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers   !
139
81
! coefa coefb      ! tr ! <-- ! conditions aux limites aux                     !
140
82
!  (nfabor,*)      !    !     !    faces de bord                               !
141
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
142
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
143
 
! ra(*)            ! ra ! --- ! main real work array                           !
144
83
!__________________!____!_____!________________________________________________!
145
84
 
146
85
!     TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU)
149
88
!            --- tableau de travail
150
89
!===============================================================================
151
90
 
 
91
!===============================================================================
 
92
! Module files
 
93
!===============================================================================
 
94
 
 
95
use paramx
 
96
use numvar
 
97
use optcal
 
98
use cstphy
 
99
use cstnum
 
100
use entsor
 
101
use parall
 
102
use period
 
103
use ppppar
 
104
use ppthch
 
105
use coincl
 
106
use cpincl
 
107
use ppincl
 
108
use mesh
 
109
 
 
110
!===============================================================================
 
111
 
152
112
implicit none
153
113
 
154
 
!===============================================================================
155
 
! Common blocks
156
 
!===============================================================================
157
 
 
158
 
include "paramx.h"
159
 
include "numvar.h"
160
 
include "optcal.h"
161
 
include "cstphy.h"
162
 
include "cstnum.h"
163
 
include "entsor.h"
164
 
include "parall.h"
165
 
include "period.h"
166
 
include "ppppar.h"
167
 
include "ppthch.h"
168
 
include "coincl.h"
169
 
include "cpincl.h"
170
 
include "ppincl.h"
171
 
 
172
 
!===============================================================================
173
 
 
174
 
integer          idbia0 , idbra0
175
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
176
 
integer          nfml   , nprfml
177
 
integer          nnod   , lndfac , lndfbr , ncelbr
178
 
integer          nvar   , nscal  , nphas
179
 
integer          nideve , nrdeve , nituse , nrtuse
180
 
 
181
 
integer          ifacel(2,nfac) , ifabor(nfabor)
182
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
183
 
integer          iprfml(nfml,nprfml), maxelt, lstelt(maxelt)
184
 
integer          ipnfac(nfac+1), nodfac(lndfac)
185
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
186
 
integer          idevel(nideve), ituser(nituse), ia(*)
187
 
 
188
 
double precision xyzcen(ndim,ncelet)
189
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
190
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
191
 
double precision xyznod(ndim,nnod), volume(ncelet)
 
114
integer          nvar   , nscal
 
115
 
 
116
 
192
117
double precision dt(ncelet), rtp(ncelet,*), propce(ncelet,*)
193
118
double precision propfa(nfac,*), propfb(nfabor,*)
194
119
double precision coefa(nfabor,*), coefb(nfabor,*)
195
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
196
120
 
197
121
 
198
122
! Local variables
199
123
 
200
124
character*80     chaine
201
 
integer          idebia, idebra
202
 
integer          iel, mode, igg, iphas, izone
203
 
integer          iscal, ivar, ii, idimte, itenso
 
125
integer          iel, mode, igg, izone
 
126
integer          iscal, ivar, ii
204
127
double precision hinit, coefg(ngazgm), hair, tinitk
205
128
double precision sommqf, sommqt, sommq, tentm, fmelm
206
129
double precision valmax, valmin, xkent, xeent, d2s3
219
142
 
220
143
ipass = ipass + 1
221
144
 
222
 
idebia = idbia0
223
 
idebra = idbra0
224
145
 
225
146
do igg = 1, ngazgm
226
147
  coefg(igg) = zero
227
148
enddo
228
149
 
229
 
iphas  = 1
230
 
 
231
150
d2s3 = 2.d0/3.d0
232
151
 
233
152
!===============================================================================
243
162
  if ( ipass.eq.1 ) then
244
163
 
245
164
! ----- Temperature du melange : air a TINITK
246
 
    tinitk = t0(iphas)
 
165
    tinitk = t0
247
166
 
248
167
! ----- Enthalpie de l'air a TINITK
249
168
    if ( ippmod(icoebu).eq.1 .or. ippmod(icoebu).eq.3 ) then
273
192
 
274
193
! ---- TURBULENCE
275
194
 
276
 
      if (itytur(iphas).eq.2) then
277
 
 
278
 
        rtp(iel,ik(iphas))  = xkent
279
 
        rtp(iel,iep(iphas)) = xeent
280
 
 
281
 
      elseif (itytur(iphas).eq.3) then
282
 
 
283
 
        rtp(iel,ir11(iphas)) = d2s3*xkent
284
 
        rtp(iel,ir22(iphas)) = d2s3*xkent
285
 
        rtp(iel,ir33(iphas)) = d2s3*xkent
286
 
        rtp(iel,ir12(iphas)) = 0.d0
287
 
        rtp(iel,ir13(iphas)) = 0.d0
288
 
        rtp(iel,ir23(iphas)) = 0.d0
289
 
        rtp(iel,iep(iphas))  = xeent
290
 
 
291
 
      elseif (iturb(iphas).eq.50) then
292
 
 
293
 
        rtp(iel,ik(iphas))   = xkent
294
 
        rtp(iel,iep(iphas))  = xeent
295
 
        rtp(iel,iphi(iphas)) = d2s3
296
 
        rtp(iel,ifb(iphas))  = 0.d0
297
 
 
298
 
      elseif (iturb(iphas).eq.60) then
299
 
 
300
 
        rtp(iel,ik(iphas))   = xkent
301
 
        rtp(iel,iomg(iphas)) = xeent/cmu/xkent
 
195
      if (itytur.eq.2) then
 
196
 
 
197
        rtp(iel,ik)  = xkent
 
198
        rtp(iel,iep) = xeent
 
199
 
 
200
      elseif (itytur.eq.3) then
 
201
 
 
202
        rtp(iel,ir11) = d2s3*xkent
 
203
        rtp(iel,ir22) = d2s3*xkent
 
204
        rtp(iel,ir33) = d2s3*xkent
 
205
        rtp(iel,ir12) = 0.d0
 
206
        rtp(iel,ir13) = 0.d0
 
207
        rtp(iel,ir23) = 0.d0
 
208
        rtp(iel,iep)  = xeent
 
209
 
 
210
      elseif (iturb.eq.50) then
 
211
 
 
212
        rtp(iel,ik)   = xkent
 
213
        rtp(iel,iep)  = xeent
 
214
        rtp(iel,iphi) = d2s3
 
215
        rtp(iel,ifb)  = 0.d0
 
216
 
 
217
      elseif (iturb.eq.60) then
 
218
 
 
219
        rtp(iel,ik)   = xkent
 
220
        rtp(iel,iomg) = xeent/cmu/xkent
 
221
 
 
222
      elseif(iturb.eq.70) then
 
223
 
 
224
        rtp(iel,inusa) = cmu*xkent**2/xeent
302
225
 
303
226
      endif
304
227
 
340
263
      tentm = sommqt / sommq
341
264
    else
342
265
      fmelm = zero
343
 
      tentm = t0(iphas)
 
266
      tentm = t0
344
267
    endif
345
268
 
346
269
! ----- Enthalpie du melange HINIT
381
304
 
382
305
    call usebui                                                   &
383
306
    !==========
384
 
 ( idebia , idebra ,                                              &
385
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
386
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
387
 
   nvar   , nscal  , nphas  ,                                     &
388
 
   nideve , nrdeve , nituse , nrtuse ,                            &
389
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , &
390
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
391
 
   idevel , ituser , ia     ,                                     &
392
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
393
 
   dt     , rtp    , propce , propfa , propfb , coefa  , coefb  , &
394
 
   rdevel , rtuser , ra     )
 
307
 ( nvar   , nscal  ,                                              &
 
308
   dt     , rtp    , propce , propfa , propfb , coefa  , coefb  )
395
309
 
396
310
! ----- En periodique et en parallele,
397
311
!       il faut echanger ces initialisations (qui sont en fait dans RTPA)
398
312
 
399
 
!     Parallele
400
 
    if(irangp.ge.0) then
401
 
      call parcom(rtp(1,isca(iygfm)))
402
 
      !==========
403
 
      if ( ippmod(icoebu).eq.2 .or. ippmod(icoebu).eq.3 ) then
404
 
        call parcom(rtp(1,isca(ifm  )))
405
 
        !==========
406
 
      endif
407
 
      if ( ippmod(icoebu).eq.1 .or. ippmod(icoebu).eq.3 ) then
408
 
        call parcom(rtp(1,isca(ihm  )))
409
 
        !==========
410
 
      endif
411
 
    endif
412
 
 
413
 
!     Periodique
414
 
    if(iperio.eq.1) then
415
 
      idimte = 0
416
 
      itenso = 0
417
 
      ivar   = isca(iygfm)
418
 
      call percom                                                 &
419
 
      !==========
420
 
      ( idimte , itenso ,                                         &
421
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar),                    &
422
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar),                    &
423
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar))
424
 
      if ( ippmod(icoebu).eq.2 .or. ippmod(icoebu).eq.3 ) then
425
 
        idimte = 0
426
 
        itenso = 0
427
 
        ivar   = isca(ifm  )
428
 
        call percom                                               &
429
 
        !==========
430
 
      ( idimte , itenso ,                                         &
431
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar),                    &
432
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar),                    &
433
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar))
434
 
      endif
435
 
      if ( ippmod(icoebu).eq.1 .or. ippmod(icoebu).eq.3 ) then
436
 
        idimte = 0
437
 
        itenso = 0
438
 
        ivar   = isca(ihm  )
439
 
        call percom                                               &
440
 
        !==========
441
 
      ( idimte , itenso ,                                         &
442
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar),                    &
443
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar),                    &
444
 
        rtp(1,ivar), rtp(1,ivar), rtp(1,ivar))
 
313
    if (irangp.ge.0.or.iperio.eq.1) then
 
314
      call synsca(rtp(1,isca(iygfm)))
 
315
      !==========
 
316
      if ( ippmod(icoebu).eq.2 .or. ippmod(icoebu).eq.3 ) then
 
317
        call synsca(rtp(1,isca(ifm)))
 
318
        !==========
 
319
      endif
 
320
      if ( ippmod(icoebu).eq.1 .or. ippmod(icoebu).eq.3 ) then
 
321
        call synsca(rtp(1,isca(ihm)))
 
322
        !==========
445
323
      endif
446
324
    endif
447
325