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