← Back to branch summary

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

« back to all changes in this revision

Viewing changes to users/atmo/usatcl.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:
users/atmo/usatcl.f90
2
2
 
3
3
!VERS
4
4
 
5
 
 
6
 
!     This file is part of the Code_Saturne Kernel, element of the
7
 
!     Code_Saturne CFD tool.
8
 
 
9
 
!     Copyright (C) 1998-2009 EDF S.A., France
10
 
 
11
 
!     contact: saturne-support@edf.fr
12
 
 
13
 
!     The Code_Saturne Kernel is free software; you can redistribute it
14
 
!     and/or modify it under the terms of the GNU General Public License
15
 
!     as published by the Free Software Foundation; either version 2 of
16
 
!     the License, or (at your option) any later version.
17
 
 
18
 
!     The Code_Saturne Kernel is distributed in the hope that it will be
19
 
!     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
20
 
!     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21
 
!     GNU General Public License for more details.
22
 
 
23
 
!     You should have received a copy of the GNU General Public License
24
 
!     along with the Code_Saturne Kernel; if not, write to the
25
 
!     Free Software Foundation, Inc.,
26
 
!     51 Franklin St, Fifth Floor,
27
 
!     Boston, MA  02110-1301  USA
 
5
! This file is part of Code_Saturne, a general-purpose CFD tool.
 
6
!
 
7
! Copyright (C) 1998-2011 EDF S.A.
 
8
!
 
9
! This program is free software; you can redistribute it and/or modify it under
 
10
! the terms of the GNU General Public License as published by the Free Software
 
11
! Foundation; either version 2 of the License, or (at your option) any later
 
12
! version.
 
13
!
 
14
! This program is distributed in the hope that it will be useful, but WITHOUT
 
15
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 
16
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 
17
! details.
 
18
!
 
19
! You should have received a copy of the GNU General Public License along with
 
20
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
 
21
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
28
22
 
29
23
!-------------------------------------------------------------------------------
30
24
 
31
25
subroutine usatcl &
32
26
!================
33
27
 
34
 
 ( idbia0 , idbra0 ,                                              &
35
 
   ndim   , ncelet , ncel   , nfac   , nfabor , nfml   , nprfml , &
36
 
   nnod   , lndfac , lndfbr , ncelbr ,                            &
37
 
   nvar   , nscal  , nphas  ,                                     &
38
 
   nideve , nrdeve , nituse , nrtuse ,                            &
39
 
   ifacel , ifabor , ifmfbr , ifmcel , iprfml , maxelt , lstelt , &
40
 
   ipnfac , nodfac , ipnfbr , nodfbr ,                            &
 
28
 ( nvar   , nscal  ,                                              &
41
29
   icodcl , itrifb , itypfb , izfppp ,                            &
42
 
   idevel , ituser , ia     ,                                     &
43
 
   xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , &
44
30
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
45
 
   coefa  , coefb  , rcodcl ,                                     &
46
 
   w1     , w2     , w3     , w4     , w5     , w6     , coefu  , &
47
 
   rdevel , rtuser , ra     )
 
31
   coefa  , coefb  , rcodcl )
48
32
 
49
33
!===============================================================================
50
34
! Purpose:
81
65
!__________________.____._____.________________________________________________.
82
66
! name             !type!mode ! role                                           !
83
67
!__________________!____!_____!________________________________________________!
84
 
! idbia0           ! i  ! <-- ! number of first free position in ia            !
85
 
! idbra0           ! i  ! <-- ! number of first free position in ra            !
86
 
! ndim             ! i  ! <-- ! spatial dimension                              !
87
 
! ncelet           ! i  ! <-- ! number of extended (real + ghost) cells        !
88
 
! ncel             ! i  ! <-- ! number of cells                                !
89
 
! nfac             ! i  ! <-- ! number of interior faces                       !
90
 
! nfabor           ! i  ! <-- ! number of boundary faces                       !
91
 
! nfml             ! i  ! <-- ! number of families (group classes)             !
92
 
! nprfml           ! i  ! <-- ! number of properties per family (group class)  !
93
 
