1!-------------------------------------------------------------------------------
2
3!VERS
4
5! This file is part of Code_Saturne, a general-purpose CFD tool.
6!
7! Copyright (C) 1998-2021 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.
22
23!-------------------------------------------------------------------------------
24
25!===============================================================================
26! Function:
27! ---------
28
29!> \file cs_user_boundary_conditions-atmospheric.f90
30!>
31!> Atmospheric example of cs_user_boundary_conditions.f90 subroutine
32!>
33!-------------------------------------------------------------------------------
34
35!-------------------------------------------------------------------------------
36! Arguments
37!______________________________________________________________________________.
38!  mode           name          role                                           !
39!______________________________________________________________________________!
40!> \param[in]     nvar          total number of variables
41!> \param[in]     nscal         total number of scalars
42!> \param[out]    icodcl        boundary condition code:
43!>                               - 1 Dirichlet
44!>                               - 2 Radiative outlet
45!>                               - 3 Neumann
46!>                               - 4 sliding and
47!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
48!>                               - 5 smooth wall and
49!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
50!>                               - 6 rough wall and
51!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
52!>                               - 9 free inlet/outlet
53!>                                 (input mass flux blocked to 0)
54!>                               - 13 Dirichlet for the advection operator and
55!>                                    Neumann for the diffusion operator
56!> \param[in]     itrifb        indirection for boundary faces ordering
57!> \param[in,out] itypfb        boundary face types
58!> \param[in,out] izfppp        boundary face zone number
59!> \param[in]     dt            time step (per cell)
60!> \param[in,out] rcodcl        boundary condition values:
61!>                               - rcodcl(1) value of the dirichlet
62!>                               - rcodcl(2) value of the exterior exchange
63!>                                 coefficient (infinite if no exchange)
64!>                               - rcodcl(3) value flux density
65!>                                 (negative if gain) in w/m2
66!>                                 -# for the velocity \f$ (\mu+\mu_T)
67!>                                    \gradt \, \vect{u} \cdot \vect{n}  \f$
68!>                                 -# for the pressure \f$ \Delta t
69!>                                    \grad P \cdot \vect{n}  \f$
70!>                                 -# for a scalar \f$ cp \left( K +
71!>                                     \dfrac{K_T}{\sigma_T} \right)
72!>                                     \grad T \cdot \vect{n} \f$
73!_______________________________________________________________________________
74
75subroutine cs_f_user_boundary_conditions &
76 ( nvar   , nscal  ,                                              &
77   icodcl , itrifb , itypfb , izfppp ,                            &
78   dt     ,                                                       &
79   rcodcl )
80
81!===============================================================================
82
83!===============================================================================
84! Module files
85!===============================================================================
86
87use paramx
88use numvar
89use optcal
90use cstphy
91use cstnum
92use entsor
93use parall
94use period
95use ppppar
96use ppthch
97use coincl
98use cpincl
99use ppincl
100use ppcpfu
101use atchem
102use atincl
103use atsoil
104use ctincl
105use cs_fuel_incl
106use mesh
107use field
108use turbomachinery
109use iso_c_binding
110use cs_c_bindings
111
112!===============================================================================
113
114implicit none
115
116! Arguments
117
118integer          nvar   , nscal
119
120integer          icodcl(nfabor,nvar)
121integer          itrifb(nfabor), itypfb(nfabor)
122integer          izfppp(nfabor)
123
124double precision dt(ncelet)
125double precision rcodcl(nfabor,nvar,3)
126
127! Local variables
128
129!< [loc_var_dec]
130integer          ifac, ii
131integer          izone
132integer          ilelt, nlelt
133integer          f_id_rough, f_id_t_rough
134double precision d2s3
135double precision zref, xuref
136double precision ustar, rugd, rugt
137double precision zent, xuent, xvent
138double precision xkent, xeent
139
140integer, allocatable, dimension(:) :: lstelt
141double precision, dimension(:), pointer :: bpro_roughness
142double precision, dimension(:), pointer :: bpro_roughness_t
143!< [loc_var_dec]
144
145!===============================================================================
146
147!===============================================================================
148! Initialization
149!===============================================================================
150
151allocate(lstelt(nfabor))  ! temporary array for boundary faces selection
152
153d2s3 = 2.d0/3.d0
154
155! Paremeters for the analytical rough wall law (neutral)
156zref = 10.d0
157xuref = 10.d0
158rugd = 0.1d0
159rugt = rugd
160
161!===============================================================================
162! Assign boundary conditions to boundary faces here
163
164! For each subset:
165! - use selection criteria to filter boundary faces of a given subset
166! - loop on faces from a subset
167!   - set the boundary condition for each face
168!===============================================================================
169
170! --- For boundary faces of color 11,
171!       assign an inlet boundary condition for all phases prescribed from the
172!       meteo profile with automatic choice between inlet/ outlet according to
173!       the meteo profile
174
175!< [example_1]
176call getfbr('11',nlelt,lstelt)
177
178! Get a new zone number (1 <= izone <= nozppm)
179izone = maxval(izfppp) + 1
180
181do ilelt = 1, nlelt
182
183  ifac = lstelt(ilelt)
184
185  ! - Zone to which the face belongs
186  izfppp(ifac) = izone
187
188  ! - Boundary conditions are prescribed from the meteo profile
189  iprofm(izone) = 1
190
191  ! - Chemical boundary conditions are prescribed from the chemistry profile
192  iprofc(izone) = 1
193
194  ! - boundary condition type can be set to ientre or i_convective_inlet
195
196  itypfb(ifac) = ientre
197
198  ! - automatic determination of type (inlet/outlet) according to sign of
199  !   mass flux
200
201  iautom(ifac) = 1
202
203enddo
204!< [example_1]
205
206
207! ---For boundary faces of color 21,
208!     assign an inlet boundary condition for all phases prescribed from the
209!     meteo profile
210
211!< [example_2]
212call getfbr('21',nlelt,lstelt)
213
214! Get a new zone number (1 <= izone <= nozppm)
215izone = maxval(izfppp) + 1
216
217do ilelt = 1, nlelt
218
219  ifac = lstelt(ilelt)
220
221  ! - Zone to which the face belongs
222  izfppp(ifac) = izone
223
224  ! - Boundary conditions are prescribed from the meteo profile
225  iprofm(izone) = 1
226
227  ! - Chemical boundary conditions are prescribed from the chemistry profile
228  iprofc(izone) = 1
229
230  ! - Assign inlet boundary conditions
231  itypfb(ifac) = ientre
232
233enddo
234!< [example_2]
235
236! --- For boundary faces of color 31,
237!       assign an inlet boundary condition for all phases prescribed from the
238!       meteo profile except for dynamical variables which are prescribed with
239!       a rough log law.
240
241!< [example_3]
242call getfbr('31',nlelt,lstelt)
243
244! Get a new zone number (1 <= izone <= nozppm)
245izone = maxval(izfppp) + 1
246
247do ilelt = 1, nlelt
248
249  ifac = lstelt(ilelt)
250
251  ! - Zone to which the face belongs
252  izfppp(ifac) = izone
253
254  ! - Boundary conditions are prescribed from the meteo profile
255  iprofm(izone) = 1
256
257  ! - Chemical boundary conditions are prescribed from the chemistry profile
258  iprofc(izone) = 1
259
260  ! - Dynamical variables are prescribed with a rough log law
261  zent=cdgfbo(3,ifac)
262
263  ustar = xkappa*xuref/log((zref+rugd)/rugd)
264  xuent = ustar/xkappa*log((zent+rugd)/rugd)
265  xvent = 0.d0
266  xkent = ustar**2/sqrt(cmu)
267  xeent = ustar**3/xkappa/(zent+rugd)
268
269  itypfb(ifac) = ientre
270
271  rcodcl(ifac,iu,1) = xuent
272  rcodcl(ifac,iv,1) = xvent
273  rcodcl(ifac,iw,1) = 0.d0
274
275  ! itytur is a flag equal to iturb/10
276  if    (itytur.eq.2) then
277
278    rcodcl(ifac,ik,1)  = xkent
279    rcodcl(ifac,iep,1) = xeent
280
281  elseif(itytur.eq.3) then
282
283    rcodcl(ifac,ir11,1) = d2s3*xkent
284    rcodcl(ifac,ir22,1) = d2s3*xkent
285    rcodcl(ifac,ir33,1) = d2s3*xkent
286    rcodcl(ifac,ir12,1) = 0.d0
287    rcodcl(ifac,ir13,1) = 0.d0
288    rcodcl(ifac,ir23,1) = 0.d0
289    rcodcl(ifac,iep,1)  = xeent
290
291  elseif(iturb.eq.50) then
292
293    rcodcl(ifac,ik,1)   = xkent
294    rcodcl(ifac,iep,1)  = xeent
295    rcodcl(ifac,iphi,1) = d2s3
296    rcodcl(ifac,ifb,1)  = 0.d0
297
298  elseif(iturb.eq.60) then
299
300    rcodcl(ifac,ik,1)   = xkent
301    rcodcl(ifac,iomg,1) = xeent/cmu/xkent
302
303  elseif(iturb.eq.70) then
304
305    rcodcl(ifac,inusa,1) = cmu*xkent**2/xeent
306
307  endif
308
309enddo
310
311!< [example_3]
312! --- Prescribe at boundary faces of color '12' an outlet for all phases
313call getfbr('12', nlelt, lstelt)
314
315! Get a new zone number (1 <= izone <= nozppm)
316izone = maxval(izfppp) + 1
317
318do ilelt = 1, nlelt
319
320  ifac = lstelt(ilelt)
321
322  ! - Zone to which the zone belongs
323  izfppp(ifac) = izone
324
325  ! Outlet: zero flux for velocity and temperature, prescribed pressure
326  !         Note that the pressure will be set to P0 at the first
327  !         free outlet face (isolib)
328
329  itypfb(ifac)   = isolib
330
331enddo
332
333! --- Prescribe at boundary faces of color 15 a rough wall
334!< [example_4]
335call field_get_id_try("boundary_roughness", f_id_rough)
336if (f_id_rough.ge.0) call field_get_val_s(f_id_rough, bpro_roughness)
337call field_get_id_try("boundary_thermal_roughness", f_id_t_rough)
338if (f_id_t_rough.ge.0) call field_get_val_s(f_id_t_rough, bpro_roughness_t)
339
340call getfbr('15', nlelt, lstelt)
341
342! Get a new zone number (1 <= izone <= nozppm)
343izone = maxval(izfppp) + 1
344
345do ilelt = 1, nlelt
346
347  ifac = lstelt(ilelt)
348
349  ! Wall: zero flow (zero flux for pressure)
350  !       rough friction for velocities (+ turbulent variables)
351
352  ! - Zone to which the zone belongs
353  izfppp(ifac) = izone
354
355  itypfb(ifac)   = iparug
356
357  ! Roughness for velocity: rugd
358  bpro_roughness(ifac) = rugd
359
360  ! Roughness for scalars (if required)
361  if (f_id_rough.ge.0) &
362    bpro_roughness_t(ifac) = 0.01d0
363
364  ! By default zero flux for scalars
365
366enddo
367!< [example_4]
368
369! --- Prescribe at boundary faces of color 4 a symmetry for all phases
370!< [example_5]
371call getfbr('4', nlelt, lstelt)
372
373! Get a new zone number (1 <= izone <= nozppm)
374izone = maxval(izfppp) + 1
375
376do ilelt = 1, nlelt
377
378  ifac = lstelt(ilelt)
379
380  ! - Zone to which the zone belongs
381  izfppp(ifac) = izone
382
383  itypfb(ifac)   = isymet
384
385enddo
386!< [example_5]
387
388!--------
389! Formats
390!--------
391
392!----
393! End
394!----
395
396deallocate(lstelt)  ! temporary array for boundary faces selection
397
398return
399end subroutine cs_f_user_boundary_conditions
400