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