! nnod             ! i  ! <-- ! number of vertices                             !
94
 
! lndfac           ! i  ! <-- ! size of nodfac indexed array                   !
95
 
! lndfbr           ! i  ! <-- ! size of nodfbr indexed array                   !
96
 
! ncelbr           ! i  ! <-- ! number of cells with faces on boundary         !
97
68
! nvar             ! i  ! <-- ! total number of variables                      !
98
69
! nscal            ! i  ! <-- ! total number of scalars                        !
99
 
! nphas            ! i  ! <-- ! number of phases                               !
100
 
! nideve, nrdeve   ! i  ! <-- ! sizes of idevel and rdevel arrays              !
101
 
! nituse, nrtuse   ! i  ! <-- ! sizes of ituser and rtuser arrays              !
102
 
! ifacel(2, nfac)  ! ia ! <-- ! interior faces -> cells connectivity           !
103
 
! ifabor(nfabor)   ! ia ! <-- ! boundary faces -> cells connectivity           !
104
 
! ifmfbr(nfabor)   ! ia ! <-- ! boundary face family numbers                   !
105
 
! ifmcel(ncelet)   ! ia ! <-- ! cell family numbers                            !
106
 
! iprfml           ! ia ! <-- ! property numbers per family                    !
107
 
!  (nfml, nprfml)  !    !     !                                                !
108
 
! maxelt           ! i  ! <-- ! max number of cells and faces (int/boundary)   !
109
 
! lstelt(maxelt)   ! ia ! --- ! work array                                     !
110
 
! ipnfac(nfac+1)   ! ia ! <-- ! interior faces -> vertices index (optional)    !
111
 
! nodfac(lndfac)   ! ia ! <-- ! interior faces -> vertices list (optional)     !
112
 
! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional)    !
113
 
! nodfbr(lndfbr)   ! ia ! <-- ! boundary faces -> vertices list (optional)     !
114
70
! icodcl           ! ia ! --> ! boundary condition code                        !
115
71
!  (nfabor, nvar)  !    !     ! = 1  -> Dirichlet                              !
116
72
!                  !    !     ! = 2  -> flux density                           !
120
76
!                  !    !     ! = 9  -> free inlet/outlet (velocity)           !
121
77
!                  !    !     !         inflowing possibly blocked             !
122
78
! itrifb           ! ia ! <-- ! indirection for boundary faces ordering        !
123
 
!  (nfabor, nphas) !    !     !                                                !
124
79
! itypfb           ! ia ! --> ! boundary face types                            !
125
 
!  (nfabor, nphas) !    !     !                                                !
126
80
! izfppp(nfabor)   ! te ! --> ! boundary face zone number                      !
127
 
! idevel(nideve)   ! ia ! <-> ! integer work array for temporary development   !
128
 
! ituser(nituse)   ! ia ! <-> ! user-reserved integer work array               !
129
 
! ia(*)            ! ia ! --- ! main integer work array                        !
130
 
! xyzcen           ! ra ! <-- ! cell centers                                   !
131
 
!  (ndim, ncelet)  !    !     !                                                !
132
 
! surfac           ! ra ! <-- ! interior faces surface vectors                 !
133
 
!  (ndim, nfac)    !    !     !                                                !
134
 
! surfbo           ! ra ! <-- ! boundary faces surface vectors                 !
135
 
!  (ndim, nfabor)  !    !     !                                                !
136
 
! cdgfac           ! ra ! <-- ! interior faces centers of gravity              !
137
 
!  (ndim, nfac)    !    !     !                                                !
138
 
! cdgfbo           ! ra ! <-- ! boundary faces centers of gravity              !
139
 
!  (ndim, nfabor)  !    !     !                                                !
140
 
! xyznod           ! ra ! <-- ! vertex coordinates (optional)                  !
141
 
!  (ndim, nnod)    !    !     !                                                !
142
 
