ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
kw.F90
Go to the documentation of this file.
1 module kw
2 
3 contains
4 
5  subroutine kw_block(resOnly)
6  !
7  ! kw solves the transport equations for the standard and
8  ! modified k-omega turbulence models in a decoupled manner
9  ! using a diagonal dominant adi-scheme.
10  !
11  use constants
12  use blockpointers, only: il, jl, kl
13  use inputtimespectral
14  use iteration
15  use utils, only: setpointers
16  use turbutils, only: kweddyviscosity
17  use turbbcroutines, only: bcturbtreatment, applyallturbbcthisblock
18  implicit none
19  !
20  ! Subroutine argument.
21  !
22  logical, intent(in) :: resOnly
23 
24  ! Set the arrays for the boundary condition treatment.
25 
26  call bcturbtreatment
27 
28  ! Solve the transport equations for k and omega.
29 
30  call kwsolve(resonly)
31 
32  ! The eddy viscosity and the boundary conditions are only
33  ! applied if an actual update has been computed in kwSolve.
34 
35  if (.not. resonly) then
36 
37  ! Compute the corresponding eddy viscosity.
38 
39  call kweddyviscosity(2, il, 2, jl, 2, kl)
40 
41  ! Set the halo values for the turbulent variables.
42  ! We are on the finest mesh, so the second layer of halo
43  ! cells must be computed as well.
44 
45  call applyallturbbcthisblock(.true.)
46 
47  end if
48 
49  end subroutine kw_block
50 
51  subroutine kwsolve(resOnly)
52  !
53  ! kwSolve solves the k-omega transport equations of both
54  ! the original and modified k-omega models
55  ! in a coupled manner using a diagonal dominant ADI-scheme.
56  !
57  use blockpointers
58  use constants
59  use flowvarrefstate
60  use inputiteration
61  use inputphysics
62  use paramturb
63  use turbmod, only: prod, dvt, sct, kwcd, sig1, sig2, vort
64  use turbutils, only: prodsmag2, prodwmag2, prodkatolaunder, &
65  turbadvection, unsteadyturbterm, tdia3, kwcdterm
66  use turbcurvefits, only: curvetupyp
67  implicit none
68  !
69  ! Subroutine arguments.
70  !
71  logical, intent(in) :: resOnly
72  !
73  ! Local variables.
74  !
75  integer(kind=intType) :: i, j, k, nn
76 
77  real(kind=realtype) :: rkwgam1
78  real(kind=realtype) :: rhoi, ss, spk, sdk
79  real(kind=realtype) :: voli, volmi, volpi
80  real(kind=realtype) :: xm, ym, zm, xp, yp, zp, xa, ya, za
81  real(kind=realtype) :: ttm, ttp, mulm, mulp, muem, muep
82  real(kind=realtype) :: c1m, c1p, c10, c2m, c2p, c20
83  real(kind=realtype) :: b1, b2, c1, c2, d1, d2
84  real(kind=realtype) :: qs, uu, um, up, factor, utau, rblank
85 
86  real(kind=realtype), dimension(itu1:itu2) :: tup
87 
88  real(kind=realtype), dimension(2:il, 2:jl, 2:kl, 2, 2) :: qq
89  real(kind=realtype), dimension(2, 2:max(il, jl, kl)) :: bb, dd, ff
90  real(kind=realtype), dimension(2, 2, 2:max(il, jl, kl)) :: cc
91 
92  real(kind=realtype), dimension(:, :, :), pointer :: ddw, ww, ddvt
93  real(kind=realtype), dimension(:, :), pointer :: rrlv
94  real(kind=realtype), dimension(:, :), pointer :: dd2wall
95 
96  logical, dimension(2:jl, 2:kl), target :: flagI2, flagIl
97  logical, dimension(2:il, 2:kl), target :: flagJ2, flagJl
98  logical, dimension(2:il, 2:jl), target :: flagK2, flagKl
99 
100  logical, dimension(:, :), pointer :: flag
101 
102  ! Set model constants
103 
104  rkwgam1 = rkwbeta1 / rkwbetas - rkwsigw1 * rkwk * rkwk / sqrt(rkwbetas)
105  sig1 = rkwsigk1
106  sig2 = rkwsigw1
107 
108  ! Set the pointer for dvt in dw, such that the code is more
109  ! readable. Also set the pointers for the production term,
110  ! vorticity and the cross diffusion term.
111 
112  dvt => scratch(1:, 1:, 1:, idvt:)
113  prod => scratch(1:, 1:, 1:, iprod)
114  vort => prod
115  kwcd => scratch(1:, 1:, 1:, icd)
116  !
117  ! Production term.
118  !
119  select case (turbprod)
120  case (strain)
121  call prodsmag2
122 
123  case (vorticity)
124  call prodwmag2
125 
126  case (katolaunder)
127  call prodkatolaunder
128 
129  end select
130  !
131  ! Source terms.
132  ! Determine the source term and its derivative w.r.t. k and
133  ! omega for all internal cells of the block.
134  !
135  do k = 2, kl
136  do j = 2, jl
137  do i = 2, il
138 
139  ! Compute the source terms for both the k and the omega
140  ! equation. Note that prod contains the unscaled production
141  ! term. Furthermore the production term of k is limited to
142  ! a certain times the destruction term.
143 
144  rhoi = one / w(i, j, k, irho)
145  ss = prod(i, j, k)
146  spk = rev(i, j, k) * ss * rhoi
147  sdk = rkwbetas * w(i, j, k, itu1) * w(i, j, k, itu2)
148  spk = min(spk, pklim * sdk)
149 
150  dvt(i, j, k, 1) = spk - sdk
151  dvt(i, j, k, 2) = rkwgam1 * ss - rkwbeta1 * w(i, j, k, itu2)**2
152 
153  ! Compute the source term jacobian. Note that only the
154  ! destruction terms are linearized to increase the diagonal
155  ! dominance of the matrix. Furthermore minus the source
156  ! term jacobian is stored.
157 
158  qq(i, j, k, 1, 1) = rkwbetas * w(i, j, k, itu2)
159  ! qq(i,j,k,1,2) = rkwBetas*w(i,j,k,itu1)
160  qq(i, j, k, 1, 2) = zero
161  qq(i, j, k, 2, 1) = zero
162  qq(i, j, k, 2, 2) = two * rkwbeta1 * w(i, j, k, itu2)
163 
164  end do
165  end do
166  end do
167  !
168  ! Compute the cross-diffusion term for the modified version of
169  ! the k-omega model. It should cure the free-stream dependency
170  ! of the original model.
171  !
172  if (turbmodel == komegamodified) then
173 
174  ! Compute the cross diffusion term.
175 
176  call kwcdterm
177 
178  ! Multiply it with the correct constant and take it only into
179  ! account if it is positive.
180 
181  do k = 2, kl
182  do j = 2, jl
183  do i = 2, il
184  dvt(i, j, k, 2) = dvt(i, j, k, 2) &
185  + rkwsigd1 * max(zero, kwcd(i, j, k))
186  end do
187  end do
188  end do
189 
190  end if
191  !
192  ! Advection and unsteady terms.
193  !
194  nn = itu1 - 1
195  call turbadvection(2_inttype, 2_inttype, nn, qq)
196 
197  call unsteadyturbterm(2_inttype, 2_inttype, nn, qq)
198  !
199  ! Viscous terms in k-direction.
200  !
201  do k = 2, kl
202  do j = 2, jl
203  do i = 2, il
204 
205  ! Compute the metrics in zeta-direction, i.e. along the
206  ! line k = constant.
207 
208  voli = one / vol(i, j, k)
209  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
210  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
211 
212  xm = sk(i, j, k - 1, 1) * volmi
213  ym = sk(i, j, k - 1, 2) * volmi
214  zm = sk(i, j, k - 1, 3) * volmi
215  xp = sk(i, j, k, 1) * volpi
216  yp = sk(i, j, k, 2) * volpi
217  zp = sk(i, j, k, 3) * volpi
218 
219  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
220  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
221  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
222  ttm = xm * xa + ym * ya + zm * za
223  ttp = xp * xa + yp * ya + zp * za
224 
225  ! Computation of the viscous terms in zeta-direction; note
226  ! that cross-derivatives are neglected, i.e. the mesh is
227  ! assumed to be orthogonal.
228  ! The second derivative in zeta-direction is constructed as
229  ! the central difference of the first order derivatives, i.e.
230  ! d^2/dzeta^2 = d/dzeta (d/dzeta k+1/2 - d/dzeta k-1/2).
231  ! In this way the metric as well as the varying viscosity
232  ! can be taken into account; the latter appears inside the
233  ! d/dzeta derivative. The whole term is divided by rho to
234  ! obtain the diffusion term for k and omega.
235 
236  rhoi = one / w(i, j, k, irho)
237  mulm = half * (rlv(i, j, k - 1) + rlv(i, j, k))
238  mulp = half * (rlv(i, j, k + 1) + rlv(i, j, k))
239  muem = half * (rev(i, j, k - 1) + rev(i, j, k))
240  muep = half * (rev(i, j, k + 1) + rev(i, j, k))
241 
242  c1m = ttm * (mulm + sig1 * muem) * rhoi
243  c1p = ttp * (mulp + sig1 * muep) * rhoi
244  c10 = c1m + c1p
245 
246  c2m = ttm * (mulm + sig2 * muem) * rhoi
247  c2p = ttp * (mulp + sig2 * muep) * rhoi
248  c20 = c2m + c2p
249 
250  ! Update the residual for this cell and store the possible
251  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
252 
253  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i, j, k - 1, itu1) &
254  - c10 * w(i, j, k, itu1) + c1p * w(i, j, k + 1, itu1)
255  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i, j, k - 1, itu2) &
256  - c20 * w(i, j, k, itu2) + c2p * w(i, j, k + 1, itu2)
257 
258  b1 = -c1m
259  c1 = c10
260  d1 = -c1p
261 
262  b2 = -c2m
263  c2 = c20
264  d2 = -c2p
265 
266  ! Update the central jacobian. For nonboundary cells this
267  ! is simply c1 and c2. For boundary cells this is slightly
268  ! more complicated, because the boundary conditions are
269  ! treated implicitly and the off-diagonal terms b1, b2 and
270  ! d1, d2 must be taken into account.
271  ! The boundary conditions are only treated implicitly if
272  ! the diagonal dominance of the matrix is increased.
273 
274  if (k == 2) then
275  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
276  - b1 * max(bmtk1(i, j, itu1, itu1), zero)
277  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmtk1(i, j, itu1, itu2)
278  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmtk1(i, j, itu2, itu1)
279  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
280  - b2 * max(bmtk1(i, j, itu2, itu2), zero)
281  else if (k == kl) then
282  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
283  - d1 * max(bmtk2(i, j, itu1, itu1), zero)
284  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmtk2(i, j, itu1, itu2)
285  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmtk2(i, j, itu2, itu1)
286  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
287  - d2 * max(bmtk2(i, j, itu2, itu2), zero)
288  else
289  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
290  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
291  end if
292 
293  end do
294  end do
295  end do
296  !
297  ! Viscous terms in j-direction.
298  !
299  do k = 2, kl
300  do j = 2, jl
301  do i = 2, il
302 
303  ! Compute the metrics in eta-direction, i.e. along the
304  ! line j = constant.
305 
306  voli = one / vol(i, j, k)
307  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
308  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
309 
310  xm = sj(i, j - 1, k, 1) * volmi
311  ym = sj(i, j - 1, k, 2) * volmi
312  zm = sj(i, j - 1, k, 3) * volmi
313  xp = sj(i, j, k, 1) * volpi
314  yp = sj(i, j, k, 2) * volpi
315  zp = sj(i, j, k, 3) * volpi
316 
317  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
318  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
319  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
320  ttm = xm * xa + ym * ya + zm * za
321  ttp = xp * xa + yp * ya + zp * za
322 
323  ! Computation of the viscous terms in eta-direction; note
324  ! that cross-derivatives are neglected, i.e. the mesh is
325  ! assumed to be orthogonal.
326  ! The second derivative in eta-direction is constructed as
327  ! the central difference of the first order derivatives, i.e.
328  ! d^2/deta^2 = d/deta (d/deta j+1/2 - d/deta j-1/2).
329  ! In this way the metric as well as the varying viscosity
330  ! can be taken into account; the latter appears inside the
331  ! d/deta derivative. The whole term is divided by rho to
332  ! obtain the diffusion term for k and omega.
333 
334  rhoi = one / w(i, j, k, irho)
335  mulm = half * (rlv(i, j - 1, k) + rlv(i, j, k))
336  mulp = half * (rlv(i, j + 1, k) + rlv(i, j, k))
337  muem = half * (rev(i, j - 1, k) + rev(i, j, k))
338  muep = half * (rev(i, j + 1, k) + rev(i, j, k))
339 
340  c1m = ttm * (mulm + sig1 * muem) * rhoi
341  c1p = ttp * (mulp + sig1 * muep) * rhoi
342  c10 = c1m + c1p
343 
344  c2m = ttm * (mulm + sig2 * muem) * rhoi
345  c2p = ttp * (mulp + sig2 * muep) * rhoi
346  c20 = c2m + c2p
347 
348  ! Update the residual for this cell and store the possible
349  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
350 
351  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i, j - 1, k, itu1) &
352  - c10 * w(i, j, k, itu1) + c1p * w(i, j + 1, k, itu1)
353  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i, j - 1, k, itu2) &
354  - c20 * w(i, j, k, itu2) + c2p * w(i, j + 1, k, itu2)
355 
356  b1 = -c1m
357  c1 = c10
358  d1 = -c1p
359 
360  b2 = -c2m
361  c2 = c20
362  d2 = -c2p
363 
364  ! Update the central jacobian. For nonboundary cells this
365  ! is simply c1 and c2. For boundary cells this is slightly
366  ! more complicated, because the boundary conditions are
367  ! treated implicitly and the off-diagonal terms b1, b2 and
368  ! d1, d2 must be taken into account.
369  ! The boundary conditions are only treated implicitly if
370  ! the diagonal dominance of the matrix is increased.
371 
372  if (j == 2) then
373  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
374  - b1 * max(bmtj1(i, k, itu1, itu1), zero)
375  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmtj1(i, k, itu1, itu2)
376  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmtj1(i, k, itu2, itu1)
377  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
378  - b2 * max(bmtj1(i, k, itu2, itu2), zero)
379  else if (j == jl) then
380  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
381  - d1 * max(bmtj2(i, k, itu1, itu1), zero)
382  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmtj2(i, k, itu1, itu2)
383  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmtj2(i, k, itu2, itu1)
384  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
385  - d2 * max(bmtj2(i, k, itu2, itu2), zero)
386  else
387  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
388  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
389  end if
390 
391  end do
392  end do
393  end do
394  !
395  ! Viscous terms in i-direction.
396  !
397  do k = 2, kl
398  do j = 2, jl
399  do i = 2, il
400 
401  ! Compute the metrics in xi-direction, i.e. along the
402  ! line i = constant.
403 
404  voli = one / vol(i, j, k)
405  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
406  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
407 
408  xm = si(i - 1, j, k, 1) * volmi
409  ym = si(i - 1, j, k, 2) * volmi
410  zm = si(i - 1, j, k, 3) * volmi
411  xp = si(i, j, k, 1) * volpi
412  yp = si(i, j, k, 2) * volpi
413  zp = si(i, j, k, 3) * volpi
414 
415  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
416  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
417  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
418  ttm = xm * xa + ym * ya + zm * za
419  ttp = xp * xa + yp * ya + zp * za
420 
421  ! Computation of the viscous terms in xi-direction; note
422  ! that cross-derivatives are neglected, i.e. the mesh is
423  ! assumed to be orthogonal.
424  ! The second derivative in xi-direction is constructed as
425  ! the central difference of the first order derivatives, i.e.
426  ! d^2/dxi^2 = d/dxi (d/dxi i+1/2 - d/dxi i-1/2).
427  ! In this way the metric as well as the varying viscosity
428  ! can be taken into account; the latter appears inside the
429  ! d/dxi derivative. The whole term is divided by rho to
430  ! obtain the diffusion term for k and omega.
431 
432  rhoi = one / w(i, j, k, irho)
433  mulm = half * (rlv(i - 1, j, k) + rlv(i, j, k))
434  mulp = half * (rlv(i + 1, j, k) + rlv(i, j, k))
435  muem = half * (rev(i - 1, j, k) + rev(i, j, k))
436  muep = half * (rev(i + 1, j, k) + rev(i, j, k))
437 
438  c1m = ttm * (mulm + sig1 * muem) * rhoi
439  c1p = ttp * (mulp + sig1 * muep) * rhoi
440  c10 = c1m + c1p
441 
442  c2m = ttm * (mulm + sig2 * muem) * rhoi
443  c2p = ttp * (mulp + sig2 * muep) * rhoi
444  c20 = c2m + c2p
445 
446  ! Update the residual for this cell and store the possible
447  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
448 
449  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i - 1, j, k, itu1) &
450  - c10 * w(i, j, k, itu1) + c1p * w(i + 1, j, k, itu1)
451  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i - 1, j, k, itu2) &
452  - c20 * w(i, j, k, itu2) + c2p * w(i + 1, j, k, itu2)
453 
454  b1 = -c1m
455  c1 = c10
456  d1 = -c1p
457 
458  b2 = -c2m
459  c2 = c20
460  d2 = -c2p
461 
462  ! Update the central jacobian. For nonboundary cells this
463  ! is simply c1 and c2. For boundary cells this is slightly
464  ! more complicated, because the boundary conditions are
465  ! treated implicitly and the off-diagonal terms b1, b2 and
466  ! d1, d2 must be taken into account.
467  ! The boundary conditions are only treated implicitly if
468  ! the diagonal dominance of the matrix is increased.
469 
470  if (i == 2) then
471  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
472  - b1 * max(bmti1(j, k, itu1, itu1), zero)
473  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmti1(j, k, itu1, itu2)
474  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmti1(j, k, itu2, itu1)
475  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
476  - b2 * max(bmti1(j, k, itu2, itu2), zero)
477  else if (i == il) then
478  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
479  - d1 * max(bmti2(j, k, itu1, itu1), zero)
480  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmti2(j, k, itu1, itu2)
481  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmti2(j, k, itu2, itu1)
482  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
483  - d2 * max(bmti2(j, k, itu2, itu2), zero)
484  else
485  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
486  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
487  end if
488 
489  end do
490  end do
491  end do
492 
493  ! Multiply the residual by the volume and store this in dw; this
494  ! is done for monitoring reasons only. The multiplication with the
495  ! volume is present to be consistent with the flow residuals; also
496  ! the negative value is taken, again to be consistent with the
497  ! flow equations. Also multiply by iblank so that no updates occur
498  ! in holes or the overset boundary.
499 
500  do k = 2, kl
501  do j = 2, jl
502  do i = 2, il
503  rblank = real(iblank(i, j, k), realtype)
504  dw(i, j, k, itu1) = -volref(i, j, k) * dvt(i, j, k, 1) * rblank
505  dw(i, j, k, itu2) = -volref(i, j, k) * dvt(i, j, k, 2) * rblank
506  end do
507  end do
508  end do
509 
510  ! Initialize the wall function flags to .false.
511 
512  flagi2 = .false.
513  flagil = .false.
514  flagj2 = .false.
515  flagjl = .false.
516  flagk2 = .false.
517  flagkl = .false.
518 
519  ! Modify the rhs of the 1st internal cell, if wall functions
520  ! are used; their value is determined by the table.
521 
522  testwallfunctions: if (wallfunctions) then
523 
524  bocos: do nn = 1, nviscbocos
525 
526  ! Determine the block face on which the subface is located
527  ! and set some variables. As flag points to the entire array
528  ! flagI2, etc., its starting indices are the starting indices
529  ! of its target and not 1.
530 
531  select case (bcfaceid(nn))
532  case (imin)
533  flag => flagi2
534  ddw => dw(2, 1:, 1:, 1:); ddvt => dvt(2, 1:, 1:, 1:)
535  ww => w(2, 1:, 1:, 1:); rrlv => rlv(2, 1:, 1:)
536  dd2wall => d2wall(2, :, :)
537 
538  case (imax)
539  flag => flagil
540  ddw => dw(il, 1:, 1:, 1:); ddvt => dvt(il, 1:, 1:, 1:)
541  ww => w(il, 1:, 1:, 1:); rrlv => rlv(il, 1:, 1:)
542  dd2wall => d2wall(il, :, :)
543 
544  case (jmin)
545  flag => flagj2
546  ddw => dw(1:, 2, 1:, 1:); ddvt => dvt(1:, 2, 1:, 1:)
547  ww => w(1:, 2, 1:, 1:); rrlv => rlv(1:, 2, 1:)
548  dd2wall => d2wall(:, 2, :)
549 
550  case (jmax)
551  flag => flagjl
552  ddw => dw(1:, jl, 1:, 1:); ddvt => dvt(1:, jl, 1:, 1:)
553  ww => w(1:, jl, 1:, 1:); rrlv => rlv(1:, jl, 1:)
554  dd2wall => d2wall(:, jl, :)
555 
556  case (kmin)
557  flag => flagk2
558  ddw => dw(1:, 1:, 2, 1:); ddvt => dvt(1:, 1:, 2, 1:)
559  ww => w(1:, 1:, 2, 1:); rrlv => rlv(1:, 1:, 2)
560  dd2wall => d2wall(:, :, 2)
561 
562  case (kmax)
563  flag => flagkl
564  ddw => dw(1:, 1:, kl, 1:); ddvt => dvt(1:, 1:, kl, 1:)
565  ww => w(1:, 1:, kl, 1:); rrlv => rlv(1:, 1:, kl)
566  dd2wall => d2wall(:, :, kl)
567 
568  end select
569 
570  ! Loop over the owned faces of this subface. Therefore the
571  ! nodal range of BCData must be used. The offset of +1 is
572  ! present, because the starting index of the cell range is
573  ! 1 larger than the starting index of the nodal range.
574 
575  do j = (bcdata(nn)%jnBeg + 1), bcdata(nn)%jnEnd
576  do i = (bcdata(nn)%inBeg + 1), bcdata(nn)%inEnd
577 
578  ! Set ddw to zero for proper monitoring of the
579  ! convergence.
580 
581  ddw(i, j, itu1) = zero
582  ddw(i, j, itu2) = zero
583 
584  ! Enforce k and omega in the 1st internal cell from
585  ! the wall function table. There is an offset of -1 in
586  ! the wall distance. Note that the offset compared to
587  ! the current value must be stored. Also note that the
588  ! curve fits contain the non-dimensional values.
589 
590  utau = viscsubface(nn)%utau(i, j)
591  yp = ww(i, j, irho) * dd2wall(i - 1, j - 1) * utau / rrlv(i, j)
592 
593  call curvetupyp(tup, yp, itu1, itu2)
594 
595  tup(itu1) = tup(itu1) * utau**2
596  tup(itu2) = tup(itu2) * utau**2 / rrlv(i, j) * ww(i, j, irho)
597 
598  ddvt(i, j, 1) = tup(itu1) - ww(i, j, itu1)
599  ddvt(i, j, 2) = tup(itu2) - ww(i, j, itu2)
600 
601  ! Set the wall flag to .true.
602 
603  flag(i, j) = .true.
604 
605  end do
606  end do
607 
608  end do bocos
609  end if testwallfunctions
610 
611  ! Return if only the residual must be computed.
612 
613  if (resonly) return
614 
615  ! For implicit relaxation take the local time step into account,
616  ! where dt is the inverse of the central jacobian times the cfl
617  ! number. The following system is solved:
618  ! (i/dt + cc + bb + dd)*dw = rhs, in which i/dt = cc/cfl. As in
619  ! the rest of the algorithm only the modified central jacobian is
620  ! used, stored it now.
621 
622  ! Compute the factor multiplying the central jacobian, which
623  ! is 1 + 1/cfl (implicit relaxation only).
624 
625  factor = one
626  if (turbrelax == turbrelaximplicit) &
627  factor = one + (one - alfaturb) / alfaturb
628 
629  do k = 2, kl
630  do j = 2, jl
631  do i = 2, il
632  qq(i, j, k, 1, 1) = factor * qq(i, j, k, 1, 1)
633  qq(i, j, k, 1, 2) = factor * qq(i, j, k, 1, 2)
634  qq(i, j, k, 2, 1) = factor * qq(i, j, k, 2, 1)
635  qq(i, j, k, 2, 2) = factor * qq(i, j, k, 2, 2)
636 
637  ! Set qq to 1 if the value is determined by the table.
638 
639  if ((i == 2 .and. flagi2(j, k)) .or. &
640  (i == il .and. flagil(j, k)) .or. &
641  (j == 2 .and. flagj2(i, k)) .or. &
642  (j == jl .and. flagjl(i, k)) .or. &
643  (k == 2 .and. flagk2(i, j)) .or. &
644  (k == kl .and. flagkl(i, j))) then
645  qq(i, j, k, 1, 1) = one
646  qq(i, j, k, 1, 2) = zero
647  qq(i, j, k, 2, 1) = zero
648  qq(i, j, k, 2, 2) = one
649  end if
650  end do
651  end do
652  end do
653 
654  ! Initialize the grid velocity to zero. This value will be used
655  ! if the block is not moving.
656 
657  qs = zero
658  !
659  ! dd-ADI step in j-direction. There is no particular reason to
660  ! start in j-direction, it just happened to be so. As we solve
661  ! in j-direction, the j-loop is the innermost loop.
662  !
663  do k = 2, kl
664  do i = 2, il
665  do j = 2, jl
666 
667  ! More or less the same code is executed here as above when
668  ! the residual was built. However, now the off-diagonal
669  ! terms for the dd-ADI must be built and stored. This could
670  ! have been done earlier, but then all the coefficients had
671  ! to be stored. To save memory, they are recomputed.
672  ! Consequently, see the j-loop to build the residual for
673  ! the comments.
674 
675  voli = one / vol(i, j, k)
676  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
677  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
678 
679  xm = sj(i, j - 1, k, 1) * volmi
680  ym = sj(i, j - 1, k, 2) * volmi
681  zm = sj(i, j - 1, k, 3) * volmi
682  xp = sj(i, j, k, 1) * volpi
683  yp = sj(i, j, k, 2) * volpi
684  zp = sj(i, j, k, 3) * volpi
685 
686  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
687  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
688  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
689  ttm = xm * xa + ym * ya + zm * za
690  ttp = xp * xa + yp * ya + zp * za
691 
692  ! Off-diagonal terms due to the diffusion terms
693  ! in j-direction.
694 
695  rhoi = one / w(i, j, k, irho)
696  mulm = half * (rlv(i, j - 1, k) + rlv(i, j, k))
697  mulp = half * (rlv(i, j + 1, k) + rlv(i, j, k))
698  muem = half * (rev(i, j - 1, k) + rev(i, j, k))
699  muep = half * (rev(i, j + 1, k) + rev(i, j, k))
700 
701  c1m = ttm * (mulm + sig1 * muem) * rhoi
702  c1p = ttp * (mulp + sig1 * muep) * rhoi
703 
704  c2m = ttm * (mulm + sig2 * muem) * rhoi
705  c2p = ttp * (mulp + sig2 * muep) * rhoi
706 
707  bb(1, j) = -c1m
708  dd(1, j) = -c1p
709  bb(2, j) = -c2m
710  dd(2, j) = -c2p
711 
712  ! Compute the grid velocity if present.
713  ! It is taken as the average of j and j-1,
714 
715  if (addgridvelocities) &
716  qs = half * (sfacej(i, j, k) + sfacej(i, j - 1, k)) * voli
717 
718  ! Off-diagonal terms due to the advection term in
719  ! j-direction. First order approximation.
720 
721  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
722  um = zero
723  up = zero
724  if (uu < zero) um = uu
725  if (uu > zero) up = uu
726 
727  bb(1, j) = bb(1, j) - up
728  dd(1, j) = dd(1, j) + um
729  bb(2, j) = bb(2, j) - up
730  dd(2, j) = dd(2, j) + um
731 
732  ! Store the central jacobian and rhs in cc and ff.
733  ! Multiply the off-diagonal terms and rhs by the iblank
734  ! value so the update determined for iblank = 0 is zero.
735 
736  rblank = real(iblank(i, j, k), realtype)
737 
738  cc(1, 1, j) = qq(i, j, k, 1, 1)
739  cc(1, 2, j) = qq(i, j, k, 1, 2) * rblank
740  cc(2, 1, j) = qq(i, j, k, 2, 1) * rblank
741  cc(2, 2, j) = qq(i, j, k, 2, 2)
742 
743  ff(1, j) = dvt(i, j, k, 1) * rblank
744  ff(2, j) = dvt(i, j, k, 2) * rblank
745 
746  bb(:, j) = bb(:, j) * rblank
747  dd(:, j) = dd(:, j) * rblank
748 
749  ! Set off diagonal terms to zero if wall function are used.
750 
751  if ((i == 2 .and. flagi2(j, k)) .or. &
752  (i == il .and. flagil(j, k)) .or. &
753  (j == 2 .and. flagj2(i, k)) .or. &
754  (j == jl .and. flagjl(i, k)) .or. &
755  (k == 2 .and. flagk2(i, j)) .or. &
756  (k == kl .and. flagkl(i, j))) then
757  bb(1, j) = zero
758  dd(1, j) = zero
759  bb(2, j) = zero
760  dd(2, j) = zero
761  end if
762 
763  end do
764 
765  ! Solve the tri-diagonal system in j-direction.
766 
767  call tdia3(2_inttype, jl, bb, cc, dd, ff)
768 
769  ! Determine the new rhs for the next direction.
770 
771  do j = 2, jl
772  dvt(i, j, k, 1) = qq(i, j, k, 1, 1) * ff(1, j) + qq(i, j, k, 1, 2) * ff(2, j)
773  dvt(i, j, k, 2) = qq(i, j, k, 2, 1) * ff(1, j) + qq(i, j, k, 2, 2) * ff(2, j)
774  end do
775 
776  end do
777  end do
778  !
779  ! dd-ADI step in i-direction. As we solve in i-direction, the
780  ! i-loop is the innermost loop.
781  !
782  do k = 2, kl
783  do j = 2, jl
784  do i = 2, il
785 
786  ! More or less the same code is executed here as above when
787  ! the residual was built. However, now the off-diagonal
788  ! terms for the dd-ADI must be built and stored. This could
789  ! have been done earlier, but then all the coefficients had
790  ! to be stored. To save memory, they are recomputed.
791  ! Consequently, see the i-loop to build the residual for
792  ! the comments.
793 
794  voli = one / vol(i, j, k)
795  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
796  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
797 
798  xm = si(i - 1, j, k, 1) * volmi
799  ym = si(i - 1, j, k, 2) * volmi
800  zm = si(i - 1, j, k, 3) * volmi
801  xp = si(i, j, k, 1) * volpi
802  yp = si(i, j, k, 2) * volpi
803  zp = si(i, j, k, 3) * volpi
804 
805  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
806  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
807  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
808  ttm = xm * xa + ym * ya + zm * za
809  ttp = xp * xa + yp * ya + zp * za
810 
811  ! Off-diagonal terms due to the diffusion terms
812  ! in i-direction.
813 
814  rhoi = one / w(i, j, k, irho)
815  mulm = half * (rlv(i - 1, j, k) + rlv(i, j, k))
816  mulp = half * (rlv(i + 1, j, k) + rlv(i, j, k))
817  muem = half * (rev(i - 1, j, k) + rev(i, j, k))
818  muep = half * (rev(i + 1, j, k) + rev(i, j, k))
819 
820  c1m = ttm * (mulm + sig1 * muem) * rhoi
821  c1p = ttp * (mulp + sig1 * muep) * rhoi
822 
823  c2m = ttm * (mulm + sig2 * muem) * rhoi
824  c2p = ttp * (mulp + sig2 * muep) * rhoi
825 
826  bb(1, i) = -c1m
827  dd(1, i) = -c1p
828  bb(2, i) = -c2m
829  dd(2, i) = -c2p
830 
831  ! Compute the grid velocity if present.
832  ! It is taken as the average of i and i-1,
833 
834  if (addgridvelocities) &
835  qs = half * (sfacei(i, j, k) + sfacei(i - 1, j, k)) * voli
836 
837  ! Off-diagonal terms due to the advection term in
838  ! i-direction. First order approximation.
839 
840  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
841  um = zero
842  up = zero
843  if (uu < zero) um = uu
844  if (uu > zero) up = uu
845 
846  bb(1, i) = bb(1, i) - up
847  dd(1, i) = dd(1, i) + um
848  bb(2, i) = bb(2, i) - up
849  dd(2, i) = dd(2, i) + um
850 
851  ! Store the central jacobian and rhs in cc and ff.
852  ! Multiply the off-diagonal terms and rhs by the iblank
853  ! value so the update determined for iblank = 0 is zero.
854 
855  rblank = real(iblank(i, j, k), realtype)
856 
857  cc(1, 1, i) = qq(i, j, k, 1, 1)
858  cc(1, 2, i) = qq(i, j, k, 1, 2) * rblank
859  cc(2, 1, i) = qq(i, j, k, 2, 1) * rblank
860  cc(2, 2, i) = qq(i, j, k, 2, 2)
861 
862  ff(1, i) = dvt(i, j, k, 1) * rblank
863  ff(2, i) = dvt(i, j, k, 2) * rblank
864 
865  bb(:, i) = bb(:, i) * rblank
866  dd(:, i) = dd(:, i) * rblank
867 
868  ! Set off diagonal terms to zero if wall function are used.
869 
870  if ((i == 2 .and. flagi2(j, k)) .or. &
871  (i == il .and. flagil(j, k)) .or. &
872  (j == 2 .and. flagj2(i, k)) .or. &
873  (j == jl .and. flagjl(i, k)) .or. &
874  (k == 2 .and. flagk2(i, j)) .or. &
875  (k == kl .and. flagkl(i, j))) then
876  bb(1, i) = zero
877  dd(1, i) = zero
878  bb(2, i) = zero
879  dd(2, i) = zero
880  end if
881 
882  end do
883 
884  ! Solve the tri-diagonal system in i-direction.
885 
886  call tdia3(2_inttype, il, bb, cc, dd, ff)
887 
888  ! Determine the new rhs for the next direction.
889 
890  do i = 2, il
891  dvt(i, j, k, 1) = qq(i, j, k, 1, 1) * ff(1, i) + qq(i, j, k, 1, 2) * ff(2, i)
892  dvt(i, j, k, 2) = qq(i, j, k, 2, 1) * ff(1, i) + qq(i, j, k, 2, 2) * ff(2, i)
893  end do
894 
895  end do
896  end do
897  !
898  ! dd-adi step in k-direction. As we solve in k-direction, the
899  ! k-loop is the innermost loop.
900  !
901  do j = 2, jl
902  do i = 2, il
903  do k = 2, kl
904 
905  ! More or less the same code is executed here as above when
906  ! the residual was built. However, now the off-diagonal
907  ! terms for the dd-ADI must be built and stored. This could
908  ! have been done earlier, but then all the coefficients had
909  ! to be stored. To save memory, they are recomputed.
910  ! Consequently, see the k-loop to build the residual for
911  ! the comments.
912 
913  voli = one / vol(i, j, k)
914  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
915  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
916 
917  xm = sk(i, j, k - 1, 1) * volmi
918  ym = sk(i, j, k - 1, 2) * volmi
919  zm = sk(i, j, k - 1, 3) * volmi
920  xp = sk(i, j, k, 1) * volpi
921  yp = sk(i, j, k, 2) * volpi
922  zp = sk(i, j, k, 3) * volpi
923 
924  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
925  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
926  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
927  ttm = xm * xa + ym * ya + zm * za
928  ttp = xp * xa + yp * ya + zp * za
929 
930  ! Off-diagonal terms due to the diffusion terms
931  ! in k-direction.
932 
933  rhoi = one / w(i, j, k, irho)
934  mulm = half * (rlv(i, j, k - 1) + rlv(i, j, k))
935  mulp = half * (rlv(i, j, k + 1) + rlv(i, j, k))
936  muem = half * (rev(i, j, k - 1) + rev(i, j, k))
937  muep = half * (rev(i, j, k + 1) + rev(i, j, k))
938 
939  c1m = ttm * (mulm + sig1 * muem) * rhoi
940  c1p = ttp * (mulp + sig1 * muep) * rhoi
941 
942  c2m = ttm * (mulm + sig2 * muem) * rhoi
943  c2p = ttp * (mulp + sig2 * muep) * rhoi
944 
945  bb(1, k) = -c1m
946  dd(1, k) = -c1p
947  bb(2, k) = -c2m
948  dd(2, k) = -c2p
949 
950  ! Compute the grid velocity if present.
951  ! It is taken as the average of k and k-1,
952 
953  if (addgridvelocities) &
954  qs = half * (sfacek(i, j, k) + sfacek(i, j, k - 1)) * voli
955 
956  ! Off-diagonal terms due to the advection term in
957  ! k-direction. First order approximation.
958 
959  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
960  um = zero
961  up = zero
962  if (uu < zero) um = uu
963  if (uu > zero) up = uu
964 
965  bb(1, k) = bb(1, k) - up
966  dd(1, k) = dd(1, k) + um
967  bb(2, k) = bb(2, k) - up
968  dd(2, k) = dd(2, k) + um
969 
970  ! Store the central jacobian and rhs in cc and ff.
971  ! Multiply the off-diagonal terms and rhs by the iblank
972  ! value so the update determined for iblank = 0 is zero.
973 
974  rblank = real(iblank(i, j, k), realtype)
975 
976  cc(1, 1, k) = qq(i, j, k, 1, 1)
977  cc(1, 2, k) = qq(i, j, k, 1, 2) * rblank
978  cc(2, 1, k) = qq(i, j, k, 2, 1) * rblank
979  cc(2, 2, k) = qq(i, j, k, 2, 2)
980 
981  ff(1, k) = dvt(i, j, k, 1) * rblank
982  ff(2, k) = dvt(i, j, k, 2) * rblank
983 
984  bb(:, k) = bb(:, k) * rblank
985  dd(:, k) = dd(:, k) * rblank
986 
987  ! Set off diagonal terms to zero if wall function are used.
988 
989  if ((i == 2 .and. flagi2(j, k)) .or. &
990  (i == il .and. flagil(j, k)) .or. &
991  (j == 2 .and. flagj2(i, k)) .or. &
992  (j == jl .and. flagjl(i, k)) .or. &
993  (k == 2 .and. flagk2(i, j)) .or. &
994  (k == kl .and. flagkl(i, j))) then
995  bb(1, k) = zero
996  dd(1, k) = zero
997  bb(2, k) = zero
998  dd(2, k) = zero
999  end if
1000 
1001  end do
1002 
1003  ! Solve the tri-diagonal system in k-direction.
1004 
1005  call tdia3(2_inttype, kl, bb, cc, dd, ff)
1006 
1007  ! Store the update in dvt.
1008 
1009  do k = 2, kl
1010  dvt(i, j, k, 1) = ff(1, k)
1011  dvt(i, j, k, 2) = ff(2, k)
1012  end do
1013 
1014  end do
1015  end do
1016  !
1017  ! Update the turbulent variables. For explicit relaxation the
1018  ! update must be relaxed; for implicit relaxation this has been
1019  ! done via the time step.
1020  !
1021  factor = one
1022  if (turbrelax == turbrelaxexplicit) factor = alfaturb
1023 
1024  do k = 2, kl
1025  do j = 2, jl
1026  do i = 2, il
1027  w(i, j, k, itu1) = w(i, j, k, itu1) + factor * dvt(i, j, k, 1)
1028  w(i, j, k, itu1) = max(w(i, j, k, itu1), zero)
1029 
1030  w(i, j, k, itu2) = w(i, j, k, itu2) + factor * dvt(i, j, k, 2)
1031  w(i, j, k, itu2) = max(w(i, j, k, itu2), 1.e-5_realtype * winf(itu2))
1032  end do
1033  end do
1034  end do
1035 
1036  end subroutine kwsolve
1037 end module kw
bcdata
Definition: BCData.F90:1
blockpointers::bmtk2
real(kind=realtype), dimension(:, :, :, :), pointer bmtk2
Definition: blockPointers.F90:162
blockpointers::sfacek
real(kind=realtype), dimension(:, :, :), pointer sfacek
Definition: blockPointers.F90:119
blockpointers::addgridvelocities
logical addgridvelocities
Definition: blockPointers.F90:117
blockpointers::jl
integer(kind=inttype) jl
Definition: blockPointers.F90:38
blockpointers::bmti1
real(kind=realtype), dimension(:, :, :, :), pointer bmti1
Definition: blockPointers.F90:157
blockpointers::nviscbocos
integer(kind=inttype) nviscbocos
Definition: blockPointers.F90:47
blockpointers::bmtj1
real(kind=realtype), dimension(:, :, :, :), pointer bmtj1
Definition: blockPointers.F90:159
blockpointers::bmti2
real(kind=realtype), dimension(:, :, :, :), pointer bmti2
Definition: blockPointers.F90:158
blockpointers::w
real(kind=realtype), dimension(:, :, :, :), pointer w
Definition: blockPointers.F90:120
blockpointers::scratch
real(kind=realtype), dimension(:, :, :, :), pointer scratch
Definition: blockPointers.F90:129
blockpointers::sfacei
real(kind=realtype), dimension(:, :, :), pointer sfacei
Definition: blockPointers.F90:119
blockpointers::viscsubface
type(viscsubfacetype), dimension(:), pointer viscsubface
Definition: blockPointers.F90:85
blockpointers::d2wall
real(kind=realtype), dimension(:, :, :), pointer d2wall
Definition: blockPointers.F90:154
blockpointers::iblank
integer(kind=inttype), dimension(:, :, :), pointer iblank
Definition: blockPointers.F90:71
blockpointers::rlv
real(kind=realtype), dimension(:, :, :), pointer rlv
Definition: blockPointers.F90:125
blockpointers::bcfaceid
integer(kind=inttype), dimension(:), pointer bcfaceid
Definition: blockPointers.F90:50
blockpointers::si
real(kind=realtype), dimension(:, :, :, :), pointer si
Definition: blockPointers.F90:96
blockpointers::volref
real(kind=realtype), dimension(:, :, :), pointer volref
Definition: blockPointers.F90:98
blockpointers::sj
real(kind=realtype), dimension(:, :, :, :), pointer sj
Definition: blockPointers.F90:96
blockpointers::rev
real(kind=realtype), dimension(:, :, :), pointer rev
Definition: blockPointers.F90:125
blockpointers::bmtj2
real(kind=realtype), dimension(:, :, :, :), pointer bmtj2
Definition: blockPointers.F90:160
blockpointers::dw
real(kind=realtype), dimension(:, :, :, :), pointer dw
Definition: blockPointers.F90:128
blockpointers::sk
real(kind=realtype), dimension(:, :, :, :), pointer sk
Definition: blockPointers.F90:96
blockpointers::vol
real(kind=realtype), dimension(:, :, :), pointer vol
Definition: blockPointers.F90:97
blockpointers::sfacej
real(kind=realtype), dimension(:, :, :), pointer sfacej
Definition: blockPointers.F90:119
blockpointers::bmtk1
real(kind=realtype), dimension(:, :, :, :), pointer bmtk1
Definition: blockPointers.F90:161
blockpointers::kl
integer(kind=inttype) kl
Definition: blockPointers.F90:38
blockpointers::il
integer(kind=inttype) il
Definition: blockPointers.F90:38
constants::strain
integer(kind=inttype), parameter strain
Definition: constants.F90:137
constants::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::imax
integer(kind=inttype), parameter imax
Definition: constants.F90:293
constants::idvt
integer, parameter idvt
Definition: constants.F90:51
constants::kmin
integer(kind=inttype), parameter kmin
Definition: constants.F90:296
constants::jmax
integer(kind=inttype), parameter jmax
Definition: constants.F90:295
constants::itu2
integer, parameter itu2
Definition: constants.F90:42
constants::icd
integer, parameter icd
Definition: constants.F90:56
constants::itu1
integer, parameter itu1
Definition: constants.F90:40
constants::irho
integer, parameter irho
Definition: constants.F90:34
constants::vorticity
integer(kind=inttype), parameter vorticity
Definition: constants.F90:137
constants::komegamodified
integer(kind=inttype), parameter komegamodified
Definition: constants.F90:128
constants::ivx
integer, parameter ivx
Definition: constants.F90:35
constants::iprod
integer, parameter iprod
Definition: constants.F90:55
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
constants::turbrelaximplicit
integer(kind=inttype), parameter turbrelaximplicit
Definition: constants.F90:223
constants::turbrelaxexplicit
integer(kind=inttype), parameter turbrelaxexplicit
Definition: constants.F90:223
constants::katolaunder
integer(kind=inttype), parameter katolaunder
Definition: constants.F90:137
constants::imin
integer(kind=inttype), parameter imin
Definition: constants.F90:292
constants::ivz
integer, parameter ivz
Definition: constants.F90:37
constants::two
real(kind=realtype), parameter two
Definition: constants.F90:73
constants::kmax
integer(kind=inttype), parameter kmax
Definition: constants.F90:297
constants::ivy
integer, parameter ivy
Definition: constants.F90:36
constants::jmin
integer(kind=inttype), parameter jmin
Definition: constants.F90:294
flowvarrefstate::winf
real(kind=realtype), dimension(:), allocatable winf
Definition: flowVarRefState.F90:74
inputiteration::alfaturb
real(kind=realtype) alfaturb
Definition: inputParam.F90:276
inputiteration::turbrelax
integer(kind=inttype) turbrelax
Definition: inputParam.F90:269
inputphysics::turbprod
integer(kind=inttype) turbprod
Definition: inputParam.F90:584
inputphysics::turbmodel
integer(kind=inttype) turbmodel
Definition: inputParam.F90:584
inputphysics::wallfunctions
logical wallfunctions
Definition: inputParam.F90:589
inputphysics::pklim
real(kind=realtype) pklim
Definition: inputParam.F90:596
kw
Definition: kw.F90:1
kw::kw_block
subroutine kw_block(resOnly)
Definition: kw.F90:6
kw::kwsolve
subroutine kwsolve(resOnly)
Definition: kw.F90:52
paramturb::rkwsigd1
real(kind=realtype), parameter rkwsigd1
Definition: paramTurb.F90:34
paramturb::rkwbetas
real(kind=realtype), parameter rkwbetas
Definition: paramTurb.F90:36
paramturb::rkwsigk1
real(kind=realtype), parameter rkwsigk1
Definition: paramTurb.F90:32
paramturb::rkwk
real(kind=realtype), parameter rkwk
Definition: paramTurb.F90:31
paramturb::rkwsigw1
real(kind=realtype), parameter rkwsigw1
Definition: paramTurb.F90:33
paramturb::rkwbeta1
real(kind=realtype), parameter rkwbeta1
Definition: paramTurb.F90:35
turbbcroutines::bcturbtreatment
subroutine bcturbtreatment
Definition: turbBCRoutines.F90:663
turbbcroutines::applyallturbbcthisblock
subroutine applyallturbbcthisblock(secondHalo)
Definition: turbBCRoutines.F90:50
turbcurvefits::curvetupyp
subroutine curvetupyp(tup, yp, ntu1, ntu2)
Definition: turbCurveFits.F90:104
turbmod::vort
real(kind=realtype), dimension(:, :, :), pointer vort
Definition: turbMod.F90:34
turbmod::sig2
real(kind=realtype) sig2
Definition: turbMod.F90:16
turbmod::kwcd
real(kind=realtype), dimension(:, :, :), pointer kwcd
Definition: turbMod.F90:37
turbmod::prod
real(kind=realtype), dimension(:, :, :), pointer prod
Definition: turbMod.F90:35
turbmod::sig1
real(kind=realtype) sig1
Definition: turbMod.F90:16
turbmod::dvt
real(kind=realtype), dimension(:, :, :, :), pointer dvt
Definition: turbMod.F90:33
turbmod::sct
real(kind=realtype), dimension(:, :, :), pointer sct
Definition: turbMod.F90:39
turbutils::unsteadyturbterm
subroutine unsteadyturbterm(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:411
turbutils::kwcdterm
subroutine kwcdterm
Definition: turbUtils.F90:2069
turbutils::prodwmag2
subroutine prodwmag2
Definition: turbUtils.F90:240
turbutils::tdia3
subroutine tdia3(nb, ne, l, c, u, r)
Definition: turbUtils.F90:1561
turbutils::turbadvection
subroutine turbadvection(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:826
turbutils::prodsmag2
subroutine prodsmag2
Definition: turbUtils.F90:130
turbutils::prodkatolaunder
subroutine prodkatolaunder
Definition: turbUtils.F90:6
turbutils::kweddyviscosity
subroutine kweddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:713
utils
Definition: utils.F90:1
utils::setpointers
subroutine setpointers(nn, mm, ll)
Definition: utils.F90:3237
Generated by doxygen 1.9.1