ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
turbUtils.F90
Go to the documentation of this file.
1 module turbutils
2 
3 contains
4 
5  subroutine prodkatolaunder
6  !
7  ! prodKatoLaunder computes the turbulent production term using
8  ! the Kato-Launder formulation.
9  !
10  use constants
11  use blockpointers, only: nx, ny, nz, il, jl, kl, w, si, sj, sk, vol, sectionid, scratch
12  use flowvarrefstate, only: timeref
13  use section, only: sections
14  use turbmod, only: prod
15  implicit none
16  !
17  ! Local variables.
18  !
19  integer(kind=intType) :: i, j, k, ii
20 
21  real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
22  real(kind=realtype) :: qxx, qyy, qzz, qxy, qxz, qyz, sijsij
23  real(kind=realtype) :: oxy, oxz, oyz, oijoij
24  real(kind=realtype) :: fact, omegax, omegay, omegaz
25 
26  ! Determine the non-dimensional wheel speed of this block.
27  ! The vorticity term, which appears in Kato-Launder is of course
28  ! not frame invariant. To approximate frame invariance the wheel
29  ! speed should be substracted from oxy, oxz and oyz, which results
30  ! in the vorticity in the rotating frame. However some people
31  ! claim that the absolute vorticity should be used to obtain the
32  ! best results. In that omega should be set to zero.
33 
34  omegax = timeref * sections(sectionid)%rotRate(1)
35  omegay = timeref * sections(sectionid)%rotRate(2)
36  omegaz = timeref * sections(sectionid)%rotRate(3)
37 
38  ! Loop over the cell centers of the given block. It may be more
39  ! efficient to loop over the faces and to scatter the gradient,
40  ! but in that case the gradients for u, v and w must be stored.
41  ! In the current approach no extra memory is needed.
42 #ifdef TAPENADE_REVERSE
43  !$AD II-LOOP
44  do ii = 0, nx * ny * nz - 1
45  i = mod(ii, nx) + 2
46  j = mod(ii / nx, ny) + 2
47  k = ii / (nx * ny) + 2
48 #else
49  do k = 2, kl
50  do j = 2, jl
51  do i = 2, il
52 #endif
53 
54  ! Compute the gradient of u in the cell center. Use is made
55  ! of the fact that the surrounding normals sum up to zero,
56  ! such that the cell i,j,k does not give a contribution.
57  ! The gradient is scaled by a factor 2*vol.
58 
59  uux = w(i + 1, j, k, ivx) * si(i, j, k, 1) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 1) &
60  + w(i, j + 1, k, ivx) * sj(i, j, k, 1) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 1) &
61  + w(i, j, k + 1, ivx) * sk(i, j, k, 1) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 1)
62  uuy = w(i + 1, j, k, ivx) * si(i, j, k, 2) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 2) &
63  + w(i, j + 1, k, ivx) * sj(i, j, k, 2) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 2) &
64  + w(i, j, k + 1, ivx) * sk(i, j, k, 2) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 2)
65  uuz = w(i + 1, j, k, ivx) * si(i, j, k, 3) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 3) &
66  + w(i, j + 1, k, ivx) * sj(i, j, k, 3) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 3) &
67  + w(i, j, k + 1, ivx) * sk(i, j, k, 3) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 3)
68 
69  ! Idem for the gradient of v.
70 
71  vvx = w(i + 1, j, k, ivy) * si(i, j, k, 1) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 1) &
72  + w(i, j + 1, k, ivy) * sj(i, j, k, 1) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 1) &
73  + w(i, j, k + 1, ivy) * sk(i, j, k, 1) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 1)
74  vvy = w(i + 1, j, k, ivy) * si(i, j, k, 2) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 2) &
75  + w(i, j + 1, k, ivy) * sj(i, j, k, 2) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 2) &
76  + w(i, j, k + 1, ivy) * sk(i, j, k, 2) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 2)
77  vvz = w(i + 1, j, k, ivy) * si(i, j, k, 3) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 3) &
78  + w(i, j + 1, k, ivy) * sj(i, j, k, 3) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 3) &
79  + w(i, j, k + 1, ivy) * sk(i, j, k, 3) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 3)
80 
81  ! And for the gradient of w.
82 
83  wwx = w(i + 1, j, k, ivz) * si(i, j, k, 1) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 1) &
84  + w(i, j + 1, k, ivz) * sj(i, j, k, 1) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 1) &
85  + w(i, j, k + 1, ivz) * sk(i, j, k, 1) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 1)
86  wwy = w(i + 1, j, k, ivz) * si(i, j, k, 2) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 2) &
87  + w(i, j + 1, k, ivz) * sj(i, j, k, 2) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 2) &
88  + w(i, j, k + 1, ivz) * sk(i, j, k, 2) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 2)
89  wwz = w(i + 1, j, k, ivz) * si(i, j, k, 3) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 3) &
90  + w(i, j + 1, k, ivz) * sj(i, j, k, 3) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 3) &
91  + w(i, j, k + 1, ivz) * sk(i, j, k, 3) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 3)
92 
93  ! Compute the strain and vorticity terms. The multiplication
94  ! is present to obtain the correct gradients. Note that
95  ! the wheel speed is substracted from the vorticity terms.
96 
97  fact = half / vol(i, j, k)
98 
99  qxx = fact * uux
100  qyy = fact * vvy
101  qzz = fact * wwz
102 
103  qxy = fact * half * (uuy + vvx)
104  qxz = fact * half * (uuz + wwx)
105  qyz = fact * half * (vvz + wwy)
106 
107  oxy = fact * half * (vvx - uuy) - omegaz
108  oxz = fact * half * (uuz - wwx) - omegay
109  oyz = fact * half * (wwy - vvz) - omegax
110 
111  ! Compute the summation of the strain and vorticity tensors.
112 
113  sijsij = two * (qxy**2 + qxz**2 + qyz**2) &
114  + qxx**2 + qyy**2 + qzz**2
115  oijoij = two * (oxy**2 + oxz**2 + oyz**2)
116 
117  ! Compute the production term.
118 
119  scratch(i, j, k, iprod) = two * sqrt(sijsij * oijoij)
120 #ifdef TAPENADE_REVERSE
121  end do
122 #else
123  end do
124  end do
125  end do
126 #endif
127  end subroutine prodkatolaunder
128 
129  subroutine prodsmag2
130  !
131  ! prodSmag2 computes the term:
132  ! 2*sij*sij - 2/3 div(u)**2 with sij=0.5*(duidxj+dujdxi)
133  ! which is used for the turbulence equations.
134  ! It is assumed that the pointer prod, stored in turbMod, is
135  ! already set to the correct entry.
136  !
137  use constants
138  use blockpointers, only: nx, ny, nz, il, jl, kl, w, si, sj, sk, vol, sectionid, scratch
139  implicit none
140  !
141  ! Local parameter
142  !
143  real(kind=realtype), parameter :: f23 = two * third
144  !
145  ! Local variables.
146  !
147  integer(kind=intType) :: i, j, k, ii
148  real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
149  real(kind=realtype) :: div2, fact, sxx, syy, szz, sxy, sxz, syz
150 
151  ! Loop over the cell centers of the given block. It may be more
152  ! efficient to loop over the faces and to scatter the gradient,
153  ! but in that case the gradients for u, v and w must be stored.
154  ! In the current approach no extra memory is needed.
155 
156 #ifdef TAPENADE_REVERSE
157  !$AD II-LOOP
158  do ii = 0, nx * ny * nz - 1
159  i = mod(ii, nx) + 2
160  j = mod(ii / nx, ny) + 2
161  k = ii / (nx * ny) + 2
162 #else
163  do k = 2, kl
164  do j = 2, jl
165  do i = 2, il
166 #endif
167 
168  ! Compute the gradient of u in the cell center. Use is made
169  ! of the fact that the surrounding normals sum up to zero,
170  ! such that the cell i,j,k does not give a contribution.
171  ! The gradient is scaled by the factor 2*vol.
172 
173  uux = w(i + 1, j, k, ivx) * si(i, j, k, 1) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 1) &
174  + w(i, j + 1, k, ivx) * sj(i, j, k, 1) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 1) &
175  + w(i, j, k + 1, ivx) * sk(i, j, k, 1) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 1)
176  uuy = w(i + 1, j, k, ivx) * si(i, j, k, 2) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 2) &
177  + w(i, j + 1, k, ivx) * sj(i, j, k, 2) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 2) &
178  + w(i, j, k + 1, ivx) * sk(i, j, k, 2) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 2)
179  uuz = w(i + 1, j, k, ivx) * si(i, j, k, 3) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 3) &
180  + w(i, j + 1, k, ivx) * sj(i, j, k, 3) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 3) &
181  + w(i, j, k + 1, ivx) * sk(i, j, k, 3) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 3)
182 
183  ! Idem for the gradient of v.
184 
185  vvx = w(i + 1, j, k, ivy) * si(i, j, k, 1) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 1) &
186  + w(i, j + 1, k, ivy) * sj(i, j, k, 1) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 1) &
187  + w(i, j, k + 1, ivy) * sk(i, j, k, 1) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 1)
188  vvy = w(i + 1, j, k, ivy) * si(i, j, k, 2) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 2) &
189  + w(i, j + 1, k, ivy) * sj(i, j, k, 2) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 2) &
190  + w(i, j, k + 1, ivy) * sk(i, j, k, 2) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 2)
191  vvz = w(i + 1, j, k, ivy) * si(i, j, k, 3) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 3) &
192  + w(i, j + 1, k, ivy) * sj(i, j, k, 3) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 3) &
193  + w(i, j, k + 1, ivy) * sk(i, j, k, 3) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 3)
194 
195  ! And for the gradient of w.
196 
197  wwx = w(i + 1, j, k, ivz) * si(i, j, k, 1) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 1) &
198  + w(i, j + 1, k, ivz) * sj(i, j, k, 1) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 1) &
199  + w(i, j, k + 1, ivz) * sk(i, j, k, 1) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 1)
200  wwy = w(i + 1, j, k, ivz) * si(i, j, k, 2) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 2) &
201  + w(i, j + 1, k, ivz) * sj(i, j, k, 2) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 2) &
202  + w(i, j, k + 1, ivz) * sk(i, j, k, 2) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 2)
203  wwz = w(i + 1, j, k, ivz) * si(i, j, k, 3) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 3) &
204  + w(i, j + 1, k, ivz) * sj(i, j, k, 3) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 3) &
205  + w(i, j, k + 1, ivz) * sk(i, j, k, 3) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 3)
206 
207  ! Compute the components of the stress tensor.
208  ! The combination of the current scaling of the velocity
209  ! gradients (2*vol) and the definition of the stress tensor,
210  ! leads to the factor 1/(4*vol).
211 
212  fact = fourth / vol(i, j, k)
213 
214  sxx = two * fact * uux
215  syy = two * fact * vvy
216  szz = two * fact * wwz
217 
218  sxy = fact * (uuy + vvx)
219  sxz = fact * (uuz + wwx)
220  syz = fact * (vvz + wwy)
221 
222  ! Compute 2/3 * divergence of velocity squared
223 
224  div2 = f23 * (sxx + syy + szz)**2
225 
226  ! Store the square of strain as the production term.
227 
228  scratch(i, j, k, iprod) = two * (two * (sxy**2 + sxz**2 + syz**2) &
229  + sxx**2 + syy**2 + szz**2) - div2
230 #ifdef TAPENADE_REVERSE
231  end do
232 #else
233  end do
234  end do
235  end do
236 #endif
237  end subroutine prodsmag2
238 
239  subroutine prodwmag2
240  !
241  ! prodWmag2 computes the term:
242  ! 2*oij*oij with oij=0.5*(duidxj - dujdxi).
243  ! This is equal to the magnitude squared of the vorticity.
244  ! It is assumed that the pointer vort, stored in turbMod, is
245  ! already set to the correct entry.
246  !
247  use constants
248  use blockpointers, only: nx, ny, nz, il, jl, kl, w, si, sj, sk, vol, sectionid, scratch
249  use flowvarrefstate, only: timeref
250  use section, only: sections
251  implicit none
252  !
253  ! Local variables.
254  !
255  integer :: i, j, k, ii
256 
257  real(kind=realtype) :: uuy, uuz, vvx, vvz, wwx, wwy
258  real(kind=realtype) :: fact, vortx, vorty, vortz
259  real(kind=realtype) :: omegax, omegay, omegaz
260 
261  ! Determine the non-dimensional wheel speed of this block.
262 
263  omegax = timeref * sections(sectionid)%rotRate(1)
264  omegay = timeref * sections(sectionid)%rotRate(2)
265  omegaz = timeref * sections(sectionid)%rotRate(3)
266 
267  ! Loop over the cell centers of the given block. It may be more
268  ! efficient to loop over the faces and to scatter the gradient,
269  ! but in that case the gradients for u, v and w must be stored.
270  ! In the current approach no extra memory is needed.
271 #ifdef TAPENADE_REVERSE
272  !$AD II-LOOP
273  do ii = 0, nx * ny * nz - 1
274  i = mod(ii, nx) + 2
275  j = mod(ii / nx, ny) + 2
276  k = ii / (nx * ny) + 2
277 #else
278  do k = 2, kl
279  do j = 2, jl
280  do i = 2, il
281 #endif
282 
283  ! Compute the necessary derivatives of u in the cell center.
284  ! Use is made of the fact that the surrounding normals sum up
285  ! to zero, such that the cell i,j,k does not give a
286  ! contribution. The gradient is scaled by a factor 2*vol.
287 
288  uuy = w(i + 1, j, k, ivx) * si(i, j, k, 2) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 2) &
289  + w(i, j + 1, k, ivx) * sj(i, j, k, 2) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 2) &
290  + w(i, j, k + 1, ivx) * sk(i, j, k, 2) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 2)
291  uuz = w(i + 1, j, k, ivx) * si(i, j, k, 3) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 3) &
292  + w(i, j + 1, k, ivx) * sj(i, j, k, 3) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 3) &
293  + w(i, j, k + 1, ivx) * sk(i, j, k, 3) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 3)
294 
295  ! Idem for the gradient of v.
296 
297  vvx = w(i + 1, j, k, ivy) * si(i, j, k, 1) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 1) &
298  + w(i, j + 1, k, ivy) * sj(i, j, k, 1) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 1) &
299  + w(i, j, k + 1, ivy) * sk(i, j, k, 1) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 1)
300  vvz = w(i + 1, j, k, ivy) * si(i, j, k, 3) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 3) &
301  + w(i, j + 1, k, ivy) * sj(i, j, k, 3) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 3) &
302  + w(i, j, k + 1, ivy) * sk(i, j, k, 3) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 3)
303 
304  ! And for the gradient of w.
305 
306  wwx = w(i + 1, j, k, ivz) * si(i, j, k, 1) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 1) &
307  + w(i, j + 1, k, ivz) * sj(i, j, k, 1) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 1) &
308  + w(i, j, k + 1, ivz) * sk(i, j, k, 1) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 1)
309  wwy = w(i + 1, j, k, ivz) * si(i, j, k, 2) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 2) &
310  + w(i, j + 1, k, ivz) * sj(i, j, k, 2) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 2) &
311  + w(i, j, k + 1, ivz) * sk(i, j, k, 2) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 2)
312 
313  ! Compute the three components of the vorticity vector.
314  ! Substract the part coming from the rotating frame.
315 
316  fact = half / vol(i, j, k)
317 
318  vortx = fact * (wwy - vvz) - two * omegax
319  vorty = fact * (uuz - wwx) - two * omegay
320  vortz = fact * (vvx - uuy) - two * omegaz
321 
322  ! Compute the magnitude squared of the vorticity.
323 
324  scratch(i, j, k, ivort) = vortx**2 + vorty**2 + vortz**2
325 #ifdef TAPENADE_REVERSE
326  end do
327 #else
328  end do
329  end do
330  end do
331 #endif
332  end subroutine prodwmag2
333  function sanuknowneddyratio(eddyRatio, nuLam)
334  !
335  ! saNuKnownEddyRatio computes the Spalart-Allmaras transport
336  ! variable nu for the given eddy viscosity ratio.
337  !
338  use constants
339  use paramturb
340  implicit none
341  !
342  ! Function type.
343  !
344  real(kind=realtype) :: sanuknowneddyratio
345  !
346  ! Function arguments.
347  !
348  real(kind=realtype), intent(in) :: eddyratio, nulam
349  !
350  ! Local variables.
351  !
352  real(kind=realtype) :: cv13, chi, chi2, chi3, chi4, f, df, dchi
353 
354  ! Take care of the exceptional cases.
355 
356  if (eddyratio <= zero) then
357  sanuknowneddyratio = zero
358  return
359  end if
360 
361  ! Set the value of cv1^3, which is the constant appearing in the
362  ! sa function fv1 to compute the eddy viscosity
363 
364  cv13 = rsacv1**3
365 
366  ! Determine the value of chi, which is given by the quartic
367  ! polynomial chi^4 - ratio*(chi^3 + cv1^3) = 0.
368  ! First determine the start value, depending on the eddyRatio.
369 
370  if (eddyratio < 1.e-4_realtype) then
371  chi = 0.5_realtype
372  else if (eddyratio < 1.0_realtype) then
373  chi = 5.0_realtype
374  else if (eddyratio < 10.0_realtype) then
375  chi = 10.0_realtype
376  else
377  chi = eddyratio
378  end if
379 
380  ! The actual newton algorithm.
381 
382  do
383  ! Compute the function value and the derivative.
384 
385  chi2 = chi * chi
386  chi3 = chi * chi2
387  chi4 = chi * chi3
388 
389  f = chi4 - eddyratio * (chi3 + cv13)
390  df = four * chi3 - three * eddyratio * chi2
391 
392  ! Compute the negative update and the new value of chi.
393 
394  dchi = f / df
395  chi = chi - dchi
396 
397  ! Condition to exit the loop.
398 
399  if (abs(dchi / chi) <= thresholdreal) exit
400  end do
401 
402  ! Chi is the ratio of the spalart allmaras transport variable and
403  ! the laminar viscosity. So multiply chi with the laminar viscosity
404  ! to obtain the correct value.
405 
406  sanuknowneddyratio = nulam * chi
407 
408  end function sanuknowneddyratio
409 
410  subroutine unsteadyturbterm(mAdv, nAdv, offset, qq)
411  !
412  ! unsteadyTurbTerm discretizes the time derivative of the
413  ! turbulence transport equations and add it to the residual.
414  ! As the time derivative is the same for all turbulence models,
415  ! this generic routine can be used; both the discretization of
416  ! the time derivative and its contribution to the central
417  ! jacobian are computed by this routine.
418  ! Only nAdv equations are treated, while the actual system has
419  ! size mAdv. The reason is that some equations for some
420  ! turbulence equations do not have a time derivative, e.g. the
421  ! f equation in the v2-f model. The argument offset indicates
422  ! the offset in the w vector where this subsystem starts. As a
423  ! consequence it is assumed that the indices of the current
424  ! subsystem are contiguous, e.g. if a 2*2 system is solved the
425  ! Last index in w is offset+1 and offset+2 respectively.
426  !
427  use blockpointers
428  use flowvarrefstate
429  use inputphysics
430  use inputtimespectral
431  use inputunsteady
432  use iteration
433  use section
434  use turbmod
435 
436  implicit none
437  !
438  ! Subroutine arguments.
439  !
440  integer(kind=intType), intent(in) :: mAdv, nAdv, offset
441 
442  real(kind=realtype), dimension(2:il, 2:jl, 2:kl, mAdv, mAdv), &
443  intent(inout) :: qq
444  !
445  ! Local variables.
446  !
447  integer(kind=intType) :: i, j, k, ii, jj, nn
448 
449  real(kind=realtype) :: oneoverdt, tmp
450 
451  ! Determine the equation mode.
452 
453  select case (equationmode)
454 
455  case (steady)
456 
457  ! Steady computation. No time derivative present.
458 
459  return
460 
461  !===============================================================
462 
463  case (unsteady)
464 
465  ! The time deritvative term depends on the integration
466  ! scheme used.
467 
468  select case (timeintegrationscheme)
469 
470  case (bdf)
471 
472  ! Backward difference formula is used as time
473  ! integration scheme.
474 
475  ! Store the inverse of the physical nonDimensional
476  ! time step a bit easier.
477 
478  oneoverdt = timeref / deltat
479 
480  ! Loop over the number of turbulent transport equations.
481 
482  nadvloopunsteady: do ii = 1, nadv
483 
484  ! Store the index of the current turbulent variable in jj.
485 
486  jj = ii + offset
487 
488  ! Loop over the owned cells of this block to compute the
489  ! time derivative.
490 
491  do k = 2, kl
492  do j = 2, jl
493  do i = 2, il
494 
495  ! Initialize tmp to the value of the current
496  ! level multiplied by the corresponding coefficient
497  ! in the time integration scheme.
498 
499  tmp = coeftime(0) * w(i, j, k, jj)
500 
501  ! Loop over the old time levels and add the
502  ! corresponding contribution to tmp.
503 
504  do nn = 1, noldlevels
505  tmp = tmp + coeftime(nn) * wold(nn, i, j, k, jj)
506  end do
507 
508  ! Update the residual. Note that in the turbulent
509  ! routines the residual is defined with an opposite
510  ! sign compared to the residual of the flow equations.
511  ! Therefore the time derivative must be substracted
512  ! from dvt.
513 
514  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - oneoverdt * tmp
515 
516  ! Update the central jacobian.
517 
518  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) &
519  + coeftime(0) * oneoverdt
520  end do
521  end do
522  end do
523 
524  end do nadvloopunsteady
525 
526  !===========================================================
527 
528  case (explicitrk)
529 
530  ! Explicit time integration scheme. The time derivative
531  ! is handled differently.
532 
533  return
534 
535  end select
536 
537  !===============================================================
538 
539  case (timespectral)
540 
541  ! Time spectral method.
542 
543  ! Loop over the number of turbulent transport equations.
544 
545  nadvloopspectral: do ii = 1, nadv
546 
547  ! Store the index of the current turbulent variable in jj.
548 
549  jj = ii + offset
550 
551  ! The time derivative has been computed earlier in
552  ! unsteadyTurbSpectral and stored in entry jj of scratch.
553  ! Substract this value for all owned cells. It must be
554  ! substracted, because in the turbulent routines the
555  ! residual is defined with an opposite sign compared to
556  ! the residual of the flow equations.
557  ! Also add a term to the diagonal matrix, which corresponds
558  ! to to the contribution of the highest frequency. This is
559  ! equivalent to an explicit treatment of the time derivative
560  ! and may need to be changed.
561 
562  tmp = ntimeintervalsspectral * pi * timeref &
563  / sections(sectionid)%timePeriod
564 
565  do k = 2, kl
566  do j = 2, jl
567  do i = 2, il
568  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - dw(i, j, k, jj)
569  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) + tmp
570  end do
571  end do
572  end do
573 
574  end do nadvloopspectral
575 
576  end select
577 
578  end subroutine unsteadyturbterm
579 
580  subroutine computeeddyviscosity(includeHalos)
581  !
582  ! computeEddyViscosity computes the eddy viscosity in the
583  ! owned cell centers of the given block. It is assumed that the
584  ! pointes already point to the correct block before entering
585  ! this subroutine.
586  !
587  use constants
588  use flowvarrefstate
589  use inputphysics
590  use iteration
591  use blockpointers
592  implicit none
593 
594  ! Input Parameter
595  logical, intent(in) :: includeHalos
596 
597  !
598  ! Local variables.
599  !
600  logical :: returnImmediately
601  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
602 
603  ! Check if an immediate return can be made.
604 
605  if (eddymodel) then
606  if ((currentlevel <= groundlevel)) then
607  returnimmediately = .false.
608  else
609  returnimmediately = .true.
610  end if
611  else
612  returnimmediately = .true.
613  end if
614 
615  if (returnimmediately) return
616 
617  ! Determine the turbulence model and call the appropriate
618  ! routine to compute the eddy viscosity.
619  if (includehalos) then
620  ibeg = 1
621  iend = ie
622  jbeg = 1
623  jend = je
624  kbeg = 1
625  kend = ke
626  else
627  ibeg = 2
628  iend = il
629  jbeg = 2
630  jend = jl
631  kbeg = 2
632  kend = kl
633  end if
634 
635  select case (turbmodel)
636 
637  case (spalartallmaras, spalartallmarasedwards)
638  call saeddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
639 #ifndef USE_TAPENADE
640 
641  case (v2f)
642  call vfeddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
643  case (komegawilcox, komegamodified)
644  call kweddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
645 
646  case (mentersst)
647  call ssteddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
648 
649  case (ktau)
650  call kteddyviscosity(ibeg, iend, jbeg, jend, kbeg, kend)
651 #endif
652  end select
653 
654  end subroutine computeeddyviscosity
655 
656  subroutine saeddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
657  !
658  ! saEddyViscosity computes the eddy-viscosity according to the
659  ! Spalart-Allmaras model for the block given in blockPointers.
660  ! This routine for both the original version as well as the
661  ! modified version according to Edwards.
662  !
663  use constants
664  use blockpointers
665  use constants
666  use paramturb
667  implicit none
668  ! Input variables
669  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
670 
671  !
672  ! Local variables.
673  !
674  integer(kind=intType) :: i, j, k, ii, iSize, jSize, kSize
675  real(kind=realtype) :: chi, chi3, fv1, rnusa, cv13
676 
677  ! Store the cv1^3; cv1 is a constant of the Spalart-Allmaras model.
678 
679  cv13 = rsacv1**3
680 
681  ! Loop over the cells of this block and compute the eddy viscosity.
682  ! Do not include halo's.
683 #ifdef TAPENADE_REVERSE
684  isize = (iend - ibeg) + 1
685  jsize = (jend - jbeg) + 1
686  ksize = (kend - kbeg) + 1
687 
688  !$AD II-LOOP
689  do ii = 0, isize * jsize * ksize - 1
690  i = mod(ii, isize) + ibeg
691  j = mod(ii / isize, jsize) + jbeg
692  k = ii / ((isize * jsize)) + kbeg
693 #else
694  do k = kbeg, kend
695  do j = jbeg, jend
696  do i = ibeg, iend
697 #endif
698  rnusa = w(i, j, k, itu1) * w(i, j, k, irho)
699  chi = rnusa / rlv(i, j, k)
700  chi3 = chi**3
701  fv1 = chi3 / (chi3 + cv13)
702  rev(i, j, k) = fv1 * rnusa
703 #ifdef TAPENADE_REVERSE
704  end do
705 #else
706  end do
707  end do
708  end do
709 #endif
710  end subroutine saeddyviscosity
711 
712  subroutine kweddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
713  !
714  ! kwEddyViscosity computes the eddy viscosity according to the
715  ! k-omega models (both the original Wilcox as well as the
716  ! modified version) for the block given in blockPointers.
717  !
718  use constants
719  use blockpointers
720  implicit none
721  ! Input variables
722  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
723  !
724  ! Local variables.
725  !
726  integer(kind=intType) :: i, j, k, ii, iSize, jSize, kSize
727 
728  ! Loop over the cells of this block and compute the eddy viscosity.
729  ! Do not include halo's.
730 #ifdef TAPENADE_REVERSE
731  isize = (iend - ibeg) + 1
732  jsize = (jend - jbeg) + 1
733  ksize = (kend - kbeg) + 1
734 
735  !$AD II-LOOP
736  do ii = 0, isize * jsize * ksize - 1
737  i = mod(ii, isize) + ibeg
738  j = mod(ii / isize, jsize) + jbeg
739  k = ii / ((isize * jsize)) + kbeg
740 #else
741  do k = kbeg, kend
742  do j = jbeg, jend
743  do i = ibeg, iend
744 #endif
745  rev(i, j, k) = abs(w(i, j, k, irho) * w(i, j, k, itu1) / w(i, j, k, itu2))
746 #ifdef TAPENADE_REVERSE
747  end do
748 #else
749  end do
750  end do
751  end do
752 #endif
753 
754  end subroutine kweddyviscosity
755 
756  subroutine ssteddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
757  !
758  ! SSTEddyViscosity computes the eddy viscosity according to
759  ! menter's SST variant of the k-omega turbulence model for the
760  ! block given in blockPointers.
761  !
762  use constants
763  use blockpointers
764  use paramturb
765  use turbmod
766  implicit none
767  ! Input variables
768  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
769  !
770  ! Local variables.
771  !
772  integer(kind=intType) :: i, j, k, ii, iSize, jSize, kSize
773  real(kind=realtype) :: t1, t2, arg2, f2, vortmag
774 
775  ! Compute the vorticity squared in the cell centers. The reason
776  ! for computing the vorticity squared is that a routine exists
777  ! for it; for the actual eddy viscosity computation the vorticity
778  ! itself is needed.
779 
780  call prodwmag2
781 
782  ! Loop over the cells of this block and compute the eddy viscosity.
783  ! Do not include halo's.
784 #ifdef TAPENADE_REVERSE
785  isize = (iend - ibeg) + 1
786  jsize = (jend - jbeg) + 1
787  ksize = (kend - kbeg) + 1
788 
789  !$AD II-LOOP
790  do ii = 0, isize * jsize * ksize - 1
791  i = mod(ii, isize) + ibeg
792  j = mod(ii / isize, jsize) + jbeg
793  k = ii / ((isize * jsize)) + kbeg
794 #else
795  do k = kbeg, kend
796  do j = jbeg, jend
797  do i = ibeg, iend
798 #endif
799  ! Compute the value of the function f2, which occurs in the
800  ! eddy-viscosity computation.
801 
802  t1 = two * sqrt(w(i, j, k, itu1)) &
803  / (0.09_realtype * w(i, j, k, itu2) * d2wall(i, j, k))
804  t2 = 500.0_realtype * rlv(i, j, k) &
805  / (w(i, j, k, irho) * w(i, j, k, itu2) * d2wall(i, j, k)**2)
806 
807  arg2 = max(t1, t2)
808  f2 = tanh(arg2**2)
809 
810  ! And compute the eddy viscosity.
811 
812  vortmag = sqrt(scratch(i, j, k, iprod))
813  rev(i, j, k) = w(i, j, k, irho) * rssta1 * w(i, j, k, itu1) &
814  / max(rssta1 * w(i, j, k, itu2), f2 * vortmag)
815 #ifdef TAPENADE_REVERSE
816  end do
817 #else
818  end do
819  end do
820  end do
821 #endif
822 
823  end subroutine ssteddyviscosity
824 
825  subroutine turbadvection(mAdv, nAdv, offset, qq)
826  !
827  ! turbAdvection discretizes the advection part of the turbulent
828  ! transport equations. As the advection part is the same for all
829  ! models, this generic routine can be used. Both the
830  ! discretization and the central jacobian are computed in this
831  ! subroutine. The former can either be 1st or 2nd order
832  ! accurate; the latter is always based on the 1st order upwind
833  ! discretization. When the discretization must be second order
834  ! accurate, the fully upwind (kappa = -1) scheme in combination
835  ! with the minmod limiter is used.
836  ! Only nAdv equations are treated, while the actual system has
837  ! size mAdv. The reason is that some equations for some
838  ! turbulence equations do not have an advection part, e.g. the
839  ! f equation in the v2-f model. The argument offset indicates
840  ! the offset in the w vector where this subsystem starts. As a
841  ! consequence it is assumed that the indices of the current
842  ! subsystem are contiguous, e.g. if a 2*2 system is solved the
843  ! Last index in w is offset+1 and offset+2 respectively.
844  !
845  use constants
846  use blockpointers, only: nx, ny, nz, il, jl, kl, vol, sfacei, sfacej, sfacek, &
847  w, si, sj, sk, addgridvelocities, bmti1, bmti2, bmtj1, bmtj2, &
848  bmtk1, bmtk2, scratch
849  use inputdiscretization, only: orderturb
850  use iteration, only: groundlevel
851  use turbmod, only: secondord
852  implicit none
853  !
854  ! Subroutine arguments.
855  !
856  integer(kind=intType), intent(in) :: nAdv, mAdv, offset
857 
858  real(kind=realtype), dimension(2:il, 2:jl, 2:kl, mAdv, mAdv), &
859  intent(inout) :: qq
860  !
861  ! Local variables.
862  !
863  integer(kind=intType) :: i, j, k, ii, jj, kk, iii
864 
865  real(kind=realtype) :: qs, voli, xa, ya, za
866  real(kind=realtype) :: uu, dwt, dwtm1, dwtp1, dwti, dwtj, dwtk
867 
868  real(kind=realtype), dimension(mAdv) :: impl
869 
870  ! Determine whether or not a second order discretization for the
871  ! advective terms must be used.
872  secondord = .false.
873  if (groundlevel == 1_inttype .and. &
874  orderturb == secondorder) secondord = .true.
875 
876  ! Initialize the grid velocity to zero. This value will be used
877  ! if the block is not moving.
878  continue
879  !$AD CHECKPOINT-START
880  qs = zero
881  !
882  ! Upwind discretization of the convective term in k (zeta)
883  ! direction. Either the 1st order upwind or the second order
884  ! fully upwind interpolation scheme, kappa = -1, is used in
885  ! combination with the minmod limiter.
886  ! The possible grid velocity must be taken into account.
887  !
888 #ifdef TAPENADE_REVERSE
889  !$AD II-LOOP
890  do iii = 0, nx * ny * nz - 1
891  i = mod(iii, nx) + 2
892  j = mod(iii / nx, ny) + 2
893  k = iii / (nx * ny) + 2
894 #else
895  do k = 2, kl
896  do j = 2, jl
897  do i = 2, il
898 #endif
899  ! Compute the grid velocity if present.
900  ! It is taken as the average of k and k-1,
901 
902  voli = half / vol(i, j, k)
903  if (addgridvelocities) &
904  qs = (sfacek(i, j, k) + sfacek(i, j, k - 1)) * voli
905 
906  ! Compute the normal velocity, where the normal direction
907  ! is taken as the average of faces k and k-1.
908 
909  xa = (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
910  ya = (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
911  za = (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
912 
913  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
914  ! This term has unit: velocity/length
915 
916  ! Determine the situation we are having here, i.e. positive
917  ! or negative normal velocity.
918 
919  velkdir: if (uu > zero) then
920 
921  ! Velocity has a component in positive k-direction.
922  ! Loop over the number of advection equations.
923  !$AD II-LOOP
924  do ii = 1, nadv
925 
926  ! Set the value of jj such that it corresponds to the
927  ! turbulent entry in w.
928 
929  jj = ii + offset
930 
931  ! Check whether a first or a second order discretization
932  ! must be used.
933 
934  if (secondord) then
935 
936  ! Second order; store the three differences for the
937  ! discretization of the derivative in k-direction.
938 
939  dwtm1 = w(i, j, k - 1, jj) - w(i, j, k - 2, jj)
940  dwt = w(i, j, k, jj) - w(i, j, k - 1, jj)
941  dwtp1 = w(i, j, k + 1, jj) - w(i, j, k, jj)
942 
943  ! Construct the derivative in this cell center. This
944  ! is the first order upwind derivative with two
945  ! nonlinear corrections.
946 
947  dwtk = dwt
948 
949  if (dwt * dwtp1 > zero) then
950  if (abs(dwt) < abs(dwtp1)) then
951  dwtk = dwtk + half * dwt
952  else
953  dwtk = dwtk + half * dwtp1
954  end if
955  end if
956 
957  if (dwt * dwtm1 > zero) then
958  if (abs(dwt) < abs(dwtm1)) then
959  dwtk = dwtk - half * dwt
960  else
961  dwtk = dwtk - half * dwtm1
962  end if
963  end if
964 
965  else
966 
967  ! 1st order upwind scheme.
968 
969  dwtk = w(i, j, k, jj) - w(i, j, k - 1, jj)
970 
971  end if
972 
973  ! Update the residual. The convective term must be
974  ! substracted, because it appears on the other side of
975  ! the equation as the source and viscous terms.
976  ! uu*dwtk = (V.dot.face_normal)*delta(nuTilde)/delta(x)
977 
978  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - uu * dwtk
979 #ifndef USE_TAPENADE
980  ! Update the central jacobian. First the term which is
981  ! always present, i.e. uu.
982 
983  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) + uu
984 
985  ! For boundary cells k == 2, the implicit treatment must
986  ! be taken into account. Note that the implicit part
987  ! is only based on the 1st order discretization.
988  ! To improve stability the diagonal term is only taken
989  ! into account when it improves stability, i.e. when
990  ! it is positive.
991 
992  if (k == 2) then
993  do kk = 1, madv
994  impl(kk) = bmtk1(i, j, jj, kk + offset)
995  end do
996 
997  impl(ii) = max(impl(ii), zero)
998 
999  do kk = 1, madv
1000  qq(i, j, k, ii, kk) = qq(i, j, k, ii, kk) + uu * impl(kk)
1001  end do
1002  end if
1003 #endif
1004 
1005  end do
1006 
1007  else velkdir
1008 
1009  ! Velocity has a component in negative k-direction.
1010  ! Loop over the number of advection equations
1011  !$AD II-LOOP
1012  do ii = 1, nadv
1013 
1014  ! Set the value of jj such that it corresponds to the
1015  ! turbulent entry in w.
1016 
1017  jj = ii + offset
1018 
1019  ! Check whether a first or a second order discretization
1020  ! must be used.
1021 
1022  if (secondord) then
1023 
1024  ! Store the three differences for the discretization of
1025  ! the derivative in k-direction.
1026 
1027  dwtm1 = w(i, j, k, jj) - w(i, j, k - 1, jj)
1028  dwt = w(i, j, k + 1, jj) - w(i, j, k, jj)
1029  dwtp1 = w(i, j, k + 2, jj) - w(i, j, k + 1, jj)
1030 
1031  ! Construct the derivative in this cell center. This is
1032  ! the first order upwind derivative with two nonlinear
1033  ! corrections.
1034 
1035  dwtk = dwt
1036 
1037  if (dwt * dwtp1 > zero) then
1038  if (abs(dwt) < abs(dwtp1)) then
1039  dwtk = dwtk - half * dwt
1040  else
1041  dwtk = dwtk - half * dwtp1
1042  end if
1043  end if
1044 
1045  if (dwt * dwtm1 > zero) then
1046  if (abs(dwt) < abs(dwtm1)) then
1047  dwtk = dwtk + half * dwt
1048  else
1049  dwtk = dwtk + half * dwtm1
1050  end if
1051  end if
1052 
1053  else
1054 
1055  ! 1st order upwind scheme.
1056 
1057  dwtk = w(i, j, k + 1, jj) - w(i, j, k, jj)
1058 
1059  end if
1060 
1061  ! Update the residual. The convective term must be
1062  ! substracted, because it appears on the other side
1063  ! of the equation as the source and viscous terms.
1064 
1065  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - uu * dwtk
1066 
1067  ! Update the central jacobian. First the term which is
1068  ! always present, i.e. -uu.
1069 #ifndef USE_TAPENADE
1070  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) - uu
1071 
1072  ! For boundary cells k == kl, the implicit treatment must
1073  ! be taken into account. Note that the implicit part
1074  ! is only based on the 1st order discretization.
1075  ! To improve stability the diagonal term is only taken
1076  ! into account when it improves stability, i.e. when
1077  ! it is positive.
1078 
1079  if (k == kl) then
1080  do kk = 1, madv
1081  impl(kk) = bmtk2(i, j, jj, kk + offset)
1082  end do
1083 
1084  impl(ii) = max(impl(ii), zero)
1085 
1086  do kk = 1, madv
1087  qq(i, j, k, ii, kk) = qq(i, j, k, ii, kk) - uu * impl(kk)
1088  end do
1089  end if
1090 #endif
1091  end do
1092 
1093  end if velkdir
1094 #ifdef TAPENADE_REVERSE
1095  end do
1096 #else
1097  end do
1098  end do
1099  end do
1100 #endif
1101  continue
1102  !$AD CHECKPOINT-END
1103  !
1104  ! Upwind discretization of the convective term in j (eta)
1105  ! direction. Either the 1st order upwind or the second order
1106  ! fully upwind interpolation scheme, kappa = -1, is used in
1107  ! combination with the minmod limiter.
1108  ! The possible grid velocity must be taken into account.
1109  !
1110  continue
1111  !$AD CHECKPOINT-START
1112  qs = zero
1113 #ifdef TAPENADE_REVERSE
1114  !$AD II-LOOP
1115  do iii = 0, nx * ny * nz - 1
1116  i = mod(iii, nx) + 2
1117  j = mod(iii / nx, ny) + 2
1118  k = iii / (nx * ny) + 2
1119 #else
1120  do k = 2, kl
1121  do j = 2, jl
1122  do i = 2, il
1123 #endif
1124 
1125  ! Compute the grid velocity if present.
1126  ! It is taken as the average of j and j-1,
1127 
1128  voli = half / vol(i, j, k)
1129  if (addgridvelocities) &
1130  qs = (sfacej(i, j, k) + sfacej(i, j - 1, k)) * voli
1131 
1132  ! Compute the normal velocity, where the normal direction
1133  ! is taken as the average of faces j and j-1.
1134 
1135  xa = (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
1136  ya = (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
1137  za = (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
1138 
1139  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1140 
1141  ! Determine the situation we are having here, i.e. positive
1142  ! or negative normal velocity.
1143 
1144  veljdir: if (uu > zero) then
1145 
1146  ! Velocity has a component in positive j-direction.
1147  ! Loop over the number of advection equations.
1148  !$AD II-LOOP
1149  do ii = 1, nadv
1150 
1151  ! Set the value of jj such that it corresponds to the
1152  ! turbulent entry in w.
1153 
1154  jj = ii + offset
1155 
1156  ! Check whether a first or a second order discretization
1157  ! must be used.
1158 
1159  if (secondord) then
1160 
1161  ! Second order; store the three differences for the
1162  ! discretization of the derivative in j-direction.
1163 
1164  dwtm1 = w(i, j - 1, k, jj) - w(i, j - 2, k, jj)
1165  dwt = w(i, j, k, jj) - w(i, j - 1, k, jj)
1166  dwtp1 = w(i, j + 1, k, jj) - w(i, j, k, jj)
1167 
1168  ! Construct the derivative in this cell center. This is
1169  ! the first order upwind derivative with two nonlinear
1170  ! corrections.
1171 
1172  dwtj = dwt
1173 
1174  if (dwt * dwtp1 > zero) then
1175  if (abs(dwt) < abs(dwtp1)) then
1176  dwtj = dwtj + half * dwt
1177  else
1178  dwtj = dwtj + half * dwtp1
1179  end if
1180  end if
1181 
1182  if (dwt * dwtm1 > zero) then
1183  if (abs(dwt) < abs(dwtm1)) then
1184  dwtj = dwtj - half * dwt
1185  else
1186  dwtj = dwtj - half * dwtm1
1187  end if
1188  end if
1189 
1190  else
1191 
1192  ! 1st order upwind scheme.
1193 
1194  dwtj = w(i, j, k, jj) - w(i, j - 1, k, jj)
1195 
1196  end if
1197 
1198  ! Update the residual. The convective term must be
1199  ! substracted, because it appears on the other side of
1200  ! the equation as the source and viscous terms.
1201 
1202  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - uu * dwtj
1203 
1204  ! Update the central jacobian. First the term which is
1205  ! always present, i.e. uu.
1206 #ifndef USE_TAPENADE
1207  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) + uu
1208 
1209  ! For boundary cells j == 2, the implicit treatment must
1210  ! be taken into account. Note that the implicit part
1211  ! is only based on the 1st order discretization.
1212  ! To improve stability the diagonal term is only taken
1213  ! into account when it improves stability, i.e. when
1214  ! it is positive.
1215 
1216  if (j == 2) then
1217  do kk = 1, madv
1218  impl(kk) = bmtj1(i, k, jj, kk + offset)
1219  end do
1220 
1221  impl(ii) = max(impl(ii), zero)
1222 
1223  do kk = 1, madv
1224  qq(i, j, k, ii, kk) = qq(i, j, k, ii, kk) + uu * impl(kk)
1225  end do
1226  end if
1227 #endif
1228  end do
1229 
1230  else veljdir
1231 
1232  ! Velocity has a component in negative j-direction.
1233  ! Loop over the number of advection equations.
1234  !$AD II-LOOP
1235  do ii = 1, nadv
1236 
1237  ! Set the value of jj such that it corresponds to the
1238  ! turbulent entry in w.
1239 
1240  jj = ii + offset
1241 
1242  ! Check whether a first or a second order discretization
1243  ! must be used.
1244 
1245  if (secondord) then
1246 
1247  ! Store the three differences for the discretization of
1248  ! the derivative in j-direction.
1249 
1250  dwtm1 = w(i, j, k, jj) - w(i, j - 1, k, jj)
1251  dwt = w(i, j + 1, k, jj) - w(i, j, k, jj)
1252  dwtp1 = w(i, j + 2, k, jj) - w(i, j + 1, k, jj)
1253 
1254  ! Construct the derivative in this cell center. This is
1255  ! the first order upwind derivative with two nonlinear
1256  ! corrections.
1257 
1258  dwtj = dwt
1259 
1260  if (dwt * dwtp1 > zero) then
1261  if (abs(dwt) < abs(dwtp1)) then
1262  dwtj = dwtj - half * dwt
1263  else
1264  dwtj = dwtj - half * dwtp1
1265  end if
1266  end if
1267 
1268  if (dwt * dwtm1 > zero) then
1269  if (abs(dwt) < abs(dwtm1)) then
1270  dwtj = dwtj + half * dwt
1271  else
1272  dwtj = dwtj + half * dwtm1
1273  end if
1274  end if
1275 
1276  else
1277 
1278  ! 1st order upwind scheme.
1279 
1280  dwtj = w(i, j + 1, k, jj) - w(i, j, k, jj)
1281 
1282  end if
1283 
1284  ! Update the residual. The convective term must be
1285  ! substracted, because it appears on the other side
1286  ! of the equation as the source and viscous terms.
1287 
1288  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - uu * dwtj
1289 
1290  ! Update the central jacobian. First the term which is
1291  ! always present, i.e. -uu.
1292 #ifndef USE_TAPENADE
1293  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) - uu
1294 
1295  ! For boundary cells j == jl, the implicit treatment must
1296  ! be taken into account. Note that the implicit part
1297  ! is only based on the 1st order discretization.
1298  ! To improve stability the diagonal term is only taken
1299  ! into account when it improves stability, i.e. when
1300  ! it is positive.
1301 
1302  if (j == jl) then
1303  do kk = 1, madv
1304  impl(kk) = bmtj2(i, k, jj, kk + offset)
1305  end do
1306 
1307  impl(ii) = max(impl(ii), zero)
1308 
1309  do kk = 1, madv
1310  qq(i, j, k, ii, kk) = qq(i, j, k, ii, kk) - uu * impl(kk)
1311  end do
1312  end if
1313 #endif
1314  end do
1315 
1316  end if veljdir
1317 #ifdef TAPENADE_REVERSE
1318  end do
1319 #else
1320  end do
1321  end do
1322  end do
1323 #endif
1324  continue
1325  !$AD CHECKPOINT-END
1326  !
1327  ! Upwind discretization of the convective term in i (xi)
1328  ! direction. Either the 1st order upwind or the second order
1329  ! fully upwind interpolation scheme, kappa = -1, is used in
1330  ! combination with the minmod limiter.
1331  ! The possible grid velocity must be taken into account.
1332  !
1333  continue
1334  !$AD CHECKPOINT-START
1335  qs = zero
1336 #ifdef TAPENADE_REVERSE
1337  !$AD II-LOOP
1338  do iii = 0, nx * ny * nz - 1
1339  i = mod(iii, nx) + 2
1340  j = mod(iii / nx, ny) + 2
1341  k = iii / (nx * ny) + 2
1342 #else
1343  do k = 2, kl
1344  do j = 2, jl
1345  do i = 2, il
1346 #endif
1347 
1348  ! Compute the grid velocity if present.
1349  ! It is taken as the average of i and i-1,
1350 
1351  voli = half / vol(i, j, k)
1352  if (addgridvelocities) &
1353  qs = (sfacei(i, j, k) + sfacei(i - 1, j, k)) * voli
1354 
1355  ! Compute the normal velocity, where the normal direction
1356  ! is taken as the average of faces i and i-1.
1357 
1358  xa = (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
1359  ya = (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
1360  za = (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
1361 
1362  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1363 
1364  ! Determine the situation we are having here, i.e. positive
1365  ! or negative normal velocity.
1366 
1367  velidir: if (uu > zero) then
1368 
1369  ! Velocity has a component in positive i-direction.
1370  ! Loop over the number of advection equations.
1371  !$AD II-LOOP
1372  do ii = 1, nadv
1373 
1374  ! Set the value of jj such that it corresponds to the
1375  ! turbulent entry in w.
1376 
1377  jj = ii + offset
1378 
1379  ! Check whether a first or a second order discretization
1380  ! must be used.
1381 
1382  if (secondord) then
1383 
1384  ! Second order; store the three differences for the
1385  ! discretization of the derivative in i-direction.
1386 
1387  dwtm1 = w(i - 1, j, k, jj) - w(i - 2, j, k, jj)
1388  dwt = w(i, j, k, jj) - w(i - 1, j, k, jj)
1389  dwtp1 = w(i + 1, j, k, jj) - w(i, j, k, jj)
1390 
1391  ! Construct the derivative in this cell center. This is
1392  ! the first order upwind derivative with two nonlinear
1393  ! corrections.
1394 
1395  dwti = dwt
1396 
1397  if (dwt * dwtp1 > zero) then
1398  if (abs(dwt) < abs(dwtp1)) then
1399  dwti = dwti + half * dwt
1400  else
1401  dwti = dwti + half * dwtp1
1402  end if
1403  end if
1404 
1405  if (dwt * dwtm1 > zero) then
1406  if (abs(dwt) < abs(dwtm1)) then
1407  dwti = dwti - half * dwt
1408  else
1409  dwti = dwti - half * dwtm1
1410  end if
1411  end if
1412 
1413  else
1414 
1415  ! 1st order upwind scheme.
1416 
1417  dwti = w(i, j, k, jj) - w(i - 1, j, k, jj)
1418 
1419  end if
1420 
1421  ! Update the residual. The convective term must be
1422  ! substracted, because it appears on the other side of
1423  ! the equation as the source and viscous terms.
1424 
1425  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - uu * dwti
1426 
1427  ! Update the central jacobian. First the term which is
1428  ! always present, i.e. uu.
1429 #ifndef USE_TAPENADE
1430  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) + uu
1431 
1432  ! For boundary cells i == 2, the implicit treatment must
1433  ! be taken into account. Note that the implicit part
1434  ! is only based on the 1st order discretization.
1435  ! To improve stability the diagonal term is only taken
1436  ! into account when it improves stability, i.e. when
1437  ! it is positive.
1438 
1439  if (i == 2) then
1440  do kk = 1, madv
1441  impl(kk) = bmti1(j, k, jj, kk + offset)
1442  end do
1443 
1444  impl(ii) = max(impl(ii), zero)
1445 
1446  do kk = 1, madv
1447  qq(i, j, k, ii, kk) = qq(i, j, k, ii, kk) + uu * impl(kk)
1448  end do
1449  end if
1450 #endif
1451  end do
1452 
1453  else velidir
1454 
1455  ! Velocity has a component in negative i-direction.
1456  ! Loop over the number of advection equations.
1457  !$AD II-LOOP
1458  do ii = 1, nadv
1459 
1460  ! Set the value of jj such that it corresponds to the
1461  ! turbulent entry in w.
1462 
1463  jj = ii + offset
1464 
1465  ! Check whether a first or a second order discretization
1466  ! must be used.
1467 
1468  if (secondord) then
1469 
1470  ! Second order; store the three differences for the
1471  ! discretization of the derivative in i-direction.
1472 
1473  dwtm1 = w(i, j, k, jj) - w(i - 1, j, k, jj)
1474  dwt = w(i + 1, j, k, jj) - w(i, j, k, jj)
1475  dwtp1 = w(i + 2, j, k, jj) - w(i + 1, j, k, jj)
1476 
1477  ! Construct the derivative in this cell center. This is
1478  ! the first order upwind derivative with two nonlinear
1479  ! corrections.
1480 
1481  dwti = dwt
1482 
1483  if (dwt * dwtp1 > zero) then
1484  if (abs(dwt) < abs(dwtp1)) then
1485  dwti = dwti - half * dwt
1486  else
1487  dwti = dwti - half * dwtp1
1488  end if
1489  end if
1490 
1491  if (dwt * dwtm1 > zero) then
1492  if (abs(dwt) < abs(dwtm1)) then
1493  dwti = dwti + half * dwt
1494  else
1495  dwti = dwti + half * dwtm1
1496  end if
1497  end if
1498 
1499  else
1500 
1501  ! 1st order upwind scheme.
1502 
1503  dwti = w(i + 1, j, k, jj) - w(i, j, k, jj)
1504 
1505  end if
1506 
1507  ! Update the residual. The convective term must be
1508  ! substracted, because it appears on the other side
1509  ! of the equation as the source and viscous terms.
1510 
1511  scratch(i, j, k, idvt + ii - 1) = scratch(i, j, k, idvt + ii - 1) - uu * dwti
1512 
1513  ! Update the central jacobian. First the term which is
1514  ! always present, i.e. -uu.
1515 #ifndef USE_TAPENADE
1516  qq(i, j, k, ii, ii) = qq(i, j, k, ii, ii) - uu
1517 
1518  ! For boundary cells i == il, the implicit treatment must
1519  ! be taken into account. Note that the implicit part
1520  ! is only based on the 1st order discretization.
1521  ! To improve stability the diagonal term is only taken
1522  ! into account when it improves stability, i.e. when
1523  ! it is positive.
1524 
1525  if (i == il) then
1526  do kk = 1, madv
1527  impl(kk) = bmti2(j, k, jj, kk + offset)
1528  end do
1529 
1530  impl(ii) = max(impl(ii), zero)
1531 
1532  do kk = 1, madv
1533  qq(i, j, k, ii, kk) = qq(i, j, k, ii, kk) - uu * impl(kk)
1534  end do
1535  end if
1536 #endif
1537  end do
1538 
1539  end if velidir
1540 #ifdef TAPENADE_REVERSE
1541  end do
1542 #else
1543  end do
1544  end do
1545  end do
1546 #endif
1547 
1548  !$AD CHECKPOINT-END
1549  continue
1550  end subroutine turbadvection
1551 
1552  ! ----------------------------------------------------------------------
1553  ! |
1554  ! No Tapenade Routine below this line |
1555  ! |
1556  ! ----------------------------------------------------------------------
1557 
1558 #ifndef USE_TAPENADE
1559 
1560  subroutine tdia3(nb, ne, l, c, u, r)
1561  !
1562  ! tdia3 solves the tridiagonal linear system (l+c+u) v = r,
1563  ! where l is the lower, c the central and u the upper diagonal.
1564  ! Every entry in the matrix is a 2x2 block matrix, i.e.
1565  ! x x x 0 0 0 ........ = c(nb) u(nb)
1566  ! x x 0 x 0 0 ........
1567  ! x 0 x x x 0 ........ = l(i) c(i) u(i)
1568  ! 0 x x x 0 x ........
1569  ! ........ x 0 x x = l(ne) c(ne)
1570  ! ........ 0 x x x
1571  ! With c = x x u,l = x 0
1572  ! x x 0 x
1573  !
1574  use constants
1575  implicit none
1576  !
1577  ! Subroutine arguments.
1578  !
1579  integer(kind=intType), intent(in) :: nb, ne
1580 
1581  real(kind=realtype), dimension(2, nb:ne), intent(inout) :: l, u, r
1582  real(kind=realtype), dimension(2, 2, nb:ne), intent(inout) :: c
1583  !
1584  ! Local variables.
1585  !
1586  integer(kind=intType) :: n
1587  real(kind=realtype) :: deti, f11, f12, f21, f22, r1
1588 
1589  ! Perform the backward sweep to eliMinate the upper diagonal uu.
1590  ! f = u(n)*c^-1(n+1),
1591  ! c'(n) = c(n) - f*l(n+1)
1592  ! r'(n) = r(n) - f*r(n+1)
1593 
1594  do n = ne - 1, nb, -1
1595 
1596  deti = one / (c(1, 1, n + 1) * c(2, 2, n + 1) - c(1, 2, n + 1) * c(2, 1, n + 1))
1597  f11 = u(1, n) * c(2, 2, n + 1) * deti
1598  f12 = -u(1, n) * c(1, 2, n + 1) * deti
1599  f21 = -u(2, n) * c(2, 1, n + 1) * deti
1600  f22 = u(2, n) * c(1, 1, n + 1) * deti
1601 
1602  c(1, 1, n) = c(1, 1, n) - f11 * l(1, n + 1)
1603  c(1, 2, n) = c(1, 2, n) - f12 * l(2, n + 1)
1604  c(2, 1, n) = c(2, 1, n) - f21 * l(1, n + 1)
1605  c(2, 2, n) = c(2, 2, n) - f22 * l(2, n + 1)
1606 
1607  r(1, n) = r(1, n) - f11 * r(1, n + 1) - f12 * r(2, n + 1)
1608  r(2, n) = r(2, n) - f21 * r(1, n + 1) - f22 * r(2, n + 1)
1609 
1610  end do
1611 
1612  ! The matrix is now in low block bi-diagonal form and can thus
1613  ! be solved be a forward sweep. The solution is stored in r.
1614 
1615  deti = one / (c(1, 1, nb) * c(2, 2, nb) - c(1, 2, nb) * c(2, 1, nb))
1616  r1 = r(1, nb)
1617  r(1, nb) = deti * (c(2, 2, nb) * r1 - c(1, 2, nb) * r(2, nb))
1618  r(2, nb) = -deti * (c(2, 1, nb) * r1 - c(1, 1, nb) * r(2, nb))
1619 
1620  do n = nb + 1, ne
1621 
1622  r(1, n) = r(1, n) - l(1, n) * r(1, n - 1)
1623  r(2, n) = r(2, n) - l(2, n) * r(2, n - 1)
1624 
1625  deti = one / (c(1, 1, n) * c(2, 2, n) - c(1, 2, n) * c(2, 1, n))
1626  r1 = r(1, n)
1627  r(1, n) = deti * (c(2, 2, n) * r1 - c(1, 2, n) * r(2, n))
1628  r(2, n) = -deti * (c(2, 1, n) * r1 - c(1, 1, n) * r(2, n))
1629 
1630  end do
1631 
1632  end subroutine tdia3
1633 
1634  subroutine vfeddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
1635  !
1636  ! vfEddyViscosity computes the eddy-viscosity according to the
1637  ! v2f turbulence model for the block given in blockPointers.
1638  ! This routine is for both the n=1 and n=6 version.
1639  !
1640  use constants
1641  use blockpointers
1642  use constants
1643  use paramturb
1644  use turbmod
1645  use inputphysics
1646  use turbcurvefits, only: curvetupyp
1647  implicit none
1648  ! Input variables
1649  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
1650  !
1651  ! Local variables.
1652  !
1653  integer(kind=intType) :: i, j, k, nn, ii, iSize, jSize, kSize
1654  real(kind=realtype) :: tke, tep, tkea, tepa, tepl, tv2, tv2a
1655  real(kind=realtype) :: yp, utau
1656 
1657  real(kind=realtype), dimension(itu1:itu5) :: tup
1658  real(kind=realtype), dimension(:, :, :), pointer :: ww
1659  real(kind=realtype), dimension(:, :), pointer :: rrlv, rrev
1660  real(kind=realtype), dimension(:, :), pointer :: dd2wall
1661 
1662  ! Compute time and length scale
1663 
1664  call vfscale
1665 
1666  ! Loop over the cells of this block and compute the eddy viscosity.
1667  ! Do not include halo's.
1668 #ifdef TAPENADE_REVERSE
1669  isize = (iend - ibeg) + 1
1670  jsize = (jend - jbeg) + 1
1671  ksize = (kend - kbeg) + 1
1672 
1673  !$AD II-LOOP
1674  do ii = 0, isize * jsize * ksize - 1
1675  i = mod(ii, isize) + ibeg
1676  j = mod(ii / isize, jsize) + jbeg
1677  k = ii / ((isize * jsize)) + kbeg
1678 #else
1679  do k = kbeg, kend
1680  do j = jbeg, jend
1681  do i = ibeg, iend
1682 #endif
1683 
1684  tke = w(i, j, k, itu1)
1685  tep = w(i, j, k, itu2)
1686  tv2 = w(i, j, k, itu3)
1687  tkea = abs(tke)
1688  tepa = abs(tep)
1689  tv2a = abs(tv2)
1690  tepl = max(tepa, rvflimite)
1691 
1692  rev(i, j, k) = rvfcmu * w(i, j, k, irho) * tv2a / tepl * sct(i, j, k)
1693 #ifdef TAPENADE_REVERSE
1694  end do
1695 #else
1696  end do
1697  end do
1698  end do
1699 #endif
1700 
1701  ! Modify the rhs of the 1st internal cell, if wall functions
1702  ! are used; their value is determined by the table.
1703 
1704  testwallfunctions: if (wallfunctions) then
1705 
1706  bocos: do nn = 1, nviscbocos
1707 
1708  ! Determine the block face on which the subface is located
1709  ! and set some variables. As flag points to the entire array
1710  ! flagI2, etc., its starting indices are the starting indices
1711  ! of its target and not 1.
1712 
1713  select case (bcfaceid(nn))
1714  case (imin)
1715  ww => w(2, 1:, 1:, 1:); rrlv => rlv(2, 1:, 1:)
1716  dd2wall => d2wall(2, :, :); rrev => rev(2, 1:, 1:)
1717 
1718  case (imax)
1719  ww => w(il, 1:, 1:, 1:); rrlv => rlv(il, 1:, 1:)
1720  dd2wall => d2wall(il, :, :); rrev => rev(il, 1:, 1:)
1721 
1722  case (jmin)
1723  ww => w(1:, 2, 1:, 1:); rrlv => rlv(1:, 2, 1:)
1724  dd2wall => d2wall(:, 2, :); rrev => rev(1:, 2, 1:)
1725 
1726  case (jmax)
1727  ww => w(1:, jl, 1:, 1:); rrlv => rlv(1:, jl, 1:)
1728  dd2wall => d2wall(:, jl, :); rrev => rev(1:, jl, 1:)
1729 
1730  case (kmin)
1731  ww => w(1:, 1:, 2, 1:); rrlv => rlv(1:, 1:, 2)
1732  dd2wall => d2wall(:, :, 2); rrev => rev(1:, 1:, 2)
1733 
1734  case (kmax)
1735  ww => w(1:, 1:, kl, 1:); rrlv => rlv(1:, 1:, kl)
1736  dd2wall => d2wall(:, :, kl); rrev => rev(1:, 1:, kl)
1737 
1738  end select
1739 
1740  ! Loop over the owned faces of this subface. Therefore the
1741  ! nodal range of bcData must be used. The offset of +1 is
1742  ! present, because the starting index of the cell range is
1743  ! 1 larger than the starting index of the nodal range.
1744 
1745  do j = (bcdata(nn)%jnBeg + 1), bcdata(nn)%jnEnd
1746  do i = (bcdata(nn)%inBeg + 1), bcdata(nn)%inEnd
1747 
1748  ! Enforce k and epsilon in the 1st internal cell from
1749  ! the wall function table. There is an offset of -1 in
1750  ! the wall distance. Note that the offset compared to
1751  ! the current value must be stored. Also note that the
1752  ! curve fits contain the non-dimensional values.
1753 
1754  utau = viscsubface(nn)%utau(i, j)
1755  yp = ww(i, j, irho) * dd2wall(i - 1, j - 1) * utau / rrlv(i, j)
1756 
1757  call curvetupyp(tup(itu5:itu5), yp, itu5, itu5)
1758  rrev(i, j) = tup(itu5) * rrlv(i, j)
1759  end do
1760  end do
1761 
1762  end do bocos
1763  end if testwallfunctions
1764 
1765  end subroutine vfeddyviscosity
1766  subroutine vfscale
1767  !
1768  ! time and length scale definition for v2f turbulence model. The
1769  ! upper bound can be switched on by setting rvfB to .true.
1770  ! The strain squared is defined as: strain2 = 2 sij sij
1771  !
1772  use constants
1773  use blockpointers
1774  use inputphysics
1775  use paramturb
1776  use turbmod
1777  implicit none
1778  !
1779  ! Local variables.
1780  !
1781  integer(kind=intType) :: i, j, k
1782 
1783  real(kind=realtype) :: sqrt3
1784  real(kind=realtype) :: tkea, tepa, tv2a, supi, rn2
1785  real(kind=realtype) :: rsct, rscl2, rnu, rstrain
1786 
1787  ! Some constants in the model.
1788 
1789  sqrt3 = sqrt(three)
1790 
1791  ! Set the pointer for dvt in dw, such that the code is more
1792  ! readable. Also set the pointers for the production term,
1793  ! vorticity, strain and the time and lenght scale of v2f.
1794 
1795  dvt => scratch(1:, 1:, 1:, idvt:)
1796  prod => scratch(1:, 1:, 1:, iprod)
1797  sct => scratch(1:, 1:, 1:, isct)
1798  scl2 => scratch(1:, 1:, 1:, iscl2)
1799  vort => prod
1800  strain2 => prod
1801  !
1802  ! Production term.
1803  !
1804  select case (turbprod)
1805  case (strain)
1806  call prodsmag2
1807 
1808  case (vorticity)
1809  call prodwmag2
1810 
1811  case (katolaunder)
1812  call prodkatolaunder
1813 
1814  end select
1815  !
1816  ! Compute the length and time scale for all internal cells.
1817  !
1818  if (rvfb) then
1819 
1820  do k = 2, kl
1821  do j = 2, jl
1822  do i = 2, il
1823 
1824  ! Compute the time and length scale with upper bound
1825 
1826  rstrain = sqrt(strain2(i, j, k))
1827  rnu = rlv(i, j, k) / w(i, j, k, irho)
1828  tkea = abs(w(i, j, k, itu1))
1829  tepa = abs(w(i, j, k, itu2))
1830  tv2a = abs(w(i, j, k, itu3))
1831  supi = tepa * tkea / max(sqrt3 * tv2a * rvfcmu * rstrain, eps)
1832  rn2 = rvfcn**2 * (rnu * tepa)**1.5_realtype
1833 
1834  rsct = max(tkea, six * sqrt(rnu * tepa))
1835  sct(i, j, k) = min(rsct, 0.6_realtype * supi)
1836  rscl2 = tkea * min(tkea**2, supi**2)
1837  scl2(i, j, k) = rvfcl**2 * max(rscl2, rn2)
1838 
1839  end do
1840  end do
1841  end do
1842 
1843  else
1844 
1845  do k = 2, kl
1846  do j = 2, jl
1847  do i = 2, il
1848 
1849  ! Compute the time and length scale without upper bound
1850 
1851  rnu = rlv(i, j, k) / w(i, j, k, irho)
1852  tkea = abs(w(i, j, k, itu1))
1853  tepa = abs(w(i, j, k, itu2))
1854  rn2 = rvfcn**2 * (rnu * tepa)**1.5_realtype
1855 
1856  rsct = max(tkea, six * sqrt(rnu * tepa))
1857  sct(i, j, k) = rsct
1858  rscl2 = tkea**3
1859  scl2(i, j, k) = rvfcl**2 * max(rscl2, rn2)
1860  end do
1861  end do
1862  end do
1863  end if
1864 
1865  end subroutine vfscale
1866 
1867  subroutine kteddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
1868  !
1869  ! ktEddyViscosity computes the eddy viscosity according to the
1870  ! k-tau turbulence model for the block given in blockPointers.
1871  !
1872  use blockpointers
1873  use constants
1874  implicit none
1875  ! Input variables
1876  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
1877  !
1878  ! Local variables.
1879  !
1880  integer(kind=intType) :: i, j, k, ii, iSize, jSize, kSize
1881 
1882  ! Loop over the cells of this block and compute the eddy viscosity.
1883  ! Do not include halo's.
1884 #ifdef TAPENADE_REVERSE
1885  isize = (iend - ibeg) + 1
1886  jsize = (jend - jbeg) + 1
1887  ksize = (kend - kbeg) + 1
1888 
1889  !$AD II-LOOP
1890  do ii = 0, isize * jsize * ksize - 1
1891  i = mod(ii, isize) + ibeg
1892  j = mod(ii / isize, jsize) + jbeg
1893  k = ii / ((isize * jsize)) + kbeg
1894 #else
1895  do k = kbeg, kend
1896  do j = jbeg, jend
1897  do i = ibeg, iend
1898 #endif
1899  rev(i, j, k) = abs(w(i, j, k, irho) * w(i, j, k, itu1) * w(i, j, k, itu2))
1900 #ifdef TAPENADE_REVERSE
1901  end do
1902 #else
1903  end do
1904  end do
1905  end do
1906 #endif
1907  end subroutine kteddyviscosity
1908 
1909 #ifndef USE_TAPENADE
1910  subroutine unsteadyturbspectral(ntu1, ntu2)
1911  use constants
1912  use blockpointers
1913  use inputphysics
1914  use inputtimespectral
1915  use iteration
1916  use utils, only: setpointers
1917  implicit none
1918  !
1919  ! Subroutine arguments.
1920  !
1921  integer(kind=intType), intent(in) :: ntu1, ntu2
1922  !
1923  ! Local variables.
1924  !
1925  integer(kind=intType) :: nn, sps
1926 
1927  ! Return immediately if not the time spectral equations are to
1928  ! be solved.
1929 
1930  if (equationmode /= timespectral) return
1931 
1932  ! Loop over the number of spectral modes and local blocks.
1933 
1934  spectralloop: do sps = 1, ntimeintervalsspectral
1935  domains: do nn = 1, ndom
1936 
1937  ! Set the pointers for this block.
1938 
1939  call setpointers(nn, currentlevel, sps)
1940 
1941  call unsteadyturbspectral_block(ntu1, ntu2, nn, sps)
1942  end do domains
1943  end do spectralloop
1944  end subroutine unsteadyturbspectral
1945 #endif
1946 
1947  subroutine unsteadyturbspectral_block(ntu1, ntu2, nn, sps)
1948  !
1949  ! unsteadyTurbSpectral determines the spectral time derivative
1950  ! for all owned cells. This routine is called before the actual
1951  ! solve routines, such that the treatment is identical for all
1952  ! spectral solutions. The results is stored in the corresponding
1953  ! entry in dw.
1954  !
1955  use constants
1956  use blockpointers
1957  use inputphysics
1958  use inputtimespectral
1959  use iteration
1960  implicit none
1961  !
1962  ! Subroutine arguments.
1963  !
1964  integer(kind=intType), intent(in) :: ntu1, ntu2, nn, sps
1965  !
1966  ! Local variables.
1967  !
1968  integer(kind=intType) :: ii, mm, i, j, k
1969 
1970  ! Return immediately if not the time spectral equations are to
1971  ! be solved.
1972 
1973  if (equationmode /= timespectral) return
1974 
1975  ! Loop over the number of turbulent transport equations.
1976 
1977  nadvloop: do ii = ntu1, ntu2
1978 
1979  ! Initialize the time derivative to zero for the owned
1980  ! cell centers.
1981 
1982  do k = 2, kl
1983  do j = 2, jl
1984  do i = 2, il
1985  dw(i, j, k, ii) = zero
1986  end do
1987  end do
1988  end do
1989 
1990  ! Loop over the number of terms which contribute to the
1991  ! time derivative.
1992 
1993  do mm = 1, ntimeintervalsspectral
1994 
1995  ! Add the contribution to the time derivative for
1996  ! all owned cells.
1997 
1998  do k = 2, kl
1999  do j = 2, jl
2000  do i = 2, il
2001  dw(i, j, k, ii) = dw(i, j, k, ii) &
2002  + dscalar(sectionid, sps, mm) &
2003  * flowdoms(nn, currentlevel, mm)%w(i, j, k, ii)
2004  end do
2005  end do
2006  end do
2007 
2008  end do
2009 
2010  end do nadvloop
2011  end subroutine unsteadyturbspectral_block
2012 
2013  subroutine initkomega(pOffset)
2014  !
2015  ! initKOmega initializes the values of k and omega a bit more
2016  ! intelligently than just free-stream values.
2017  ! It is assumed that the pointers in blockPointers already
2018  ! point to the correct block on the correct level.
2019  ! The argument pOffset is present such that in case of restart
2020  ! a possible pointer offset is taken into account. For more
2021  ! details see the corresponding routines in initFlow.
2022  !
2023  use blockpointers
2024  use constants
2025  use flowvarrefstate
2026  use paramturb
2027  implicit none
2028  !
2029  ! Subroutine arguments.
2030  !
2031  integer(kind=intType), intent(in) :: pOffset
2032  !
2033  ! Local variables.
2034  !
2035  integer(kind=intType) :: i, j, k, ip, jp, kp
2036  real(kind=realtype) :: rhoi, tmp
2037 
2038  ! Loop over the owned cells of the block.
2039 
2040  do k = 2, kl
2041  kp = k + poffset
2042  do j = 2, jl
2043  jp = j + poffset
2044  do i = 2, il
2045  ip = i + poffset
2046 
2047  ! Compute a value of omega based on the wall distance.
2048 
2049  rhoi = one / w(ip, jp, kp, irho)
2050  tmp = six * rhoi * rlv(i, j, k) / (rkwbeta1 * (d2wall(i, j, k)**2))
2051 
2052  ! Initialize omega using the value just computed; make sure
2053  ! that the free stream value is the lowest possible value.
2054  ! After that initialize k using the current value of omega
2055  ! and the eddy viscosity. Again the free-stream value is
2056  ! the lower bound.
2057 
2058  w(ip, jp, kp, itu2) = max(tmp, winf(itu2))
2059  tmp = rhoi * rev(i, j, k) * w(ip, jp, kp, itu2)
2060  w(ip, jp, kp, itu1) = max(tmp, winf(itu1))
2061 
2062  end do
2063  end do
2064  end do
2065 
2066  end subroutine initkomega
2067 
2068  subroutine kwcdterm
2069  !
2070  ! kwCDterm computes the cross-diffusion term in the omega-eqn
2071  ! for the SST version as well as the modified k-omega turbulence
2072  ! model. It is assumed that the pointers in blockPointers and
2073  ! turbMod are already set.
2074  !
2075  use constants
2076  use blockpointers
2077  use turbmod
2078  implicit none
2079  !
2080  ! Local variables.
2081  !
2082  integer(kind=intType) :: i, j, k
2083  real(kind=realtype) :: kx, ky, kz, wwx, wwy, wwz
2084  real(kind=realtype) :: lnwip1, lnwim1, lnwjp1, lnwjm1
2085  real(kind=realtype) :: lnwkp1, lnwkm1
2086 
2087  ! Loop over the cell centers of the given block. It may be more
2088  ! efficient to loop over the faces and to scatter the gradient,
2089  ! but in that case the gradients for k and omega must be stored.
2090  ! In the current approach no extra memory is needed.
2091 
2092  do k = 2, kl
2093  do j = 2, jl
2094  do i = 2, il
2095 
2096  ! Compute the gradient of k in the cell center. Use is made
2097  ! of the fact that the surrounding normals sum up to zero,
2098  ! such that the cell i,j,k does not give a contribution.
2099  ! The gradient is scaled by a factor 1/2vol.
2100 
2101  kx = w(i + 1, j, k, itu1) * si(i, j, k, 1) - w(i - 1, j, k, itu1) * si(i - 1, j, k, 1) &
2102  + w(i, j + 1, k, itu1) * sj(i, j, k, 1) - w(i, j - 1, k, itu1) * sj(i, j - 1, k, 1) &
2103  + w(i, j, k + 1, itu1) * sk(i, j, k, 1) - w(i, j, k - 1, itu1) * sk(i, j, k - 1, 1)
2104  ky = w(i + 1, j, k, itu1) * si(i, j, k, 2) - w(i - 1, j, k, itu1) * si(i - 1, j, k, 2) &
2105  + w(i, j + 1, k, itu1) * sj(i, j, k, 2) - w(i, j - 1, k, itu1) * sj(i, j - 1, k, 2) &
2106  + w(i, j, k + 1, itu1) * sk(i, j, k, 2) - w(i, j, k - 1, itu1) * sk(i, j, k - 1, 2)
2107  kz = w(i + 1, j, k, itu1) * si(i, j, k, 3) - w(i - 1, j, k, itu1) * si(i - 1, j, k, 3) &
2108  + w(i, j + 1, k, itu1) * sj(i, j, k, 3) - w(i, j - 1, k, itu1) * sj(i, j - 1, k, 3) &
2109  + w(i, j, k + 1, itu1) * sk(i, j, k, 3) - w(i, j, k - 1, itu1) * sk(i, j, k - 1, 3)
2110 
2111  ! Compute the logarithm of omega in the points that
2112  ! contribute to the gradient in this cell.
2113 
2114  lnwip1 = log(abs(w(i + 1, j, k, itu2)))
2115  lnwim1 = log(abs(w(i - 1, j, k, itu2)))
2116  lnwjp1 = log(abs(w(i, j + 1, k, itu2)))
2117  lnwjm1 = log(abs(w(i, j - 1, k, itu2)))
2118  lnwkp1 = log(abs(w(i, j, k + 1, itu2)))
2119  lnwkm1 = log(abs(w(i, j, k - 1, itu2)))
2120 
2121  ! Compute the scaled gradient of ln omega.
2122 
2123  wwx = lnwip1 * si(i, j, k, 1) - lnwim1 * si(i - 1, j, k, 1) &
2124  + lnwjp1 * sj(i, j, k, 1) - lnwjm1 * sj(i, j - 1, k, 1) &
2125  + lnwkp1 * sk(i, j, k, 1) - lnwkm1 * sk(i, j, k - 1, 1)
2126  wwy = lnwip1 * si(i, j, k, 2) - lnwim1 * si(i - 1, j, k, 2) &
2127  + lnwjp1 * sj(i, j, k, 2) - lnwjm1 * sj(i, j - 1, k, 2) &
2128  + lnwkp1 * sk(i, j, k, 2) - lnwkm1 * sk(i, j, k - 1, 2)
2129  wwz = lnwip1 * si(i, j, k, 3) - lnwim1 * si(i - 1, j, k, 3) &
2130  + lnwjp1 * sj(i, j, k, 3) - lnwjm1 * sj(i, j - 1, k, 3) &
2131  + lnwkp1 * sk(i, j, k, 3) - lnwkm1 * sk(i, j, k - 1, 3)
2132 
2133  ! Compute the dot product grad k grad ln omega.
2134  ! Multiply it by the correct scaling factor and store it.
2135 
2136  kwcd(i, j, k) = fourth * (kx * wwx + ky * wwy + kz * wwz) / (vol(i, j, k)**2)
2137 
2138  end do
2139  end do
2140  end do
2141 
2142  end subroutine kwcdterm
2143 #endif
2144 end module turbutils
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::nx
integer(kind=inttype) nx
Definition: blockPointers.F90:38
blockpointers::bmti2
real(kind=realtype), dimension(:, :, :, :), pointer bmti2
Definition: blockPointers.F90:158
blockpointers::ny
integer(kind=inttype) ny
Definition: blockPointers.F90:38
blockpointers::ie
integer(kind=inttype) ie
Definition: blockPointers.F90:39
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::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::sj
real(kind=realtype), dimension(:, :, :, :), pointer sj
Definition: blockPointers.F90:96
blockpointers::sectionid
integer(kind=inttype) sectionid
Definition: blockPointers.F90:30
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::nz
integer(kind=inttype) nz
Definition: blockPointers.F90:38
blockpointers::je
integer(kind=inttype) je
Definition: blockPointers.F90:39
blockpointers::ke
integer(kind=inttype) ke
Definition: blockPointers.F90:39
blockpointers::bmtk1
real(kind=realtype), dimension(:, :, :, :), pointer bmtk1
Definition: blockPointers.F90:161
blockpointers::wold
real(kind=realtype), dimension(:, :, :, :, :), pointer wold
Definition: blockPointers.F90:121
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::spalartallmarasedwards
integer(kind=inttype), parameter spalartallmarasedwards
Definition: constants.F90:128
constants::spalartallmaras
integer(kind=inttype), parameter spalartallmaras
Definition: constants.F90:128
constants::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::three
real(kind=realtype), parameter three
Definition: constants.F90:74
constants::imax
integer(kind=inttype), parameter imax
Definition: constants.F90:293
constants::itu3
integer, parameter itu3
Definition: constants.F90:43
constants::idvt
integer, parameter idvt
Definition: constants.F90:51
constants::ivort
integer, parameter ivort
Definition: constants.F90:53
constants::kmin
integer(kind=inttype), parameter kmin
Definition: constants.F90:296
constants::four
real(kind=realtype), parameter four
Definition: constants.F90:75
constants::third
real(kind=realtype), parameter third
Definition: constants.F90:81
constants::jmax
integer(kind=inttype), parameter jmax
Definition: constants.F90:295
constants::itu2
integer, parameter itu2
Definition: constants.F90:42
constants::six
real(kind=realtype), parameter six
Definition: constants.F90:77
constants::isct
integer, parameter isct
Definition: constants.F90:58
constants::eps
real(kind=realtype), parameter eps
Definition: constants.F90:23
constants::thresholdreal
real(kind=realtype), parameter thresholdreal
Definition: constants.F90:101
constants::ktau
integer(kind=inttype), parameter ktau
Definition: constants.F90:128
constants::itu1
integer, parameter itu1
Definition: constants.F90:40
constants::itu5
integer, parameter itu5
Definition: constants.F90:45
constants::irho
integer, parameter irho
Definition: constants.F90:34
constants::vorticity
integer(kind=inttype), parameter vorticity
Definition: constants.F90:137
constants::komegawilcox
integer(kind=inttype), parameter komegawilcox
Definition: constants.F90:128
constants::timespectral
integer(kind=inttype), parameter timespectral
Definition: constants.F90:115
constants::komegamodified
integer(kind=inttype), parameter komegamodified
Definition: constants.F90:128
constants::ivx
integer, parameter ivx
Definition: constants.F90:35
constants::iprod
integer, parameter iprod
Definition: constants.F90:55
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
constants::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::mentersst
integer(kind=inttype), parameter mentersst
Definition: constants.F90:128
constants::secondorder
integer(kind=inttype), parameter secondorder
Definition: constants.F90:142
constants::kmax
integer(kind=inttype), parameter kmax
Definition: constants.F90:297
constants::ivy
integer, parameter ivy
Definition: constants.F90:36
constants::jmin
integer(kind=inttype), parameter jmin
Definition: constants.F90:294
constants::iscl2
integer, parameter iscl2
Definition: constants.F90:59
constants::v2f
integer(kind=inttype), parameter v2f
Definition: constants.F90:128
flowvarrefstate::winf
real(kind=realtype), dimension(:), allocatable winf
Definition: flowVarRefState.F90:74
flowvarrefstate::timeref
real(kind=realtype) timeref
Definition: flowVarRefState.F90:40
inputdiscretization::orderturb
integer(kind=inttype) orderturb
Definition: inputParam.F90:73
inputphysics::equationmode
integer(kind=inttype) equationmode
Definition: inputParam.F90:583
inputphysics::rvfb
logical rvfb
Definition: inputParam.F90:586
inputphysics::turbprod
integer(kind=inttype) turbprod
Definition: inputParam.F90:584
inputphysics::turbmodel
integer(kind=inttype) turbmodel
Definition: inputParam.F90:584
inputphysics::wallfunctions
logical wallfunctions
Definition: inputParam.F90:589
inputtimespectral::ntimeintervalsspectral
integer(kind=inttype) ntimeintervalsspectral
Definition: inputParam.F90:645
inputtimespectral::dscalar
real(kind=realtype), dimension(:, :, :), allocatable dscalar
Definition: inputParam.F90:659
inputunsteady::deltat
real(kind=realtype) deltat
Definition: inputParam.F90:733
inputunsteady::timeintegrationscheme
integer(kind=inttype) timeintegrationscheme
Definition: inputParam.F90:719
iteration::noldlevels
integer(kind=inttype) noldlevels
Definition: iteration.f90:79
iteration::currentlevel
integer(kind=inttype) currentlevel
Definition: iteration.f90:18
iteration::groundlevel
integer(kind=inttype) groundlevel
Definition: iteration.f90:18
iteration::coeftime
real(kind=realtype), dimension(:), allocatable coeftime
Definition: iteration.f90:80
paramturb::rssta1
real(kind=realtype), parameter rssta1
Definition: paramTurb.F90:41
paramturb::rvflimite
real(kind=realtype) rvflimite
Definition: paramTurb.F90:80
paramturb::rvfcmu
real(kind=realtype) rvfcmu
Definition: paramTurb.F90:81
paramturb::rvfcn
real(kind=realtype), parameter rvfcn
Definition: paramTurb.F90:69
paramturb::rvfcl
real(kind=realtype) rvfcl
Definition: paramTurb.F90:80
paramturb::rsacv1
real(kind=realtype), parameter rsacv1
Definition: paramTurb.F90:17
paramturb::rkwbeta1
real(kind=realtype), parameter rkwbeta1
Definition: paramTurb.F90:35
section::sections
type(sectiontype), dimension(:), allocatable sections
Definition: section.f90:46
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::scl2
real(kind=realtype), dimension(:, :, :), pointer scl2
Definition: turbMod.F90:40
turbmod::kwcd
real(kind=realtype), dimension(:, :, :), pointer kwcd
Definition: turbMod.F90:37
turbmod::prod
real(kind=realtype), dimension(:, :, :), pointer prod
Definition: turbMod.F90:35
turbmod::secondord
logical secondord
Definition: turbMod.F90:15
turbmod::strain2
real(kind=realtype), dimension(:, :, :), pointer strain2
Definition: turbMod.F90:41
turbmod::dvt
real(kind=realtype), dimension(:, :, :, :), pointer dvt
Definition: turbMod.F90:33
turbmod::sct
real(kind=realtype), dimension(:, :, :), pointer sct
Definition: turbMod.F90:39
turbutils::vfeddyviscosity
subroutine vfeddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:1635
turbutils::vfscale
subroutine vfscale
Definition: turbUtils.F90:1767
turbutils::initkomega
subroutine initkomega(pOffset)
Definition: turbUtils.F90:2014
turbutils::unsteadyturbterm
subroutine unsteadyturbterm(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:411
turbutils::sanuknowneddyratio
real(kind=realtype) function sanuknowneddyratio(eddyRatio, nuLam)
Definition: turbUtils.F90:334
turbutils::saeddyviscosity
subroutine saeddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:657
turbutils::computeeddyviscosity
subroutine computeeddyviscosity(includeHalos)
Definition: turbUtils.F90:581
turbutils::ssteddyviscosity
subroutine ssteddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:757
turbutils::kwcdterm
subroutine kwcdterm
Definition: turbUtils.F90:2069
turbutils::prodwmag2
subroutine prodwmag2
Definition: turbUtils.F90:240
turbutils::tdia3
subroutine tdia3(nb, ne, l, c, u, r)
Definition: turbUtils.F90:1561
turbutils::turbadvection
subroutine turbadvection(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:826
turbutils::prodsmag2
subroutine prodsmag2
Definition: turbUtils.F90:130
turbutils::prodkatolaunder
subroutine prodkatolaunder
Definition: turbUtils.F90:6
turbutils::unsteadyturbspectral
subroutine unsteadyturbspectral(ntu1, ntu2)
Definition: turbUtils.F90:1911
turbutils::kweddyviscosity
subroutine kweddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:713
turbutils::unsteadyturbspectral_block
subroutine unsteadyturbspectral_block(ntu1, ntu2, nn, sps)
Definition: turbUtils.F90:1948
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