! volume(ncelet)   ! ra ! <-- ! cell volumes                                   !
143
81
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
144
82
! rtp, rtpa        ! ra ! <-- ! calculated variables at cell centers           !
145
83
!  (ncelet, *)     !    !     !  (at current and previous time steps)          !
158
96
!                  !    !     ! for velocities           ( vistl+visct)*gradu  !
159
97
!                  !    !     ! for pressure                         dt*gradp  !
160
98
!                  !    !     ! for scalars    cp*(viscls+visct/sigmas)*gradt  !
161
 
! w1,2,3,4,5,6     ! ra ! --- ! work arrays                                    !
162
 
!  (ncelet)        !    !     !  (computation of pressure gradient)            !
163
 
! coefu            ! ra ! --- ! work array                                     !
164
 
!  (nfabor, 3)     !    !     !  (computation of pressure gradient)            !
165
 
! rdevel(nrdeve)   ! ra ! <-> ! real work array for temporary development      !
166
 
! rtuser(nrtuse)   ! ra ! <-> ! user-reserved real work array                  !
167
 
! ra(*)            ! ra ! --- ! main real work array                           !
168
99
!__________________!____!_____!________________________________________________!
169
100
 
170
101
!     Type: i (integer), r (real), s (string), a (array), l (logical),
172
103
!     mode: <-- input, --> output, <-> modifies data, --- work array
173
104
!===============================================================================
174
105
 
 
106
!===============================================================================
 
107
! Module files
 
108
!===============================================================================
 
109
 
 
110
use paramx
 
111
use numvar
 
112
use optcal
 
113
use cstphy
 
114
use cstnum
 
115
use entsor
 
116
use parall
 
117
use period
 
118
use ihmpre
 
119
use ppppar
 
120
use atincl
 
121
use mesh
 
122
 
 
123
!===============================================================================
 
124
 
175
125
implicit none
176
126
 
177
 
!===============================================================================
178
 
! Common blocks
179
 
!===============================================================================
180
 
 
181
 
include "paramx.h"
182
 
include "pointe.h"
183
 
include "numvar.h"
184
 
include "optcal.h"
185
 
include "cstphy.h"
186
 
include "cstnum.h"
187
 
include "entsor.h"
188
 
include "parall.h"
189
 
include "period.h"
190
 
include "ihmpre.h"
191
 
include "ppppar.h"
192
 
include "atincl.h"
193
 
 
194
 
!===============================================================================
195
 
 
196
127
! Arguments
197
128
 
198
 
integer          idbia0 , idbra0
199
 
integer          ndim   , ncelet , ncel   , nfac   , nfabor
200
 
integer          nfml   , nprfml
201
 
integer          nnod   , lndfac , lndfbr , ncelbr
202
 
integer          nvar   , nscal  , nphas
203
 
integer          nideve , nrdeve , nituse , nrtuse
 
129
integer          nvar   , nscal
204
130
 
205
 
integer          ifacel(2,nfac) , ifabor(nfabor)
206
 
integer          ifmfbr(nfabor) , ifmcel(ncelet)
207
 
integer          iprfml(nfml,nprfml), maxelt, lstelt(maxelt)
208
 
integer          ipnfac(nfac+1), nodfac(lndfac)
209
 
integer          ipnfbr(nfabor+1), nodfbr(lndfbr)
210
131
integer          icodcl(nfabor,nvar)
211
 
integer          itrifb(nfabor,nphas), itypfb(nfabor,nphas)
 
132
integer          itrifb(nfabor), itypfb(nfabor)
212
133
integer          izfppp(nfabor)
213
 
integer          idevel(nideve), ituser(nituse), ia(*)
214
134
 
215
 
double precision xyzcen(ndim,ncelet)
216
 
double precision surfac(ndim,nfac), surfbo(ndim,nfabor)
217
 
double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor)
218
 
double precision xyznod(ndim,nnod), volume(ncelet)
219
135
double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
220
136
double precision propce(ncelet,*)
221
137
double precision propfa(nfac,*), propfb(nfabor,*)
222
138
double precision coefa(nfabor,*), coefb(nfabor,*)
223
139
double precision rcodcl(nfabor,nvar,3)
224
 
double precision w1(ncelet),w2(ncelet),w3(ncelet)
225
 
