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