ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
vf.F90
Go to the documentation of this file.
1 module vf
2 contains
3 
4  subroutine vf_block(resOnly)
5  !
6  ! vf solves the transport equations for the v2-f model
7  ! in a coupled manner using a diagonal dominant ADI-scheme.
8  !
9  use constants
10  use blockpointers, only: il, jl, kl
11  use inputtimespectral
12  use iteration
13  use utils, only: setpointers
14  use turbutils, only: vfeddyviscosity, vfscale, unsteadyturbterm
15  use turbbcroutines, only: bcturbtreatment, applyallturbbcthisblock
16 
17  implicit none
18  !
19  ! Subroutine argument.
20  !
21  logical, intent(in) :: resOnly
22 
23  ! Set the arrays for the boundary condition treatment.
24 
25  call bcturbtreatment
26 
27  ! Compute time and length scale
28 
29  call vfscale
30 
31  ! Solve the transport equations for k and epsilon.
32 
33  call kesolve(resonly)
34 
35  ! Solve the transport equation for v2 and the elliptic
36  ! equation for f.
37 
38  call vfsolve(resonly)
39 
40  ! The eddy viscosity and the boundary conditions are only
41  ! applied if actual updates have been computed in keSolve
42  ! and vfSolve.
43 
44  if (.not. resonly) then
45 
46  ! Compute the corresponding eddy viscosity.
47 
48  call vfeddyviscosity(2, il, 2, jl, 2, kl)
49 
50  ! Set the halo values for the turbulent variables.
51  ! We are on the finest mesh, so the second layer of halo
52  ! cells must be computed as well.
53 
54  call applyallturbbcthisblock(.true.)
55 
56  end if
57 
58  end subroutine vf_block
59 
60  subroutine vfsolve(resOnly)
61  !
62  ! vfSolve solves the v2 transport equation and the
63  ! f elliptic relaxation equation of the v2-f model
64  ! in a coupled manner using a diagonal dominant ADI-scheme.
65  !
66  use blockpointers
67  use constants
68  use flowvarrefstate
69  use inputiteration
70  use inputphysics
71  use paramturb
72  use turbmod, only: prod, dvt, sig1, sig2, sct, scl2
73  use turbutils, only: turbadvection, unsteadyturbterm, tdia3
74  use turbcurvefits, only: curvetupyp
75  implicit none
76  !
77  ! Subroutine arguments.
78  !
79  logical, intent(in) :: resOnly
80  !
81  ! Local variables.
82  !
83  integer(kind=intType) :: i, j, k, nn
84 
85  real(kind=realtype) :: rhoi, ss, spk, ff1, ff2, ff3, ss1, ss2, ss3
86  real(kind=realtype) :: voli, volmi, volpi
87  real(kind=realtype) :: xm, ym, zm, xp, yp, zp, xa, ya, za
88  real(kind=realtype) :: ttm, ttp, mulm, mulp, muem, muep
89  real(kind=realtype) :: c1m, c1p, c10, c2m, c2p, c20
90  real(kind=realtype) :: b1, b2, c1, c2, d1, d2
91  real(kind=realtype) :: qs, uu, um, up, utau
92  real(kind=realtype) :: tke, tep, tv2, tf2, tkea, tepa, tv2a
93  real(kind=realtype) :: tkel, tv2l, sle2i, stei
94  real(kind=realtype) :: rsct, rnu
95  real(kind=realtype) :: tu12, tu22, tu32, tu42, tu52
96  real(kind=realtype) :: rnu23, dtu23, rblank
97 
98  real(kind=realtype), dimension(itu1:itu5) :: tup
99 
100  real(kind=realtype), dimension(2:il, 2:jl, 2:kl, 2, 2) :: qq
101  real(kind=realtype), dimension(2, 2:max(il, jl, kl)) :: bb, dd, ff
102  real(kind=realtype), dimension(2, 2, 2:max(il, jl, kl)) :: cc
103 
104  real(kind=realtype), dimension(:, :, :), pointer :: dw2, dvt2, w2, w3
105  real(kind=realtype), dimension(:, :), pointer :: rlv2, rlv3
106  real(kind=realtype), dimension(:, :), pointer :: rev2, rev3
107  real(kind=realtype), dimension(:, :), pointer :: d2wall2, d2wall3
108 
109  logical, dimension(2:jl, 2:kl), target :: flagI2, flagIl
110  logical, dimension(2:il, 2:kl), target :: flagJ2, flagJl
111  logical, dimension(2:il, 2:jl), target :: flagK2, flagKl
112 
113  logical, dimension(:, :), pointer :: flag
114 
115  ! Set model constants
116 
117  sig1 = rvfsigv1
118  sig2 = one
119  !
120  ! Source terms.
121  ! Determine the source term and its derivative w.r.t. v2 and f
122  ! for all internal cells of the block.
123  !
124  do k = 2, kl
125  do j = 2, jl
126  do i = 2, il
127 
128  ! Compute the source terms for both the k and the epsilon
129  ! equation. Note that dw(i,j,k,iprod) contains the unscaled
130  ! production term.
131 
132  rhoi = one / w(i, j, k, irho)
133  tke = w(i, j, k, itu1)
134  tep = w(i, j, k, itu2)
135  tv2 = w(i, j, k, itu3)
136  tf2 = w(i, j, k, itu4)
137  tkea = abs(tke)
138  tepa = abs(tep)
139  tv2a = abs(tv2)
140  tkel = max(tkea, rvflimitk)
141  tv2l = max(tv2a, rvflimitk * 0.666666666_realtype)
142  stei = tepa / sct(i, j, k)
143  sle2i = tepa**2 / scl2(i, j, k)
144 
145  if (rvfn == 6) then
146  rnu = rlv(i, j, k) * rhoi
147  rsct = max(tkea, 6.*sqrt(rnu * tepa))
148  stei = tepa / rsct
149 
150  ! rn2 = rvfCn**2*(rnu*tepa)**1.5
151  ! sle2i= rvfCl**2*max(tkea**3,rn2)
152  end if
153 
154  ss = prod(i, j, k)
155  spk = rev(i, j, k) * ss * rhoi
156  spk = min(spk, pklim * tepa)
157 
158  ff1 = -tepa / tkel
159  ff2 = tkel
160  ff3 = zero
161 
162  ss1 = -(rvfc1 - 1.) / tkel * stei * sle2i
163  ss2 = -sle2i
164  ss3 = (rvfc1 - 1.) * 2./3.*stei * sle2i + rvfc2 * spk / tkel * sle2i
165 
166  if (rvfn == 6) then
167  ff1 = ff1 - 5.0 * tepa / tkel
168  ss1 = ss1 + 5.0 / tkel * stei * sle2i
169  end if
170 
171  dvt(i, j, k, 1) = ff1 * tv2 + ff2 * tf2 + ff3
172  dvt(i, j, k, 2) = ss1 * tv2 + ss2 * tf2 + ss3
173 
174  ! Compute the source term jacobian. Note that only the
175  ! destruction terms are linearized to increase the diagonal
176  ! dominance of the matrix. Furthermore minus the source
177  ! term jacobian is stored.
178 
179  qq(i, j, k, 1, 1) = -ff1
180  qq(i, j, k, 1, 2) = -ff2
181  qq(i, j, k, 2, 1) = -ss1
182  qq(i, j, k, 2, 2) = -ss2
183 
184  end do
185  end do
186  end do
187  !
188  ! Advection and unsteady terms.
189  !
190  nn = itu3 - 1
191  call turbadvection(2_inttype, 1_inttype, nn, qq)
192 
193  call unsteadyturbterm(2_inttype, 1_inttype, nn, qq)
194  !
195  ! Viscous terms in k-direction.
196  !
197  do k = 2, kl
198  do j = 2, jl
199  do i = 2, il
200 
201  ! Compute the metrics in zeta-direction, i.e. along the
202  ! line k = constant.
203 
204  voli = one / vol(i, j, k)
205  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
206  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
207 
208  xm = sk(i, j, k - 1, 1) * volmi
209  ym = sk(i, j, k - 1, 2) * volmi
210  zm = sk(i, j, k - 1, 3) * volmi
211  xp = sk(i, j, k, 1) * volpi
212  yp = sk(i, j, k, 2) * volpi
213  zp = sk(i, j, k, 3) * volpi
214 
215  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
216  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
217  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
218  ttm = xm * xa + ym * ya + zm * za
219  ttp = xp * xa + yp * ya + zp * za
220 
221  ! Computation of the viscous terms in zeta-direction; note
222  ! that cross-derivatives are neglected, i.e. the mesh is
223  ! assumed to be orthogonal.
224  ! The second derivative in zeta-direction is constructed as
225  ! the central difference of the first order derivatives, i.e.
226  ! d^2/dzeta^2 = d/dzeta (d/dzeta k+1/2 - d/dzeta k-1/2).
227  ! In this way the metric as well as the varying viscosity
228  ! can be taken into account; the latter appears inside the
229  ! d/dzeta derivative. The whole term is divided by rho to
230  ! obtain the diffusion term for v2 and f.
231 
232  rhoi = one / w(i, j, k, irho)
233  mulm = half * (rlv(i, j, k - 1) + rlv(i, j, k))
234  mulp = half * (rlv(i, j, k + 1) + rlv(i, j, k))
235  muem = half * (rev(i, j, k - 1) + rev(i, j, k))
236  muep = half * (rev(i, j, k + 1) + rev(i, j, k))
237 
238  c1m = ttm * (mulm + sig1 * muem) * rhoi
239  c1p = ttp * (mulp + sig1 * muep) * rhoi
240  c10 = c1m + c1p
241 
242  c2m = ttm
243  c2p = ttp
244  c20 = c2m + c2p
245 
246  ! Update the residual for this cell and store the possible
247  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
248 
249  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i, j, k - 1, itu3) &
250  - c10 * w(i, j, k, itu3) + c1p * w(i, j, k + 1, itu3)
251  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i, j, k - 1, itu4) &
252  - c20 * w(i, j, k, itu4) + c2p * w(i, j, k + 1, itu4)
253 
254  b1 = -c1m
255  c1 = c10
256  d1 = -c1p
257 
258  b2 = -c2m
259  c2 = c20
260  d2 = -c2p
261 
262  ! Update the central jacobian. For nonboundary cells this
263  ! is simply c1 and c2. For boundary cells this is slightly
264  ! more complicated, because the boundary conditions are
265  ! treated implicitly and the off-diagonal terms b1, b2 and
266  ! d1, d2 must be taken into account.
267  ! The boundary conditions are only treated implicitly if
268  ! the diagonal dominance of the matrix is increased.
269 
270  if (k == 2) then
271  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
272  - b1 * max(bmtk1(i, j, itu3, itu3), zero)
273  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmtk1(i, j, itu3, itu4)
274  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmtk1(i, j, itu4, itu3)
275  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
276  - b2 * max(bmtk1(i, j, itu4, itu4), zero)
277  else if (k == kl) then
278  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
279  - d1 * max(bmtk2(i, j, itu3, itu3), zero)
280  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmtk2(i, j, itu3, itu4)
281  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmtk2(i, j, itu4, itu3)
282  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
283  - d2 * max(bmtk2(i, j, itu4, itu4), zero)
284  else
285  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
286  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
287  end if
288 
289  end do
290  end do
291  end do
292  !
293  ! Viscous terms in j-direction.
294  !
295  do k = 2, kl
296  do j = 2, jl
297  do i = 2, il
298 
299  ! Compute the metrics in eta-direction, i.e. along the
300  ! line j = constant.
301 
302  voli = one / vol(i, j, k)
303  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
304  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
305 
306  xm = sj(i, j - 1, k, 1) * volmi
307  ym = sj(i, j - 1, k, 2) * volmi
308  zm = sj(i, j - 1, k, 3) * volmi
309  xp = sj(i, j, k, 1) * volpi
310  yp = sj(i, j, k, 2) * volpi
311  zp = sj(i, j, k, 3) * volpi
312 
313  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
314  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
315  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
316  ttm = xm * xa + ym * ya + zm * za
317  ttp = xp * xa + yp * ya + zp * za
318 
319  ! Computation of the viscous terms in eta-direction; note
320  ! that cross-derivatives are neglected, i.e. the mesh is
321  ! assumed to be orthogonal.
322  ! The second derivative in eta-direction is constructed as
323  ! the central difference of the first order derivatives, i.e.
324  ! d^2/deta^2 = d/deta (d/deta j+1/2 - d/deta j-1/2).
325  ! In this way the metric as well as the varying viscosity
326  ! can be taken into account; the latter appears inside the
327  ! d/deta derivative. The whole term is divided by rho to
328  ! obtain the diffusion term for v2 and f.
329 
330  rhoi = one / w(i, j, k, irho)
331  mulm = half * (rlv(i, j - 1, k) + rlv(i, j, k))
332  mulp = half * (rlv(i, j + 1, k) + rlv(i, j, k))
333  muem = half * (rev(i, j - 1, k) + rev(i, j, k))
334  muep = half * (rev(i, j + 1, k) + rev(i, j, k))
335 
336  c1m = ttm * (mulm + sig1 * muem) * rhoi
337  c1p = ttp * (mulp + sig1 * muep) * rhoi
338  c10 = c1m + c1p
339 
340  c2m = ttm
341  c2p = ttp
342  c20 = c2m + c2p
343 
344  ! Update the residual for this cell and store the possible
345  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
346 
347  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i, j - 1, k, itu3) &
348  - c10 * w(i, j, k, itu3) + c1p * w(i, j + 1, k, itu3)
349  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i, j - 1, k, itu4) &
350  - c20 * w(i, j, k, itu4) + c2p * w(i, j + 1, k, itu4)
351 
352  b1 = -c1m
353  c1 = c10
354  d1 = -c1p
355 
356  b2 = -c2m
357  c2 = c20
358  d2 = -c2p
359 
360  ! Update the central jacobian. For nonboundary cells this
361  ! is simply c1 and c2. For boundary cells this is slightly
362  ! more complicated, because the boundary conditions are
363  ! treated implicitly and the off-diagonal terms b1, b2 and
364  ! d1, d2 must be taken into account.
365  ! The boundary conditions are only treated implicitly if
366  ! the diagonal dominance of the matrix is increased.
367 
368  if (j == 2) then
369  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
370  - b1 * max(bmtj1(i, k, itu3, itu3), zero)
371  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmtj1(i, k, itu3, itu4)
372  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmtj1(i, k, itu4, itu3)
373  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
374  - b2 * max(bmtj1(i, k, itu4, itu4), zero)
375  else if (j == jl) then
376  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
377  - d1 * max(bmtj2(i, k, itu3, itu3), zero)
378  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmtj2(i, k, itu3, itu4)
379  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmtj2(i, k, itu4, itu3)
380  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
381  - d2 * max(bmtj2(i, k, itu4, itu4), zero)
382  else
383  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
384  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
385  end if
386 
387  end do
388  end do
389  end do
390  !
391  ! Viscous terms in i-direction.
392  !
393  do k = 2, kl
394  do j = 2, jl
395  do i = 2, il
396 
397  ! Compute the metrics in xi-direction, i.e. along the
398  ! line i = constant.
399 
400  voli = one / vol(i, j, k)
401  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
402  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
403 
404  xm = si(i - 1, j, k, 1) * volmi
405  ym = si(i - 1, j, k, 2) * volmi
406  zm = si(i - 1, j, k, 3) * volmi
407  xp = si(i, j, k, 1) * volpi
408  yp = si(i, j, k, 2) * volpi
409  zp = si(i, j, k, 3) * volpi
410 
411  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
412  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
413  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
414  ttm = xm * xa + ym * ya + zm * za
415  ttp = xp * xa + yp * ya + zp * za
416 
417  ! Computation of the viscous terms in xi-direction; note
418  ! that cross-derivatives are neglected, i.e. the mesh is
419  ! assumed to be orthogonal.
420  ! The second derivative in xi-direction is constructed as
421  ! the central difference of the first order derivatives, i.e.
422  ! d^2/dxi^2 = d/dxi (d/dxi i+1/2 - d/dxi i-1/2).
423  ! In this way the metric as well as the varying viscosity
424  ! can be taken into account; the latter appears inside the
425  ! d/dxi derivative. The whole term is divided by rho to
426  ! obtain the diffusion term for v2 and f.
427 
428  rhoi = one / w(i, j, k, irho)
429  mulm = half * (rlv(i - 1, j, k) + rlv(i, j, k))
430  mulp = half * (rlv(i + 1, j, k) + rlv(i, j, k))
431  muem = half * (rev(i - 1, j, k) + rev(i, j, k))
432  muep = half * (rev(i + 1, j, k) + rev(i, j, k))
433 
434  c1m = ttm * (mulm + sig1 * muem) * rhoi
435  c1p = ttp * (mulp + sig1 * muep) * rhoi
436  c10 = c1m + c1p
437 
438  c2m = ttm
439  c2p = ttp
440  c20 = c2m + c2p
441 
442  ! Update the residual for this cell and store the possible
443  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
444 
445  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i - 1, j, k, itu3) &
446  - c10 * w(i, j, k, itu3) + c1p * w(i + 1, j, k, itu3)
447  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i - 1, j, k, itu4) &
448  - c20 * w(i, j, k, itu4) + c2p * w(i + 1, j, k, itu4)
449 
450  b1 = -c1m
451  c1 = c10
452  d1 = -c1p
453 
454  b2 = -c2m
455  c2 = c20
456  d2 = -c2p
457 
458  ! Update the central jacobian. For nonboundary cells this
459  ! is simply c1 and c2. For boundary cells this is slightly
460  ! more complicated, because the boundary conditions are
461  ! treated implicitly and the off-diagonal terms b1, b2 and
462  ! d1, d2 must be taken into account.
463  ! The boundary conditions are only treated implicitly if
464  ! the diagonal dominance of the matrix is increased.
465 
466  if (i == 2) then
467  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
468  - b1 * max(bmti1(j, k, itu3, itu3), zero)
469  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmti1(j, k, itu3, itu4)
470  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmti1(j, k, itu4, itu3)
471  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
472  - b2 * max(bmti1(j, k, itu4, itu4), zero)
473  else if (i == il) then
474  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
475  - d1 * max(bmti2(j, k, itu3, itu3), zero)
476  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmti2(j, k, itu3, itu4)
477  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmti2(j, k, itu4, itu3)
478  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
479  - d2 * max(bmti2(j, k, itu4, itu4), zero)
480  else
481  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
482  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
483  end if
484 
485  end do
486  end do
487  end do
488 
489  ! Multiply the residual by the volume and store this in dw; this
490  ! is done for monitoring reasons only. The multiplication with the
491  ! volume is present to be consistent with the flow residuals; also
492  ! the negative value is taken, again to be consistent with the
493  ! flow equations. Also multiply by iblank so that no updates occur
494  ! in holes or the overset boundary.
495 
496  do k = 2, kl
497  do j = 2, jl
498  do i = 2, il
499  rblank = real(iblank(i, j, k), realtype)
500  dw(i, j, k, itu3) = -volref(i, j, k) * dvt(i, j, k, 1) * rblank
501  dw(i, j, k, itu4) = -volref(i, j, k) * dvt(i, j, k, 2) * rblank
502  end do
503  end do
504  end do
505 
506  ! Initialize the wall function flags to .false.
507 
508  flagi2 = .false.
509  flagil = .false.
510  flagj2 = .false.
511  flagjl = .false.
512  flagk2 = .false.
513  flagkl = .false.
514 
515  ! Modify the rhs of the 1st internal cell, if wall functions
516  ! are used; their value is determined by the table.
517 
518  testwallfunctions: if (wallfunctions) then
519 
520  bocos: do nn = 1, nviscbocos
521 
522  ! Determine the block face on which the subface is located
523  ! and set some variables. As flag points to the entire array
524  ! flagI2, etc., its starting indices are the starting indices
525  ! of its target and not 1.
526 
527  select case (bcfaceid(nn))
528  case (imin)
529  flag => flagi2
530  dw2 => dw(2, 1:, 1:, 1:); dvt2 => dvt(2, 1:, 1:, 1:)
531  w2 => w(2, 1:, 1:, 1:); rlv2 => rlv(2, 1:, 1:)
532  w3 => w(3, 1:, 1:, 1:); rlv3 => rlv(3, 1:, 1:)
533  d2wall2 => d2wall(2, :, :); rev2 => rev(2, 1:, 1:)
534  d2wall3 => d2wall(3, :, :); rev3 => rev(3, 1:, 1:)
535 
536  case (imax)
537  flag => flagil
538  dw2 => dw(il, 1:, 1:, 1:); dvt2 => dvt(il, 1:, 1:, 1:)
539  w2 => w(il, 1:, 1:, 1:); rlv2 => rlv(il, 1:, 1:)
540  w3 => w(il - 1, 1:, 1:, 1:); rlv3 => rlv(il - 1, 1:, 1:)
541  d2wall2 => d2wall(il, :, :); rev2 => rev(il, 1:, 1:)
542  d2wall3 => d2wall(il - 1, :, :); rev3 => rev(il - 1, 1:, 1:)
543 
544  case (jmin)
545  flag => flagj2
546  dw2 => dw(1:, 2, 1:, 1:); dvt2 => dvt(1:, 2, 1:, 1:)
547  w2 => w(1:, 2, 1:, 1:); rlv2 => rlv(1:, 2, 1:)
548  w3 => w(1:, 3, 1:, 1:); rlv3 => rlv(1:, 3, 1:)
549  d2wall2 => d2wall(:, 2, :); rev2 => rev(1:, 2, 1:)
550  d2wall3 => d2wall(:, 3, :); rev3 => rev(1:, 3, 1:)
551 
552  case (jmax)
553  flag => flagjl
554  dw2 => dw(1:, jl, 1:, 1:); dvt2 => dvt(1:, jl, 1:, 1:)
555  w2 => w(1:, jl, 1:, 1:); rlv2 => rlv(1:, jl, 1:)
556  w3 => w(1:, jl - 1, 1:, 1:); rlv3 => rlv(1:, jl - 1, 1:)
557  d2wall2 => d2wall(:, jl, :); rev2 => rev(1:, jl, 1:)
558  d2wall3 => d2wall(:, jl - 1, :); rev3 => rev(1:, jl - 1, 1:)
559 
560  case (kmin)
561  flag => flagk2
562  dw2 => dw(1:, 1:, 2, 1:); dvt2 => dvt(1:, 1:, 2, 1:)
563  w2 => w(1:, 1:, 2, 1:); rlv2 => rlv(1:, 1:, 2)
564  w3 => w(1:, 1:, 3, 1:); rlv3 => rlv(1:, 1:, 3)
565  d2wall2 => d2wall(:, :, 2); rev2 => rev(1:, 1:, 2)
566  d2wall3 => d2wall(:, :, 3); rev3 => rev(1:, 1:, 3)
567 
568  case (kmax)
569  flag => flagkl
570  dw2 => dw(1:, 1:, kl, 1:); dvt2 => dvt(1:, 1:, kl, 1:)
571  w2 => w(1:, 1:, kl, 1:); rlv2 => rlv(1:, 1:, kl)
572  w3 => w(1:, 1:, kl - 1, 1:); rlv3 => rlv(1:, 1:, kl - 1)
573  d2wall2 => d2wall(:, :, kl); rev2 => rev(1:, 1:, kl)
574  d2wall3 => d2wall(:, :, kl - 1); rev3 => rev(1:, 1:, kl - 1)
575 
576  end select
577 
578  ! Loop over the owned faces of this subface. Therefore the
579  ! nodal range of BCData must be used. The offset of +1 is
580  ! present, because the starting index of the cell range is
581  ! 1 larger than the starting index of the nodal range.
582 
583  do j = (bcdata(nn)%jnBeg + 1), bcdata(nn)%jnEnd
584  do i = (bcdata(nn)%inBeg + 1), bcdata(nn)%inEnd
585 
586  ! Enforce v and f in the 1st internal cell from
587  ! the wall function table. There is an offset of -1 in
588  ! the wall distance. Note that the offset compared to
589  ! the current value must be stored. Also note that the
590  ! curve fits contain the non-dimensional values.
591 
592  utau = viscsubface(nn)%utau(i, j)
593  yp = w2(i, j, irho) * d2wall2(i - 1, j - 1) * utau / rlv2(i, j)
594 
595  ! Set dw2 to zero for proper monitoring of the
596  ! convergence.
597 
598  dw2(i, j, itu3) = zero
599  dw2(i, j, itu4) = zero
600 
601  ! Get table values
602 
603  call curvetupyp(tup(itu3:itu4), yp, itu3, itu4)
604  tu32 = tup(itu3) * utau**2
605  tu42 = tup(itu4) * utau**2 / rlv2(i, j) * w2(i, j, irho)
606 
607  ! Compute f from balance
608 
609  if (rvfn .eq. 1) then
610  call curvetupyp(tup(itu1:itu2), yp, itu1, itu2)
611  tu12 = tup(itu1) * utau**2
612  tu22 = tup(itu2) * utau**4 / rlv2(i, j) * w2(i, j, irho)
613  call curvetupyp(tup(itu5:itu5), yp, itu5, itu5)
614  tu52 = tup(itu5) * rlv2(i, j)
615  dtu23 = (w3(i, j, itu3) - tu32) &
616  / (d2wall3(i - 1, j - 1) - d2wall2(i - 1, j - 1))
617  rnu23 = half * ((tu52 + rlv2(i, j)) / w2(i, j, irho) + &
618  (rev3(i, j) + rlv3(i, j)) / w3(i, j, irho))
619  tu42 = (tu22 * tu32 / tu12 - rnu23 * dtu23 &
620  / (two * d2wall2(i - 1, j - 1))) / tu12
621  end if
622 
623  ! Set rhs to turbulence variables divide by betaTurb
624  ! because the update is scaled by betaTurb.
625  ! (see end of routine)
626 
627  dvt2(i, j, 1) = (tu32 - w2(i, j, itu3)) / betaturb
628  dvt2(i, j, 2) = (tu42 - w2(i, j, itu4)) / betaturb
629  if (rvfn .eq. 1) dvt2(i, j, 2) = (tu42 - w2(i, j, itu4)) * 0.01
630 
631  ! Set the wall flag to .true.
632 
633  flag(i, j) = .true.
634 
635  end do
636  end do
637 
638  end do bocos
639  end if testwallfunctions
640 
641  ! Return if only the residual must be computed.
642 
643  if (resonly) return
644 
645  do k = 2, kl
646  do j = 2, jl
647  do i = 2, il
648 
649  ! Set qq to 1 if the value is determined by the table.
650 
651  if ((i == 2 .and. flagi2(j, k)) .or. &
652  (i == il .and. flagil(j, k)) .or. &
653  (j == 2 .and. flagj2(i, k)) .or. &
654  (j == jl .and. flagjl(i, k)) .or. &
655  (k == 2 .and. flagk2(i, j)) .or. &
656  (k == kl .and. flagkl(i, j))) then
657  qq(i, j, k, 1, 1) = one
658  qq(i, j, k, 1, 2) = zero
659  qq(i, j, k, 2, 1) = zero
660  qq(i, j, k, 2, 2) = one
661  end if
662  end do
663  end do
664  end do
665 
666  ! Initialize the grid velocity to zero. This value will be used
667  ! if the block is not moving.
668 
669  qs = zero
670  !
671  ! dd-ADI step in j-direction. There is no particular reason to
672  ! start in j-direction, it just happened to be so. As we solve
673  ! in j-direction, the j-loop is the innermost loop.
674  !
675  do k = 2, kl
676  do i = 2, il
677  do j = 2, jl
678 
679  ! More or less the same code is executed here as above when
680  ! the residual was built. However, now the off-diagonal
681  ! terms for the dd-ADI must be built and stored. This could
682  ! have been done earlier, but then all the coefficients had
683  ! to be stored. To save memory, they are recomputed.
684  ! Consequently, see the j-loop to build the residual for
685  ! the comments.
686 
687  voli = one / vol(i, j, k)
688  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
689  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
690 
691  xm = sj(i, j - 1, k, 1) * volmi
692  ym = sj(i, j - 1, k, 2) * volmi
693  zm = sj(i, j - 1, k, 3) * volmi
694  xp = sj(i, j, k, 1) * volpi
695  yp = sj(i, j, k, 2) * volpi
696  zp = sj(i, j, k, 3) * volpi
697 
698  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
699  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
700  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
701  ttm = xm * xa + ym * ya + zm * za
702  ttp = xp * xa + yp * ya + zp * za
703 
704  ! Off-diagonal terms due to the diffusion terms
705  ! in j-direction.
706 
707  rhoi = one / w(i, j, k, irho)
708  mulm = half * (rlv(i, j - 1, k) + rlv(i, j, k))
709  mulp = half * (rlv(i, j + 1, k) + rlv(i, j, k))
710  muem = half * (rev(i, j - 1, k) + rev(i, j, k))
711  muep = half * (rev(i, j + 1, k) + rev(i, j, k))
712 
713  c1m = ttm * (mulm + sig1 * muem) * rhoi
714  c1p = ttp * (mulp + sig1 * muep) * rhoi
715 
716  c2m = ttm
717  c2p = ttp
718 
719  bb(1, j) = -c1m
720  dd(1, j) = -c1p
721  bb(2, j) = -c2m
722  dd(2, j) = -c2p
723 
724  ! Compute the grid velocity if present.
725  ! It is taken as the average of j and j-1,
726 
727  if (addgridvelocities) &
728  qs = half * (sfacej(i, j, k) + sfacej(i, j - 1, k)) * voli
729 
730  ! Off-diagonal terms due to the advection term in
731  ! j-direction. First order approximation.
732 
733  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
734  um = zero
735  up = zero
736  if (uu < zero) um = uu
737  if (uu > zero) up = uu
738 
739  bb(1, j) = bb(1, j) - up
740  dd(1, j) = dd(1, j) + um
741  bb(2, j) = bb(2, j)
742  dd(2, j) = dd(2, j)
743 
744  ! Store the central jacobian and rhs in cc and ff.
745  ! Multiply the off-diagonal terms and rhs by the iblank
746  ! value so the update determined for iblank = 0 is zero.
747 
748  rblank = real(iblank(i, j, k), realtype)
749 
750  cc(1, 1, j) = qq(i, j, k, 1, 1)
751  cc(1, 2, j) = qq(i, j, k, 1, 2) * rblank
752  cc(2, 1, j) = qq(i, j, k, 2, 1) * rblank
753  cc(2, 2, j) = qq(i, j, k, 2, 2)
754 
755  ff(1, j) = dvt(i, j, k, 1) * rblank
756  ff(2, j) = dvt(i, j, k, 2) * rblank
757 
758  bb(:, j) = bb(:, j) * rblank
759  dd(:, j) = dd(:, j) * rblank
760 
761  ! Set off diagonal terms to zero if wall function are used.
762 
763  if ((i == 2 .and. flagi2(j, k)) .or. &
764  (i == il .and. flagil(j, k)) .or. &
765  (j == 2 .and. flagj2(i, k)) .or. &
766  (j == jl .and. flagjl(i, k)) .or. &
767  (k == 2 .and. flagk2(i, j)) .or. &
768  (k == kl .and. flagkl(i, j))) then
769  bb(1, j) = zero
770  dd(1, j) = zero
771  bb(2, j) = zero
772  dd(2, j) = zero
773  end if
774 
775  end do
776 
777  ! Solve the tri-diagonal system in j-direction.
778 
779  call tdia3(2_inttype, jl, bb, cc, dd, ff)
780 
781  ! Determine the new rhs for the next direction.
782 
783  do j = 2, jl
784  dvt(i, j, k, 1) = qq(i, j, k, 1, 1) * ff(1, j) + qq(i, j, k, 1, 2) * ff(2, j)
785  dvt(i, j, k, 2) = qq(i, j, k, 2, 1) * ff(1, j) + qq(i, j, k, 2, 2) * ff(2, j)
786  end do
787 
788  end do
789  end do
790  !
791  ! dd-ADI step in i-direction. As we solve in i-direction, the
792  ! i-loop is the innermost loop.
793  !
794  do k = 2, kl
795  do j = 2, jl
796  do i = 2, il
797 
798  ! More or less the same code is executed here as above when
799  ! the residual was built. However, now the off-diagonal
800  ! terms for the dd-ADI must be built and stored. This could
801  ! have been done earlier, but then all the coefficients had
802  ! to be stored. To save memory, they are recomputed.
803  ! Consequently, see the i-loop to build the residual for
804  ! the comments.
805 
806  voli = one / vol(i, j, k)
807  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
808  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
809 
810  xm = si(i - 1, j, k, 1) * volmi
811  ym = si(i - 1, j, k, 2) * volmi
812  zm = si(i - 1, j, k, 3) * volmi
813  xp = si(i, j, k, 1) * volpi
814  yp = si(i, j, k, 2) * volpi
815  zp = si(i, j, k, 3) * volpi
816 
817  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
818  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
819  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
820  ttm = xm * xa + ym * ya + zm * za
821  ttp = xp * xa + yp * ya + zp * za
822 
823  ! Off-diagonal terms due to the diffusion terms
824  ! in i-direction.
825 
826  rhoi = one / w(i, j, k, irho)
827  mulm = half * (rlv(i - 1, j, k) + rlv(i, j, k))
828  mulp = half * (rlv(i + 1, j, k) + rlv(i, j, k))
829  muem = half * (rev(i - 1, j, k) + rev(i, j, k))
830  muep = half * (rev(i + 1, j, k) + rev(i, j, k))
831 
832  c1m = ttm * (mulm + sig1 * muem) * rhoi
833  c1p = ttp * (mulp + sig1 * muep) * rhoi
834 
835  c2m = ttm
836  c2p = ttp
837 
838  bb(1, i) = -c1m
839  dd(1, i) = -c1p
840  bb(2, i) = -c2m
841  dd(2, i) = -c2p
842 
843  ! Compute the grid velocity if present.
844  ! It is taken as the average of i and i-1,
845 
846  if (addgridvelocities) &
847  qs = half * (sfacei(i, j, k) + sfacei(i - 1, j, k)) * voli
848 
849  ! Off-diagonal terms due to the advection term in
850  ! i-direction. First order approximation.
851 
852  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
853  um = zero
854  up = zero
855  if (uu < zero) um = uu
856  if (uu > zero) up = uu
857 
858  bb(1, i) = bb(1, i) - up
859  dd(1, i) = dd(1, i) + um
860  bb(2, i) = bb(2, i)
861  dd(2, i) = dd(2, i)
862 
863  ! Store the central jacobian and rhs in cc and ff.
864  ! Multiply the off-diagonal terms and rhs by the iblank
865  ! value so the update determined for iblank = 0 is zero.
866 
867  rblank = real(iblank(i, j, k), realtype)
868 
869  cc(1, 1, i) = qq(i, j, k, 1, 1)
870  cc(1, 2, i) = qq(i, j, k, 1, 2) * rblank
871  cc(2, 1, i) = qq(i, j, k, 2, 1) * rblank
872  cc(2, 2, i) = qq(i, j, k, 2, 2)
873 
874  ff(1, i) = dvt(i, j, k, 1) * rblank
875  ff(2, i) = dvt(i, j, k, 2) * rblank
876 
877  bb(:, i) = bb(:, i) * rblank
878  dd(:, i) = dd(:, i) * rblank
879 
880  ! Set off diagonal terms to zero if wall function are used.
881 
882  if ((i == 2 .and. flagi2(j, k)) .or. &
883  (i == il .and. flagil(j, k)) .or. &
884  (j == 2 .and. flagj2(i, k)) .or. &
885  (j == jl .and. flagjl(i, k)) .or. &
886  (k == 2 .and. flagk2(i, j)) .or. &
887  (k == kl .and. flagkl(i, j))) then
888  bb(1, i) = zero
889  dd(1, i) = zero
890  bb(2, i) = zero
891  dd(2, i) = zero
892  end if
893 
894  end do
895 
896  ! Solve the tri-diagonal system in i-direction.
897 
898  call tdia3(2_inttype, il, bb, cc, dd, ff)
899 
900  ! Determine the new rhs for the next direction.
901 
902  do i = 2, il
903  dvt(i, j, k, 1) = qq(i, j, k, 1, 1) * ff(1, i) + qq(i, j, k, 1, 2) * ff(2, i)
904  dvt(i, j, k, 2) = qq(i, j, k, 2, 1) * ff(1, i) + qq(i, j, k, 2, 2) * ff(2, i)
905  end do
906 
907  end do
908  end do
909  !
910  ! dd-ADI step in k-direction. As we solve in k-direction, the
911  ! k-loop is the innermost loop.
912  !
913  do j = 2, jl
914  do i = 2, il
915  do k = 2, kl
916 
917  ! More or less the same code is executed here as above when
918  ! the residual was built. However, now the off-diagonal
919  ! terms for the dd-ADI must be built and stored. This could
920  ! have been done earlier, but then all the coefficients had
921  ! to be stored. To save memory, they are recomputed.
922  ! Consequently, see the k-loop to build the residual for
923  ! the comments.
924 
925  voli = one / vol(i, j, k)
926  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
927  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
928 
929  xm = sk(i, j, k - 1, 1) * volmi
930  ym = sk(i, j, k - 1, 2) * volmi
931  zm = sk(i, j, k - 1, 3) * volmi
932  xp = sk(i, j, k, 1) * volpi
933  yp = sk(i, j, k, 2) * volpi
934  zp = sk(i, j, k, 3) * volpi
935 
936  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
937  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
938  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
939  ttm = xm * xa + ym * ya + zm * za
940  ttp = xp * xa + yp * ya + zp * za
941 
942  ! Off-diagonal terms due to the diffusion terms
943  ! in k-direction.
944 
945  rhoi = one / w(i, j, k, irho)
946  mulm = half * (rlv(i, j, k - 1) + rlv(i, j, k))
947  mulp = half * (rlv(i, j, k + 1) + rlv(i, j, k))
948  muem = half * (rev(i, j, k - 1) + rev(i, j, k))
949  muep = half * (rev(i, j, k + 1) + rev(i, j, k))
950 
951  c1m = ttm * (mulm + sig1 * muem) * rhoi
952  c1p = ttp * (mulp + sig1 * muep) * rhoi
953 
954  c2m = ttm
955  c2p = ttp
956 
957  bb(1, k) = -c1m
958  dd(1, k) = -c1p
959  bb(2, k) = -c2m
960  dd(2, k) = -c2p
961 
962  ! Compute the grid velocity if present.
963  ! It is taken as the average of k and k-1,
964 
965  if (addgridvelocities) &
966  qs = half * (sfacek(i, j, k) + sfacek(i, j, k - 1)) * voli
967 
968  ! Off-diagonal terms due to the advection term in
969  ! k-direction. First order approximation.
970 
971  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
972  um = zero
973  up = zero
974  if (uu < zero) um = uu
975  if (uu > zero) up = uu
976 
977  bb(1, k) = bb(1, k) - up
978  dd(1, k) = dd(1, k) + um
979  bb(2, k) = bb(2, k)
980  dd(2, k) = dd(2, k)
981 
982  ! Store the central jacobian and rhs in cc and ff.
983  ! Multiply the off-diagonal terms and rhs by the iblank
984  ! value so the update determined for iblank = 0 is zero.
985 
986  rblank = real(iblank(i, j, k), realtype)
987 
988  cc(1, 1, k) = qq(i, j, k, 1, 1)
989  cc(1, 2, k) = qq(i, j, k, 1, 2) * rblank
990  cc(2, 1, k) = qq(i, j, k, 2, 1) * rblank
991  cc(2, 2, k) = qq(i, j, k, 2, 2)
992 
993  ff(1, k) = dvt(i, j, k, 1) * rblank
994  ff(2, k) = dvt(i, j, k, 2) * rblank
995 
996  bb(:, k) = bb(:, k) * rblank
997  dd(:, k) = dd(:, k) * rblank
998 
999  ! Set off diagonal terms to zero if wall function are used.
1000 
1001  if ((i == 2 .and. flagi2(j, k)) .or. &
1002  (i == il .and. flagil(j, k)) .or. &
1003  (j == 2 .and. flagj2(i, k)) .or. &
1004  (j == jl .and. flagjl(i, k)) .or. &
1005  (k == 2 .and. flagk2(i, j)) .or. &
1006  (k == kl .and. flagkl(i, j))) then
1007  bb(1, k) = zero
1008  dd(1, k) = zero
1009  bb(2, k) = zero
1010  dd(2, k) = zero
1011  end if
1012 
1013  end do
1014 
1015  ! Solve the tri-diagonal system in k-direction.
1016 
1017  call tdia3(2_inttype, kl, bb, cc, dd, ff)
1018 
1019  ! Store the update in dvt.
1020 
1021  do k = 2, kl
1022  dvt(i, j, k, 1) = ff(1, k)
1023  dvt(i, j, k, 2) = ff(2, k)
1024  end do
1025 
1026  end do
1027  end do
1028  !
1029  ! Update the turbulent variables.
1030  !
1031  do k = 2, kl
1032  do j = 2, jl
1033  do i = 2, il
1034  w(i, j, k, itu3) = w(i, j, k, itu3) + betaturb * dvt(i, j, k, 1)
1035  w(i, j, k, itu4) = w(i, j, k, itu4) + betaturb * dvt(i, j, k, 2)
1036  end do
1037  end do
1038  end do
1039 
1040  end subroutine vfsolve
1041  subroutine kesolve(resOnly)
1042  !
1043  ! keSolve solves the k-eps transport equations of the v2-f model
1044  ! in a coupled manner using a diagonal dominant ADI-scheme.
1045  !
1046  use blockpointers
1047  use constants
1048  use flowvarrefstate
1049  use inputiteration
1050  use inputphysics
1051  use paramturb
1052  use turbmod, only: prod, dvt, sct, scl2, sig1, sig2
1053  use turbutils, only: turbadvection, unsteadyturbterm, tdia3
1054  use turbcurvefits, only: curvetupyp
1055  implicit none
1056  !
1057  ! Subroutine arguments.
1058  !
1059  logical, intent(in) :: resOnly
1060  !
1061  ! Local variables.
1062  !
1063  integer(kind=intType) :: i, j, k, nn
1064 
1065  real(kind=realtype) :: alp
1066  real(kind=realtype) :: rhoi, ss, spk, ff1, ff2, ff3
1067  real(kind=realtype) :: ss1, ss2, ss3
1068  real(kind=realtype) :: voli, volmi, volpi
1069  real(kind=realtype) :: xm, ym, zm, xp, yp, zp, xa, ya, za
1070  real(kind=realtype) :: ttm, ttp, mulm, mulp, muem, muep
1071  real(kind=realtype) :: c1m, c1p, c10, c2m, c2p, c20
1072  real(kind=realtype) :: b1, b2, c1, c2, d1, d2
1073  real(kind=realtype) :: qs, uu, um, up, utau
1074  real(kind=realtype) :: tke, tep, tv2, tkea, tepa, tv2a
1075  real(kind=realtype) :: tv2l, stei, sle2i
1076  real(kind=realtype) :: rnu23, tu12, tu22, tu52, prod2, dtu23
1077  real(kind=realtype) :: factor, rblank
1078 
1079  real(kind=realtype), dimension(itu1:itu5) :: tup
1080 
1081  real(kind=realtype), dimension(2:il, 2:jl, 2:kl, 2, 2) :: qq
1082  real(kind=realtype), dimension(2, 2:max(il, jl, kl)) :: bb, dd, ff
1083  real(kind=realtype), dimension(2, 2, 2:max(il, jl, kl)) :: cc
1084 
1085  real(kind=realtype), dimension(:, :, :), pointer :: dw2, dvt2, w2, w3
1086  real(kind=realtype), dimension(:, :), pointer :: rlv2, rlv3
1087  real(kind=realtype), dimension(:, :), pointer :: rev2, rev3
1088  real(kind=realtype), dimension(:, :), pointer :: d2wall2, d2wall3
1089 
1090  logical, dimension(2:jl, 2:kl), target :: flagI2, flagIl
1091  logical, dimension(2:il, 2:kl), target :: flagJ2, flagJl
1092  logical, dimension(2:il, 2:jl), target :: flagK2, flagKl
1093 
1094  logical, dimension(:, :), pointer :: flag
1095 
1096  ! Set model constants
1097 
1098  sig1 = rvfsigk1
1099  sig2 = rvfsige1
1100 
1101  ! Set the pointer dvt to the correct entries in dw.
1102 
1103  dvt => scratch(1:, 1:, 1:, idvt:)
1104  sct => scratch(1:, 1:, 1:, isct)
1105  scl2 => scratch(1:, 1:, 1:, iscl2)
1106  !
1107  ! Source terms.
1108  ! Determine the source term and its derivative w.r.t. k and
1109  ! epsilon for all internal cells of the block.
1110  !
1111  do k = 2, kl
1112  do j = 2, jl
1113  do i = 2, il
1114 
1115  ! Compute the source terms for both the k and the epsilon
1116  ! equation.
1117 
1118  rhoi = one / w(i, j, k, irho)
1119  tke = w(i, j, k, itu1)
1120  tep = w(i, j, k, itu2)
1121  tv2 = w(i, j, k, itu3)
1122  tkea = abs(tke)
1123  tepa = abs(tep)
1124  tv2a = abs(tv2)
1125  tv2l = max(tv2a, rvflimitk * 0.66666_realtype)
1126  stei = tepa / sct(i, j, k)
1127  sle2i = tepa**2 / scl2(i, j, k)
1128  ss = prod(i, j, k)
1129  spk = rev(i, j, k) * ss * rhoi
1130  spk = min(spk, pklim * tepa)
1131 
1132  ! alp = rvfN1A+rvfN1B &
1133  ! / (1.+(d2Wall(i,j,k)*.5*rvfCl)**2*sle2i)**4
1134  alp = rvfn1a + rvfn1b &
1135  / (1.+d2wall(i, j, k)**2*.25 * rvfcl**2 * tepa**2 / scl2(i, j, k))**4
1136  if (rvfn == 6) alp = rvfn6a * (1.+rvfn6b * sqrt(tkea / tv2l))
1137 
1138  ff1 = zero
1139  ff2 = -one
1140  ff3 = spk
1141 
1142  ss1 = zero
1143  ss2 = -rvfbeta * stei
1144  ss3 = alp * spk * stei
1145 
1146  dvt(i, j, k, 1) = ff1 * tke + ff2 * tep + ff3
1147  dvt(i, j, k, 2) = ss1 * tke + ss2 * tep + ss3
1148 
1149  ! Compute the source term jacobian. Note that only the
1150  ! destruction terms are linearized to increase the diagonal
1151  ! dominance of the matrix. Furthermore minus the source
1152  ! term jacobian is stored.
1153 
1154  qq(i, j, k, 1, 1) = -ff1
1155  qq(i, j, k, 1, 2) = -ff2
1156  qq(i, j, k, 2, 1) = -ss1
1157  qq(i, j, k, 2, 2) = -ss2
1158 
1159  end do
1160  end do
1161  end do
1162  !
1163  ! Advection and unsteady terms.
1164  !
1165  nn = itu1 - 1
1166  call turbadvection(2_inttype, 2_inttype, nn, qq)
1167 
1168  call unsteadyturbterm(2_inttype, 2_inttype, nn, qq)
1169  !
1170  ! Viscous terms in k-direction.
1171  !
1172  do k = 2, kl
1173  do j = 2, jl
1174  do i = 2, il
1175 
1176  ! Compute the metrics in zeta-direction, i.e. along the
1177  ! line k = constant.
1178 
1179  voli = one / vol(i, j, k)
1180  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
1181  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
1182 
1183  xm = sk(i, j, k - 1, 1) * volmi
1184  ym = sk(i, j, k - 1, 2) * volmi
1185  zm = sk(i, j, k - 1, 3) * volmi
1186  xp = sk(i, j, k, 1) * volpi
1187  yp = sk(i, j, k, 2) * volpi
1188  zp = sk(i, j, k, 3) * volpi
1189 
1190  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
1191  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
1192  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
1193  ttm = xm * xa + ym * ya + zm * za
1194  ttp = xp * xa + yp * ya + zp * za
1195 
1196  ! Computation of the viscous terms in zeta-direction; note
1197  ! that cross-derivatives are neglected, i.e. the mesh is
1198  ! assumed to be orthogonal.
1199  ! The second derivative in zeta-direction is constructed as
1200  ! the central difference of the first order derivatives, i.e.
1201  ! d^2/dzeta^2 = d/dzeta (d/dzeta k+1/2 - d/dzeta k-1/2).
1202  ! In this way the metric as well as the varying viscosity
1203  ! can be taken into account; the latter appears inside the
1204  ! d/dzeta derivative. The whole term is divided by rho to
1205  ! obtain the diffusion term for k and epsilon.
1206 
1207  rhoi = one / w(i, j, k, irho)
1208  mulm = half * (rlv(i, j, k - 1) + rlv(i, j, k))
1209  mulp = half * (rlv(i, j, k + 1) + rlv(i, j, k))
1210  muem = half * (rev(i, j, k - 1) + rev(i, j, k))
1211  muep = half * (rev(i, j, k + 1) + rev(i, j, k))
1212 
1213  c1m = ttm * (mulm + sig1 * muem) * rhoi
1214  c1p = ttp * (mulp + sig1 * muep) * rhoi
1215  c10 = c1m + c1p
1216 
1217  c2m = ttm * (mulm + sig2 * muem) * rhoi
1218  c2p = ttp * (mulp + sig2 * muep) * rhoi
1219  c20 = c2m + c2p
1220 
1221  ! Update the residual for this cell and store the possible
1222  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
1223 
1224  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i, j, k - 1, itu1) &
1225  - c10 * w(i, j, k, itu1) + c1p * w(i, j, k + 1, itu1)
1226  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i, j, k - 1, itu2) &
1227  - c20 * w(i, j, k, itu2) + c2p * w(i, j, k + 1, itu2)
1228 
1229  b1 = -c1m
1230  c1 = c10
1231  d1 = -c1p
1232 
1233  b2 = -c2m
1234  c2 = c20
1235  d2 = -c2p
1236 
1237  ! Update the central jacobian. For nonboundary cells this
1238  ! is simply c1 and c2. For boundary cells this is slightly
1239  ! more complicated, because the boundary conditions are
1240  ! treated implicitly and the off-diagonal terms b1, b2 and
1241  ! d1, d2 must be taken into account.
1242  ! The boundary conditions are only treated implicitly if
1243  ! the diagonal dominance of the matrix is increased.
1244 
1245  if (k == 2) then
1246  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
1247  - b1 * max(bmtk1(i, j, itu1, itu1), zero)
1248  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmtk1(i, j, itu1, itu2)
1249  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmtk1(i, j, itu2, itu1)
1250  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
1251  - b2 * max(bmtk1(i, j, itu2, itu2), zero)
1252  else if (k == kl) then
1253  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
1254  - d1 * max(bmtk2(i, j, itu1, itu1), zero)
1255  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmtk2(i, j, itu1, itu2)
1256  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmtk2(i, j, itu2, itu1)
1257  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
1258  - d2 * max(bmtk2(i, j, itu2, itu2), zero)
1259  else
1260  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
1261  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
1262  end if
1263 
1264  end do
1265  end do
1266  end do
1267  !
1268  ! Viscous terms in j-direction.
1269  !
1270  do k = 2, kl
1271  do j = 2, jl
1272  do i = 2, il
1273 
1274  ! Compute the metrics in eta-direction, i.e. along the
1275  ! line j = constant.
1276 
1277  voli = one / vol(i, j, k)
1278  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
1279  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
1280 
1281  xm = sj(i, j - 1, k, 1) * volmi
1282  ym = sj(i, j - 1, k, 2) * volmi
1283  zm = sj(i, j - 1, k, 3) * volmi
1284  xp = sj(i, j, k, 1) * volpi
1285  yp = sj(i, j, k, 2) * volpi
1286  zp = sj(i, j, k, 3) * volpi
1287 
1288  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
1289  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
1290  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
1291  ttm = xm * xa + ym * ya + zm * za
1292  ttp = xp * xa + yp * ya + zp * za
1293 
1294  ! Computation of the viscous terms in eta-direction; note
1295  ! that cross-derivatives are neglected, i.e. the mesh is
1296  ! assumed to be orthogonal.
1297  ! The second derivative in eta-direction is constructed as
1298  ! the central difference of the first order derivatives, i.e.
1299  ! d^2/deta^2 = d/deta (d/deta j+1/2 - d/deta j-1/2).
1300  ! In this way the metric as well as the varying viscosity
1301  ! can be taken into account; the latter appears inside the
1302  ! d/deta derivative. The whole term is divided by rho to
1303  ! obtain the diffusion term for k and epsilon.
1304 
1305  rhoi = one / w(i, j, k, irho)
1306  mulm = half * (rlv(i, j - 1, k) + rlv(i, j, k))
1307  mulp = half * (rlv(i, j + 1, k) + rlv(i, j, k))
1308  muem = half * (rev(i, j - 1, k) + rev(i, j, k))
1309  muep = half * (rev(i, j + 1, k) + rev(i, j, k))
1310 
1311  c1m = ttm * (mulm + sig1 * muem) * rhoi
1312  c1p = ttp * (mulp + sig1 * muep) * rhoi
1313  c10 = c1m + c1p
1314 
1315  c2m = ttm * (mulm + sig2 * muem) * rhoi
1316  c2p = ttp * (mulp + sig2 * muep) * rhoi
1317  c20 = c2m + c2p
1318 
1319  ! Update the residual for this cell and store the possible
1320  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
1321 
1322  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i, j - 1, k, itu1) &
1323  - c10 * w(i, j, k, itu1) + c1p * w(i, j + 1, k, itu1)
1324  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i, j - 1, k, itu2) &
1325  - c20 * w(i, j, k, itu2) + c2p * w(i, j + 1, k, itu2)
1326 
1327  b1 = -c1m
1328  c1 = c10
1329  d1 = -c1p
1330 
1331  b2 = -c2m
1332  c2 = c20
1333  d2 = -c2p
1334 
1335  ! Update the central jacobian. For nonboundary cells this
1336  ! is simply c1 and c2. For boundary cells this is slightly
1337  ! more complicated, because the boundary conditions are
1338  ! treated implicitly and the off-diagonal terms b1, b2 and
1339  ! d1, d2 must be taken into account.
1340  ! The boundary conditions are only treated implicitly if
1341  ! the diagonal dominance of the matrix is increased.
1342 
1343  if (j == 2) then
1344  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
1345  - b1 * max(bmtj1(i, k, itu1, itu1), zero)
1346  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmtj1(i, k, itu1, itu2)
1347  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmtj1(i, k, itu2, itu1)
1348  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
1349  - b2 * max(bmtj1(i, k, itu2, itu2), zero)
1350  else if (j == jl) then
1351  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
1352  - d1 * max(bmtj2(i, k, itu1, itu1), zero)
1353  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmtj2(i, k, itu1, itu2)
1354  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmtj2(i, k, itu2, itu1)
1355  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
1356  - d2 * max(bmtj2(i, k, itu2, itu2), zero)
1357  else
1358  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
1359  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
1360  end if
1361 
1362  end do
1363  end do
1364  end do
1365  !
1366  ! Viscous terms in i-direction.
1367  !
1368  do k = 2, kl
1369  do j = 2, jl
1370  do i = 2, il
1371 
1372  ! Compute the metrics in xi-direction, i.e. along the
1373  ! line i = constant.
1374 
1375  voli = one / vol(i, j, k)
1376  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
1377  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
1378 
1379  xm = si(i - 1, j, k, 1) * volmi
1380  ym = si(i - 1, j, k, 2) * volmi
1381  zm = si(i - 1, j, k, 3) * volmi
1382  xp = si(i, j, k, 1) * volpi
1383  yp = si(i, j, k, 2) * volpi
1384  zp = si(i, j, k, 3) * volpi
1385 
1386  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
1387  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
1388  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
1389  ttm = xm * xa + ym * ya + zm * za
1390  ttp = xp * xa + yp * ya + zp * za
1391 
1392  ! Computation of the viscous terms in xi-direction; note
1393  ! that cross-derivatives are neglected, i.e. the mesh is
1394  ! assumed to be orthogonal.
1395  ! The second derivative in xi-direction is constructed as
1396  ! the central difference of the first order derivatives, i.e.
1397  ! d^2/dxi^2 = d/dxi (d/dxi i+1/2 - d/dxi i-1/2).
1398  ! In this way the metric as well as the varying viscosity
1399  ! can be taken into account; the latter appears inside the
1400  ! d/dxi derivative. The whole term is divided by rho to
1401  ! obtain the diffusion term for k and epsilon.
1402 
1403  rhoi = one / w(i, j, k, irho)
1404  mulm = half * (rlv(i - 1, j, k) + rlv(i, j, k))
1405  mulp = half * (rlv(i + 1, j, k) + rlv(i, j, k))
1406  muem = half * (rev(i - 1, j, k) + rev(i, j, k))
1407  muep = half * (rev(i + 1, j, k) + rev(i, j, k))
1408 
1409  c1m = ttm * (mulm + sig1 * muem) * rhoi
1410  c1p = ttp * (mulp + sig1 * muep) * rhoi
1411  c10 = c1m + c1p
1412 
1413  c2m = ttm * (mulm + sig2 * muem) * rhoi
1414  c2p = ttp * (mulp + sig2 * muep) * rhoi
1415  c20 = c2m + c2p
1416 
1417  ! Update the residual for this cell and store the possible
1418  ! coefficients for the matrix in b1, b2, c1, c2, d1 and d2.
1419 
1420  dvt(i, j, k, 1) = dvt(i, j, k, 1) + c1m * w(i - 1, j, k, itu1) &
1421  - c10 * w(i, j, k, itu1) + c1p * w(i + 1, j, k, itu1)
1422  dvt(i, j, k, 2) = dvt(i, j, k, 2) + c2m * w(i - 1, j, k, itu2) &
1423  - c20 * w(i, j, k, itu2) + c2p * w(i + 1, j, k, itu2)
1424 
1425  b1 = -c1m
1426  c1 = c10
1427  d1 = -c1p
1428 
1429  b2 = -c2m
1430  c2 = c20
1431  d2 = -c2p
1432 
1433  ! Update the central jacobian. For nonboundary cells this
1434  ! is simply c1 and c2. For boundary cells this is slightly
1435  ! more complicated, because the boundary conditions are
1436  ! treated implicitly and the off-diagonal terms b1, b2 and
1437  ! d1, d2 must be taken into account.
1438  ! The boundary conditions are only treated implicitly if
1439  ! the diagonal dominance of the matrix is increased.
1440 
1441  if (i == 2) then
1442  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
1443  - b1 * max(bmti1(j, k, itu1, itu1), zero)
1444  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - b1 * bmti1(j, k, itu1, itu2)
1445  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - b2 * bmti1(j, k, itu2, itu1)
1446  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
1447  - b2 * max(bmti1(j, k, itu2, itu2), zero)
1448  else if (i == il) then
1449  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1 &
1450  - d1 * max(bmti2(j, k, itu1, itu1), zero)
1451  qq(i, j, k, 1, 2) = qq(i, j, k, 1, 2) - d1 * bmti2(j, k, itu1, itu2)
1452  qq(i, j, k, 2, 1) = qq(i, j, k, 2, 1) - d2 * bmti2(j, k, itu2, itu1)
1453  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2 &
1454  - d2 * max(bmti2(j, k, itu2, itu2), zero)
1455  else
1456  qq(i, j, k, 1, 1) = qq(i, j, k, 1, 1) + c1
1457  qq(i, j, k, 2, 2) = qq(i, j, k, 2, 2) + c2
1458  end if
1459 
1460  end do
1461  end do
1462  end do
1463 
1464  ! Multiply the residual by the volume and store this in dw; this
1465  ! is done for monitoring reasons only. The multiplication with the
1466  ! volume is present to be consistent with the flow residuals; also
1467  ! the negative value is taken, again to be consistent with the
1468  ! flow equations. Also multiply by iblank so that no updates occur
1469  ! in holes or the overset boundary.
1470 
1471  do k = 2, kl
1472  do j = 2, jl
1473  do i = 2, il
1474  rblank = real(iblank(i, j, k), realtype)
1475  dw(i, j, k, itu1) = -volref(i, j, k) * dvt(i, j, k, 1) * rblank
1476  dw(i, j, k, itu2) = -volref(i, j, k) * dvt(i, j, k, 2) * rblank
1477  end do
1478  end do
1479  end do
1480 
1481  ! Initialize the wall function flags to .false.
1482 
1483  flagi2 = .false.
1484  flagil = .false.
1485  flagj2 = .false.
1486  flagjl = .false.
1487  flagk2 = .false.
1488  flagkl = .false.
1489 
1490  ! Modify the rhs of the 1st internal cell, if wall functions
1491  ! are used; their value is determined by the table.
1492 
1493  testwallfunctions: if (wallfunctions) then
1494 
1495  bocos: do nn = 1, nviscbocos
1496 
1497  ! Determine the block face on which the subface is located
1498  ! and set some variables. As flag points to the entire array
1499  ! flagI2, etc., its starting indices are the starting indices
1500  ! of its target and not 1.
1501 
1502  select case (bcfaceid(nn))
1503  case (imin)
1504  flag => flagi2
1505  dw2 => dw(2, 1:, 1:, 1:); dvt2 => dvt(2, 1:, 1:, 1:)
1506  w2 => w(2, 1:, 1:, 1:); rlv2 => rlv(2, 1:, 1:)
1507  w3 => w(3, 1:, 1:, 1:); rlv3 => rlv(3, 1:, 1:)
1508  d2wall2 => d2wall(2, :, :); rev2 => rev(2, 1:, 1:)
1509  d2wall3 => d2wall(3, :, :); rev3 => rev(3, 1:, 1:)
1510 
1511  case (imax)
1512  flag => flagil
1513  dw2 => dw(il, 1:, 1:, 1:); dvt2 => dvt(il, 1:, 1:, 1:)
1514  w2 => w(il, 1:, 1:, 1:); rlv2 => rlv(il, 1:, 1:)
1515  w3 => w(il - 1, 1:, 1:, 1:); rlv3 => rlv(il - 1, 1:, 1:)
1516  d2wall2 => d2wall(il, :, :); rev2 => rev(il, 1:, 1:)
1517  d2wall3 => d2wall(il - 1, :, :); rev3 => rev(il - 1, 1:, 1:)
1518 
1519  case (jmin)
1520  flag => flagj2
1521  dw2 => dw(1:, 2, 1:, 1:); dvt2 => dvt(1:, 2, 1:, 1:)
1522  w2 => w(1:, 2, 1:, 1:); rlv2 => rlv(1:, 2, 1:)
1523  w3 => w(1:, 3, 1:, 1:); rlv3 => rlv(1:, 3, 1:)
1524  d2wall2 => d2wall(:, 2, :); rev2 => rev(1:, 2, 1:)
1525  d2wall3 => d2wall(:, 3, :); rev3 => rev(1:, 3, 1:)
1526 
1527  case (jmax)
1528  flag => flagjl
1529  dw2 => dw(1:, jl, 1:, 1:); dvt2 => dvt(1:, jl, 1:, 1:)
1530  w2 => w(1:, jl, 1:, 1:); rlv2 => rlv(1:, jl, 1:)
1531  w3 => w(1:, jl - 1, 1:, 1:); rlv3 => rlv(1:, jl - 1, 1:)
1532  d2wall2 => d2wall(:, jl, :); rev2 => rev(1:, jl, 1:)
1533  d2wall3 => d2wall(:, jl - 1, :); rev3 => rev(1:, jl - 1, 1:)
1534 
1535  case (kmin)
1536  flag => flagk2
1537  dw2 => dw(1:, 1:, 2, 1:); dvt2 => dvt(1:, 1:, 2, 1:)
1538  w2 => w(1:, 1:, 2, 1:); rlv2 => rlv(1:, 1:, 2)
1539  w3 => w(1:, 1:, 3, 1:); rlv3 => rlv(1:, 1:, 3)
1540  d2wall2 => d2wall(:, :, 2); rev2 => rev(1:, 1:, 2)
1541  d2wall3 => d2wall(:, :, 3); rev3 => rev(1:, 1:, 3)
1542 
1543  case (kmax)
1544  flag => flagkl
1545  dw2 => dw(1:, 1:, kl, 1:); dvt2 => dvt(1:, 1:, kl, 1:)
1546  w2 => w(1:, 1:, kl, 1:); rlv2 => rlv(1:, 1:, kl)
1547  w3 => w(1:, 1:, kl - 1, 1:); rlv3 => rlv(1:, 1:, kl - 1)
1548  d2wall2 => d2wall(:, :, kl); rev2 => rev(1:, 1:, kl)
1549  d2wall3 => d2wall(:, :, kl - 1); rev3 => rev(1:, 1:, kl - 1)
1550 
1551  end select
1552 
1553  ! Loop over the owned faces of this subface. Therefore the
1554  ! nodal range of BCData must be used. The offset of +1 is
1555  ! present, because the starting index of the cell range is
1556  ! 1 larger than the starting index of the nodal range.
1557 
1558  do j = (bcdata(nn)%jnBeg + 1), bcdata(nn)%jnEnd
1559  do i = (bcdata(nn)%inBeg + 1), bcdata(nn)%inEnd
1560 
1561  ! Enforce k and epsilon in the 1st internal cell from
1562  ! the wall function table. There is an offset of -1 in
1563  ! the wall distance. Note that the offset compared to
1564  ! the current value must be stored. Also note that the
1565  ! curve fits contain the non-dimensional values.
1566 
1567  utau = viscsubface(nn)%utau(i, j)
1568  yp = w2(i, j, irho) * d2wall2(i - 1, j - 1) * utau / rlv2(i, j)
1569 
1570  ! Set dw2 to zero for proper monitoring of the
1571  ! convergence.
1572 
1573  dw2(i, j, itu1) = zero
1574  dw2(i, j, itu2) = zero
1575 
1576  ! Get table values
1577 
1578  call curvetupyp(tup(itu1:itu2), yp, itu1, itu2)
1579  tu12 = tup(itu1) * utau**2
1580  tu22 = tup(itu2) * utau**4 / rlv2(i, j) * w2(i, j, irho)
1581 
1582  ! Compute epsilon from balance
1583 
1584  call curvetupyp(tup(itu5:itu5), yp, itu5, itu5)
1585  tu52 = tup(itu5) * rlv2(i, j)
1586  dtu23 = (w3(i, j, itu1) - tu12) &
1587  / (d2wall3(i - 1, j - 1) - d2wall2(i - 1, j - 1))
1588  rnu23 = half * ((tu52 + rlv2(i, j)) / w2(i, j, irho) + &
1589  (rev3(i, j) + rlv3(i, j)) / w3(i, j, irho))
1590  prod2 = tu52 / w2(i, j, irho) * (utau**2 * w2(i, j, irho) &
1591  / (rlv2(i, j) + tu52))**2
1592  tu22 = prod2 + rnu23 * dtu23 / (two * d2wall2(i - 1, j - 1))
1593 
1594  ! Set rhs to turbulence variables
1595 
1596  ! dvt2(i,j,1) = (tu12 - w2(i,j,itu1))/alfaTurb
1597  ! dvt2(i,j,2) = (tu22 - w2(i,j,itu2))/alfaTurb
1598 
1599  dvt2(i, j, 1) = (tu12 - w2(i, j, itu1))
1600  dvt2(i, j, 2) = (tu22 - w2(i, j, itu2)) * 0.01
1601 
1602  ! Set the wall flag to .true.
1603 
1604  flag(i, j) = .true.
1605 
1606  end do
1607  end do
1608 
1609  end do bocos
1610  end if testwallfunctions
1611 
1612  ! Return if only the residual must be computed.
1613 
1614  if (resonly) return
1615 
1616  ! Take the local time step into account. Use characteristic
1617  ! time stepping, i.e. a matrix time step is used, where
1618  ! dt is the inverse of the central jacobian times the cfl
1619  ! number. To avoid a matrix inversion the following system
1620  ! is solved. (I/dt + cc + bb + dd)*dw = rhs. Due to the
1621  ! matrix time stepping I/dt = cc/cfl. As in the rest of the
1622  ! algorithm only the modified central jacobian is used,
1623  ! stored it now.
1624 
1625  ! Compute the factor multiplying the central jacobian, which
1626  ! is 1 + 1/cfl.
1627 
1628  factor = one
1629  if (wallfunctions) factor = one / alfaturb
1630 
1631  do k = 2, kl
1632  do j = 2, jl
1633  do i = 2, il
1634  qq(i, j, k, 1, 1) = factor * qq(i, j, k, 1, 1)
1635  qq(i, j, k, 1, 2) = factor * qq(i, j, k, 1, 2)
1636  qq(i, j, k, 2, 1) = factor * qq(i, j, k, 2, 1)
1637  qq(i, j, k, 2, 2) = factor * qq(i, j, k, 2, 2)
1638 
1639  ! Set qq to 1 if the value is determined by the table.
1640 
1641  if ((i == 2 .and. flagi2(j, k)) .or. &
1642  (i == il .and. flagil(j, k)) .or. &
1643  (j == 2 .and. flagj2(i, k)) .or. &
1644  (j == jl .and. flagjl(i, k)) .or. &
1645  (k == 2 .and. flagk2(i, j)) .or. &
1646  (k == kl .and. flagkl(i, j))) then
1647  qq(i, j, k, 1, 1) = one
1648  qq(i, j, k, 1, 2) = zero
1649  qq(i, j, k, 2, 1) = zero
1650  qq(i, j, k, 2, 2) = one
1651  end if
1652  end do
1653  end do
1654  end do
1655 
1656  ! Initialize the grid velocity to zero. This value will be used
1657  ! if the block is not moving.
1658 
1659  qs = zero
1660  !
1661  ! dd-ADI step in j-direction. There is no particular reason to
1662  ! start in j-direction, it just happened to be so. As we solve
1663  ! in j-direction, the j-loop is the innermost loop.
1664  !
1665  do k = 2, kl
1666  do i = 2, il
1667  do j = 2, jl
1668 
1669  ! More or less the same code is executed here as above when
1670  ! the residual was built. However, now the off-diagonal
1671  ! terms for the dd-ADI must be built and stored. This could
1672  ! have been done earlier, but then all the coefficients had
1673  ! to be stored. To save memory, they are recomputed.
1674  ! Consequently, see the j-loop to build the residual for
1675  ! the comments.
1676 
1677  voli = one / vol(i, j, k)
1678  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
1679  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
1680 
1681  xm = sj(i, j - 1, k, 1) * volmi
1682  ym = sj(i, j - 1, k, 2) * volmi
1683  zm = sj(i, j - 1, k, 3) * volmi
1684  xp = sj(i, j, k, 1) * volpi
1685  yp = sj(i, j, k, 2) * volpi
1686  zp = sj(i, j, k, 3) * volpi
1687 
1688  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
1689  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
1690  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
1691  ttm = xm * xa + ym * ya + zm * za
1692  ttp = xp * xa + yp * ya + zp * za
1693 
1694  ! Off-diagonal terms due to the diffusion terms
1695  ! in j-direction.
1696 
1697  rhoi = one / w(i, j, k, irho)
1698  mulm = half * (rlv(i, j - 1, k) + rlv(i, j, k))
1699  mulp = half * (rlv(i, j + 1, k) + rlv(i, j, k))
1700  muem = half * (rev(i, j - 1, k) + rev(i, j, k))
1701  muep = half * (rev(i, j + 1, k) + rev(i, j, k))
1702 
1703  c1m = ttm * (mulm + sig1 * muem) * rhoi
1704  c1p = ttp * (mulp + sig1 * muep) * rhoi
1705 
1706  c2m = ttm * (mulm + sig2 * muem) * rhoi
1707  c2p = ttp * (mulp + sig2 * muep) * rhoi
1708 
1709  bb(1, j) = -c1m
1710  dd(1, j) = -c1p
1711  bb(2, j) = -c2m
1712  dd(2, j) = -c2p
1713 
1714  ! Compute the grid velocity if present.
1715  ! It is taken as the average of j and j-1,
1716 
1717  if (addgridvelocities) &
1718  qs = half * (sfacej(i, j, k) + sfacej(i, j - 1, k)) * voli
1719 
1720  ! Off-diagonal terms due to the advection term in
1721  ! j-direction. First order approximation.
1722 
1723  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1724  um = zero
1725  up = zero
1726  if (uu < zero) um = uu
1727  if (uu > zero) up = uu
1728 
1729  bb(1, j) = bb(1, j) - up
1730  dd(1, j) = dd(1, j) + um
1731  bb(2, j) = bb(2, j) - up
1732  dd(2, j) = dd(2, j) + um
1733 
1734  ! Store the central jacobian and rhs in cc and ff.
1735  ! Multiply the off-diagonal terms and rhs by the iblank
1736  ! value so the update determined for iblank = 0 is zero.
1737 
1738  rblank = real(iblank(i, j, k), realtype)
1739 
1740  cc(1, 1, j) = qq(i, j, k, 1, 1)
1741  cc(1, 2, j) = qq(i, j, k, 1, 2) * rblank
1742  cc(2, 1, j) = qq(i, j, k, 2, 1) * rblank
1743  cc(2, 2, j) = qq(i, j, k, 2, 2)
1744 
1745  ff(1, j) = dvt(i, j, k, 1) * rblank
1746  ff(2, j) = dvt(i, j, k, 2) * rblank
1747 
1748  bb(:, j) = bb(:, j) * rblank
1749  dd(:, j) = dd(:, j) * rblank
1750 
1751  ! Set off diagonal terms to zero if wall function are used.
1752 
1753  if ((i == 2 .and. flagi2(j, k)) .or. &
1754  (i == il .and. flagil(j, k)) .or. &
1755  (j == 2 .and. flagj2(i, k)) .or. &
1756  (j == jl .and. flagjl(i, k)) .or. &
1757  (k == 2 .and. flagk2(i, j)) .or. &
1758  (k == kl .and. flagkl(i, j))) then
1759  bb(1, j) = zero
1760  dd(1, j) = zero
1761  bb(2, j) = zero
1762  dd(2, j) = zero
1763  end if
1764 
1765  end do
1766 
1767  ! Solve the tri-diagonal system in j-direction.
1768 
1769  call tdia3(2_inttype, jl, bb, cc, dd, ff)
1770 
1771  ! Determine the new rhs for the next direction.
1772 
1773  do j = 2, jl
1774  dvt(i, j, k, 1) = qq(i, j, k, 1, 1) * ff(1, j) + qq(i, j, k, 1, 2) * ff(2, j)
1775  dvt(i, j, k, 2) = qq(i, j, k, 2, 1) * ff(1, j) + qq(i, j, k, 2, 2) * ff(2, j)
1776  end do
1777 
1778  end do
1779  end do
1780  !
1781  ! dd-ADI step in i-direction. As we solve in i-direction, the
1782  ! i-loop is the innermost loop.
1783  !
1784  do k = 2, kl
1785  do j = 2, jl
1786  do i = 2, il
1787 
1788  ! More or less the same code is executed here as above when
1789  ! the residual was built. However, now the off-diagonal
1790  ! terms for the dd-ADI must be built and stored. This could
1791  ! have been done earlier, but then all the coefficients had
1792  ! to be stored. To save memory, they are recomputed.
1793  ! Consequently, see the i-loop to build the residual for
1794  ! the comments.
1795 
1796  voli = one / vol(i, j, k)
1797  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
1798  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
1799 
1800  xm = si(i - 1, j, k, 1) * volmi
1801  ym = si(i - 1, j, k, 2) * volmi
1802  zm = si(i - 1, j, k, 3) * volmi
1803  xp = si(i, j, k, 1) * volpi
1804  yp = si(i, j, k, 2) * volpi
1805  zp = si(i, j, k, 3) * volpi
1806 
1807  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
1808  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
1809  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
1810  ttm = xm * xa + ym * ya + zm * za
1811  ttp = xp * xa + yp * ya + zp * za
1812 
1813  ! Off-diagonal terms due to the diffusion terms
1814  ! in i-direction.
1815 
1816  rhoi = one / w(i, j, k, irho)
1817  mulm = half * (rlv(i - 1, j, k) + rlv(i, j, k))
1818  mulp = half * (rlv(i + 1, j, k) + rlv(i, j, k))
1819  muem = half * (rev(i - 1, j, k) + rev(i, j, k))
1820  muep = half * (rev(i + 1, j, k) + rev(i, j, k))
1821 
1822  c1m = ttm * (mulm + sig1 * muem) * rhoi
1823  c1p = ttp * (mulp + sig1 * muep) * rhoi
1824 
1825  c2m = ttm * (mulm + sig2 * muem) * rhoi
1826  c2p = ttp * (mulp + sig2 * muep) * rhoi
1827 
1828  bb(1, i) = -c1m
1829  dd(1, i) = -c1p
1830  bb(2, i) = -c2m
1831  dd(2, i) = -c2p
1832 
1833  ! Compute the grid velocity if present.
1834  ! It is taken as the average of i and i-1,
1835 
1836  if (addgridvelocities) &
1837  qs = half * (sfacei(i, j, k) + sfacei(i - 1, j, k)) * voli
1838 
1839  ! Off-diagonal terms due to the advection term in
1840  ! i-direction. First order approximation.
1841 
1842  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1843  um = zero
1844  up = zero
1845  if (uu < zero) um = uu
1846  if (uu > zero) up = uu
1847 
1848  bb(1, i) = bb(1, i) - up
1849  dd(1, i) = dd(1, i) + um
1850  bb(2, i) = bb(2, i) - up
1851  dd(2, i) = dd(2, i) + um
1852 
1853  ! Store the central jacobian and rhs in cc and ff.
1854  ! Multiply the off-diagonal terms and rhs by the iblank
1855  ! value so the update determined for iblank = 0 is zero.
1856 
1857  rblank = real(iblank(i, j, k), realtype)
1858 
1859  cc(1, 1, i) = qq(i, j, k, 1, 1)
1860  cc(1, 2, i) = qq(i, j, k, 1, 2) * rblank
1861  cc(2, 1, i) = qq(i, j, k, 2, 1) * rblank
1862  cc(2, 2, i) = qq(i, j, k, 2, 2)
1863 
1864  ff(1, i) = dvt(i, j, k, 1) * rblank
1865  ff(2, i) = dvt(i, j, k, 2) * rblank
1866 
1867  bb(:, i) = bb(:, i) * rblank
1868  dd(:, i) = dd(:, i) * rblank
1869 
1870  ! Set off diagonal terms to zero if wall function are used.
1871 
1872  if ((i == 2 .and. flagi2(j, k)) .or. &
1873  (i == il .and. flagil(j, k)) .or. &
1874  (j == 2 .and. flagj2(i, k)) .or. &
1875  (j == jl .and. flagjl(i, k)) .or. &
1876  (k == 2 .and. flagk2(i, j)) .or. &
1877  (k == kl .and. flagkl(i, j))) then
1878  bb(1, i) = zero
1879  dd(1, i) = zero
1880  bb(2, i) = zero
1881  dd(2, i) = zero
1882  end if
1883 
1884  end do
1885 
1886  ! Solve the tri-diagonal system in i-direction.
1887 
1888  call tdia3(2_inttype, il, bb, cc, dd, ff)
1889 
1890  ! Determine the new rhs for the next direction.
1891 
1892  do i = 2, il
1893  dvt(i, j, k, 1) = qq(i, j, k, 1, 1) * ff(1, i) + qq(i, j, k, 1, 2) * ff(2, i)
1894  dvt(i, j, k, 2) = qq(i, j, k, 2, 1) * ff(1, i) + qq(i, j, k, 2, 2) * ff(2, i)
1895  end do
1896 
1897  end do
1898  end do
1899  !
1900  ! dd-ADI step in k-direction. As we solve in k-direction, the
1901  ! k-loop is the innermost loop.
1902  !
1903  do j = 2, jl
1904  do i = 2, il
1905  do k = 2, kl
1906 
1907  ! More or less the same code is executed here as above when
1908  ! the residual was built. However, now the off-diagonal
1909  ! terms for the dd-ADI must be built and stored. This could
1910  ! have been done earlier, but then all the coefficients had
1911  ! to be stored. To save memory, they are recomputed.
1912  ! Consequently, see the k-loop to build the residual for
1913  ! the comments.
1914 
1915  voli = one / vol(i, j, k)
1916  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
1917  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
1918 
1919  xm = sk(i, j, k - 1, 1) * volmi
1920  ym = sk(i, j, k - 1, 2) * volmi
1921  zm = sk(i, j, k - 1, 3) * volmi
1922  xp = sk(i, j, k, 1) * volpi
1923  yp = sk(i, j, k, 2) * volpi
1924  zp = sk(i, j, k, 3) * volpi
1925 
1926  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
1927  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
1928  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
1929  ttm = xm * xa + ym * ya + zm * za
1930  ttp = xp * xa + yp * ya + zp * za
1931 
1932  ! Off-diagonal terms due to the diffusion terms
1933  ! in k-direction.
1934 
1935  rhoi = one / w(i, j, k, irho)
1936  mulm = half * (rlv(i, j, k - 1) + rlv(i, j, k))
1937  mulp = half * (rlv(i, j, k + 1) + rlv(i, j, k))
1938  muem = half * (rev(i, j, k - 1) + rev(i, j, k))
1939  muep = half * (rev(i, j, k + 1) + rev(i, j, k))
1940 
1941  c1m = ttm * (mulm + sig1 * muem) * rhoi
1942  c1p = ttp * (mulp + sig1 * muep) * rhoi
1943 
1944  c2m = ttm * (mulm + sig2 * muem) * rhoi
1945  c2p = ttp * (mulp + sig2 * muep) * rhoi
1946 
1947  bb(1, k) = -c1m
1948  dd(1, k) = -c1p
1949  bb(2, k) = -c2m
1950  dd(2, k) = -c2p
1951 
1952  ! Compute the grid velocity if present.
1953  ! It is taken as the average of k and k-1,
1954 
1955  if (addgridvelocities) &
1956  qs = half * (sfacek(i, j, k) + sfacek(i, j, k - 1)) * voli
1957 
1958  ! Off-diagonal terms due to the advection term in
1959  ! k-direction. First order approximation.
1960 
1961  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1962  um = zero
1963  up = zero
1964  if (uu < zero) um = uu
1965  if (uu > zero) up = uu
1966 
1967  bb(1, k) = bb(1, k) - up
1968  dd(1, k) = dd(1, k) + um
1969  bb(2, k) = bb(2, k) - up
1970  dd(2, k) = dd(2, k) + um
1971 
1972  ! Store the central jacobian and rhs in cc and ff.
1973  ! Multiply the off-diagonal terms and rhs by the iblank
1974  ! value so the update determined for iblank = 0 is zero.
1975 
1976  rblank = real(iblank(i, j, k), realtype)
1977 
1978  cc(1, 1, k) = qq(i, j, k, 1, 1)
1979  cc(1, 2, k) = qq(i, j, k, 1, 2) * rblank
1980  cc(2, 1, k) = qq(i, j, k, 2, 1) * rblank
1981  cc(2, 2, k) = qq(i, j, k, 2, 2)
1982 
1983  ff(1, k) = dvt(i, j, k, 1) * rblank
1984  ff(2, k) = dvt(i, j, k, 2) * rblank
1985 
1986  bb(:, k) = bb(:, k) * rblank
1987  dd(:, k) = dd(:, k) * rblank
1988 
1989  ! Set off diagonal terms to zero if wall function are used.
1990 
1991  if ((i == 2 .and. flagi2(j, k)) .or. &
1992  (i == il .and. flagil(j, k)) .or. &
1993  (j == 2 .and. flagj2(i, k)) .or. &
1994  (j == jl .and. flagjl(i, k)) .or. &
1995  (k == 2 .and. flagk2(i, j)) .or. &
1996  (k == kl .and. flagkl(i, j))) then
1997  bb(1, k) = zero
1998  dd(1, k) = zero
1999  bb(2, k) = zero
2000  dd(2, k) = zero
2001  end if
2002 
2003  end do
2004 
2005  ! Solve the tri-diagonal system in k-direction.
2006 
2007  call tdia3(2_inttype, kl, bb, cc, dd, ff)
2008 
2009  ! Store the update in dvt.
2010 
2011  do k = 2, kl
2012  dvt(i, j, k, 1) = ff(1, k)
2013  dvt(i, j, k, 2) = ff(2, k)
2014  end do
2015 
2016  end do
2017  end do
2018  !
2019  ! Update the turbulent variables.
2020  !
2021  factor = alfaturb
2022  if (wallfunctions) factor = one
2023  do k = 2, kl
2024  do j = 2, jl
2025  do i = 2, il
2026  w(i, j, k, itu1) = w(i, j, k, itu1) + factor * dvt(i, j, k, 1)
2027  w(i, j, k, itu2) = w(i, j, k, itu2) + factor * dvt(i, j, k, 2)
2028  end do
2029  end do
2030  end do
2031 
2032  end subroutine kesolve
2033 
2034 end module vf
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::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::imax
integer(kind=inttype), parameter imax
Definition: constants.F90:293
constants::itu3
integer, parameter itu3
Definition: constants.F90:43
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::isct
integer, parameter isct
Definition: constants.F90:58
constants::itu1
integer, parameter itu1
Definition: constants.F90:40
constants::itu5
integer, parameter itu5
Definition: constants.F90:45
constants::irho
integer, parameter irho
Definition: constants.F90:34
constants::ivx
integer, parameter ivx
Definition: constants.F90:35
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
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::itu4
integer, parameter itu4
Definition: constants.F90:44
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
constants::iscl2
integer, parameter iscl2
Definition: constants.F90:59
inputiteration::betaturb
real(kind=realtype) betaturb
Definition: inputParam.F90:276
inputiteration::alfaturb
real(kind=realtype) alfaturb
Definition: inputParam.F90:276
inputphysics::wallfunctions
logical wallfunctions
Definition: inputParam.F90:589
inputphysics::pklim
real(kind=realtype) pklim
Definition: inputParam.F90:596
inputphysics::rvfn
integer(kind=inttype) rvfn
Definition: inputParam.F90:585
paramturb::rvfc2
real(kind=realtype), parameter rvfc2
Definition: paramTurb.F90:64
paramturb::rvfn6a
real(kind=realtype), parameter rvfn6a
Definition: paramTurb.F90:76
paramturb::rvfsigv1
real(kind=realtype), parameter rvfsigv1
Definition: paramTurb.F90:68
paramturb::rvfn6b
real(kind=realtype), parameter rvfn6b
Definition: paramTurb.F90:77
paramturb::rvfbeta
real(kind=realtype), parameter rvfbeta
Definition: paramTurb.F90:65
paramturb::rvfcl
real(kind=realtype) rvfcl
Definition: paramTurb.F90:80
paramturb::rvfsigk1
real(kind=realtype), parameter rvfsigk1
Definition: paramTurb.F90:66
paramturb::rvfsige1
real(kind=realtype), parameter rvfsige1
Definition: paramTurb.F90:67
paramturb::rvflimitk
real(kind=realtype) rvflimitk
Definition: paramTurb.F90:80
paramturb::rvfn1a
real(kind=realtype), parameter rvfn1a
Definition: paramTurb.F90:72
paramturb::rvfn1b
real(kind=realtype), parameter rvfn1b
Definition: paramTurb.F90:73
paramturb::rvfc1
real(kind=realtype), parameter rvfc1
Definition: paramTurb.F90:63
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::sig2
real(kind=realtype) sig2
Definition: turbMod.F90:16
turbmod::scl2
real(kind=realtype), dimension(:, :, :), pointer scl2
Definition: turbMod.F90:40
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::vfeddyviscosity
subroutine vfeddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:1635
turbutils::vfscale
subroutine vfscale
Definition: turbUtils.F90:1767
turbutils::unsteadyturbterm
subroutine unsteadyturbterm(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:411
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
utils
Definition: utils.F90:1
utils::setpointers
subroutine setpointers(nn, mm, ll)
Definition: utils.F90:3237
vf
Definition: vf.F90:1
vf::kesolve
subroutine kesolve(resOnly)
Definition: vf.F90:1042
vf::vf_block
subroutine vf_block(resOnly)
Definition: vf.F90:5
vf::vfsolve
subroutine vfsolve(resOnly)
Definition: vf.F90:61
Generated by doxygen 1.9.1