double precision w4(ncelet),w5(ncelet),w6(ncelet)
226
 
double precision coefu(nfabor,ndim)
227
 
double precision rdevel(nrdeve), rtuser(nrtuse), ra(*)
228
140
 
229
141
! Local variables
230
142
 
231
 
integer          idebia, idebra
232
 
integer          ifac, iel, ii, ivar, iphas
 
143
integer          ifac, iel, ii, ivar
233
144
integer          izone
234
145
integer          ilelt, nlelt
235
146
double precision uref2, d2s3
241
152
double precision xkent, xeent
242
153
double precision tpent
243
154
 
 
155
integer, allocatable, dimension(:) :: lstelt
244
156
 
245
157
!===============================================================================
246
158
 
282
194
! 1.  Initialization
283
195
!===============================================================================
284
196
 
285
 
idebia = idbia0
286
 
idebra = idbra0
 
197
! Allocate a temporary array for boundary faces selection
 
198
allocate(lstelt(nfabor))
 
199
 
287
200
 
288
201
d2s3 = 2.d0/3.d0
289
202
 
342
255
  iprofm(izone) = 1
343
256
 
344
257
!     - Assign inlet boundary conditions
345
 
  do iphas = 1, nphas
346
 
    itypfb(ifac,iphas) = ientre
347
 
  enddo
 
258
  itypfb(ifac) = ientre
348
259
 
349
260
enddo
350
261
 
377
288
  xkent=ustar**2/sqrt(cmu)
378
289
  xeent=ustar**3/xkappa/(zent+rugd)
379
290
 
380
 
  do iphas = 1, nphas
381
 
 
382
 
    itypfb(ifac,iphas) = ientre
383
 
 
384
 
    rcodcl(ifac,iu(iphas),1) = xuent
385
 
    rcodcl(ifac,iv(iphas),1) = xvent
386
 
    rcodcl(ifac,iw(iphas),1) = 0.d0
387
 
 
388
 
    ! itytur is a flag equal to iturb/10
389
 
    if    (itytur(iphas).eq.2) then
390
 
 
391
 
      rcodcl(ifac,ik(iphas),1)  = xkent
392
 
      rcodcl(ifac,iep(iphas),1) = xeent
393
 
 
394
 
    elseif(itytur(iphas).eq.3) then
395
 
 
396
 
      rcodcl(ifac,ir11(iphas),1) = d2s3*xkent
397
 
      rcodcl(ifac,ir22(iphas),1) = d2s3*xkent
398
 
      rcodcl(ifac,ir33(iphas),1) = d2s3*xkent
399
 
      rcodcl(ifac,ir12(iphas),1) = 0.d0
400
 
      rcodcl(ifac,ir13(iphas),1) = 0.d0
401
 
      rcodcl(ifac,ir23(iphas),1) = 0.d0
402
 
      rcodcl(ifac,iep(iphas),1)  = xeent
403
 
 
404
 
    elseif(iturb(iphas).eq.50) then
405
 
 
406
 
      rcodcl(ifac,ik(iphas),1)   = xkent
407
 
      rcodcl(ifac,iep(iphas),1)  = xeent
408
 
      rcodcl(ifac,iphi(iphas),1) = d2s3
409
 
      rcodcl(ifac,ifb(iphas),1)  = 0.d0
410
 
 
411
 
    elseif(iturb(iphas).eq.60) then
412
 
 
413
 
      rcodcl(ifac,ik(iphas),1)   = xkent
414
 
      rcodcl(ifac,iomg(iphas),1) = xeent/cmu/xkent
415
 
 
416
 
    endif
417
 
 
418
 
  enddo
 
291
  itypfb(ifac) = ientre
 
292
 
 
293
  rcodcl(ifac,iu,1) = xuent
 
294
  rcodcl(ifac,iv,1) = xvent
 
295
  rcodcl(ifac,iw,1) = 0.d0
 
296
 
 
297
  ! itytur is a flag equal to iturb/10
 
298
  if    (itytur.eq.2) then
 
299
 
 
300
    rcodcl(ifac,ik,1)  = xkent
 
