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