301
    rcodcl(ifac,iep,1) = xeent
 
302
 
 
303
  elseif(itytur.eq.3) then
 
304
 
 
305
    rcodcl(ifac,ir11,1) = d2s3*xkent
 
306
    rcodcl(ifac,ir22,1) = d2s3*xkent
 
307
    rcodcl(ifac,ir33,1) = d2s3*xkent
 
308
    rcodcl(ifac,ir12,1) = 0.d0
 
309
    rcodcl(ifac,ir13,1) = 0.d0
 
310
    rcodcl(ifac,ir23,1) = 0.d0
 
311
    rcodcl(ifac,iep,1)  = xeent
 
312
 
 
313
  elseif(iturb.eq.50) then
 
314
 
 
315
    rcodcl(ifac,ik,1)   = xkent
 
316
    rcodcl(ifac,iep,1)  = xeent
 
317
    rcodcl(ifac,iphi,1) = d2s3
 
318
    rcodcl(ifac,ifb,1)  = 0.d0
 
319
 
 
320
  elseif(iturb.eq.60) then
 
321
 
 
322
    rcodcl(ifac,ik,1)   = xkent
 
323
    rcodcl(ifac,iomg,1) = xeent/cmu/xkent
 
324
 
 
325
  elseif(iturb.eq.70) then
 
326
 
 
327
    rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent
 
328
 
 
329
  endif
419
330
 
420
331
enddo
421
332
 
437
348
  !         Note that the pressure will be set to P0 at the first
438
349
  !         free outlet face (isolib)
439
350
 
440
 
  do iphas = 1, nphas
441
 
    itypfb(ifac,iphas)   = isolib
442
 
  enddo
 
351
  itypfb(ifac)   = isolib
443
352
 
444
353
enddo
445
354
 
461
370
!     - Zone to which the zone belongs
462
371
  izfppp(ifac) = izone
463
372
 
464
 
  do iphas = 1, nphas
465
 
    itypfb(ifac,iphas)   = iparug
466
 
 
467
 
!     Roughness for velocity: rugd
468
 
    rcodcl(ifac,iu(iphas),3) = rugd
469
 
 
470
 
!     Roughness for scalars (if required):
471
 
!   rcodcl(ifac,iv(iphas),3) = rugd
472
 
 
473
 
 
474
 
    if(iscalt(iphas).ne.-1) then
 
373
  itypfb(ifac)   = iparug
 
374
 
 
375
  !     Roughness for velocity: rugd
 
376
  rcodcl(ifac,iu,3) = rugd
 
377
 
 
378
  !     Roughness for scalars (if required):
 
379
  !   rcodcl(ifac,iv,3) = rugd
 
380
 
 
381
 
 
382
  if(iscalt.ne.-1) then
475
383
 
476
384
    ! If temperature prescribed to 20 with a rough wall law (scalar ii=1)
477
 
    ! (with thermal roughness specified in rcodcl(ifac,iv(iphas),3)) :
 
385
    ! (with thermal roughness specified in rcodcl(ifac,iv,3)) :
478
386
    ! ii = 1
479
387
    ! icodcl(ifac, isca(ii))    = 6
480
388
    ! rcodcl(ifac, isca(ii),1)  = 293.15d0
484
392
    ! icodcl(ifac, isca(ii))    = 3
485
393
    ! rcodcl(ifac, isca(ii), 3) = 4.D0
486
394
 
487
 
    endif
488
 
  enddo
 
395
  endif
489
396
enddo
490
397
 
491
398
! --- Prescribe at boundary faces of color 4 a symmetry for all phases
502
409
!     - Zone to which the zone belongs
503
410
  izfppp(ifac) = izone
504
411
 
505
 
  do iphas = 1, nphas
506
 
    itypfb(ifac,iphas)   = isymet
507
 
  enddo
508
 
 
 
412
  itypfb(ifac)   = isymet
509
413
 
510
414
enddo
511
415
 
517
421
! End
518
422
!----
519
423
 
 
424
! Deallocate the temporary array
 
425
deallocate(lstelt)
 
426
 
520
427
return
521
428
end subroutine