Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
sx_lambda2.f90
Go to the documentation of this file.
1! Copyright (c) 2024, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
34submodule(opr_sx) sx_lambda2
35 use math, only : pi
36 implicit none
37
38contains
39
40 module subroutine opr_sx_lambda2(lambda2, u, v, w, coef)
41 type(coef_t), intent(in) :: coef
42 type(field_t), intent(inout) :: lambda2
43 type(field_t), intent(in) :: u, v, w
44
45 associate(xh => coef%Xh, msh => coef%msh)
46 select case(xh%lx)
47 case (18)
48 call sx_lambda2_lx18(lambda2%x, u%x, v%x, w%x, &
49 xh%dx, xh%dy, xh%dz, &
50 coef%drdx, coef%dsdx, coef%dtdx, &
51 coef%drdy, coef%dsdy, coef%dtdy, &
52 coef%drdz, coef%dsdz, coef%dtdz, &
53 xh%w3, coef%B, msh%nelv)
54 case (17)
55 call sx_lambda2_lx17(lambda2%x, u%x, v%x, w%x, &
56 xh%dx, xh%dy, xh%dz, &
57 coef%drdx, coef%dsdx, coef%dtdx, &
58 coef%drdy, coef%dsdy, coef%dtdy, &
59 coef%drdz, coef%dsdz, coef%dtdz, &
60 xh%w3, coef%B, msh%nelv)
61 case (16)
62 call sx_lambda2_lx16(lambda2%x, u%x, v%x, w%x, &
63 xh%dx, xh%dy, xh%dz, &
64 coef%drdx, coef%dsdx, coef%dtdx, &
65 coef%drdy, coef%dsdy, coef%dtdy, &
66 coef%drdz, coef%dsdz, coef%dtdz, &
67 xh%w3, coef%B, msh%nelv)
68 case (15)
69 call sx_lambda2_lx15(lambda2%x, u%x, v%x, w%x, &
70 xh%dx, xh%dy, xh%dz, &
71 coef%drdx, coef%dsdx, coef%dtdx, &
72 coef%drdy, coef%dsdy, coef%dtdy, &
73 coef%drdz, coef%dsdz, coef%dtdz, &
74 xh%w3, coef%B, msh%nelv)
75 case (14)
76 call sx_lambda2_lx14(lambda2%x, u%x, v%x, w%x, &
77 xh%dx, xh%dy, xh%dz, &
78 coef%drdx, coef%dsdx, coef%dtdx, &
79 coef%drdy, coef%dsdy, coef%dtdy, &
80 coef%drdz, coef%dsdz, coef%dtdz, &
81 xh%w3, coef%B, msh%nelv)
82 case (13)
83 call sx_lambda2_lx13(lambda2%x, u%x, v%x, w%x, &
84 xh%dx, xh%dy, xh%dz, &
85 coef%drdx, coef%dsdx, coef%dtdx, &
86 coef%drdy, coef%dsdy, coef%dtdy, &
87 coef%drdz, coef%dsdz, coef%dtdz, &
88 xh%w3, coef%B, msh%nelv)
89 case (12)
90 call sx_lambda2_lx12(lambda2%x, u%x, v%x, w%x, &
91 xh%dx, xh%dy, xh%dz, &
92 coef%drdx, coef%dsdx, coef%dtdx, &
93 coef%drdy, coef%dsdy, coef%dtdy, &
94 coef%drdz, coef%dsdz, coef%dtdz, &
95 xh%w3, coef%B, msh%nelv)
96 case (11)
97 call sx_lambda2_lx11(lambda2%x, u%x, v%x, w%x, &
98 xh%dx, xh%dy, xh%dz, &
99 coef%drdx, coef%dsdx, coef%dtdx, &
100 coef%drdy, coef%dsdy, coef%dtdy, &
101 coef%drdz, coef%dsdz, coef%dtdz, &
102 xh%w3, coef%B, msh%nelv)
103 case (10)
104 call sx_lambda2_lx10(lambda2%x, u%x, v%x, w%x, &
105 xh%dx, xh%dy, xh%dz, &
106 coef%drdx, coef%dsdx, coef%dtdx, &
107 coef%drdy, coef%dsdy, coef%dtdy, &
108 coef%drdz, coef%dsdz, coef%dtdz, &
109 xh%w3, coef%B, msh%nelv)
110 case (9)
111 call sx_lambda2_lx9(lambda2%x, u%x, v%x, w%x, &
112 xh%dx, xh%dy, xh%dz, &
113 coef%drdx, coef%dsdx, coef%dtdx, &
114 coef%drdy, coef%dsdy, coef%dtdy, &
115 coef%drdz, coef%dsdz, coef%dtdz, &
116 xh%w3, coef%B, msh%nelv)
117 case (8)
118 call sx_lambda2_lx8(lambda2%x, u%x, v%x, w%x, &
119 xh%dx, xh%dy, xh%dz, &
120 coef%drdx, coef%dsdx, coef%dtdx, &
121 coef%drdy, coef%dsdy, coef%dtdy, &
122 coef%drdz, coef%dsdz, coef%dtdz, &
123 xh%w3, coef%B, msh%nelv)
124 case (7)
125 call sx_lambda2_lx7(lambda2%x, u%x, v%x, w%x, &
126 xh%dx, xh%dy, xh%dz, &
127 coef%drdx, coef%dsdx, coef%dtdx, &
128 coef%drdy, coef%dsdy, coef%dtdy, &
129 coef%drdz, coef%dsdz, coef%dtdz, &
130 xh%w3, coef%B, msh%nelv)
131 case (6)
132 call sx_lambda2_lx6(lambda2%x, u%x, v%x, w%x, &
133 xh%dx, xh%dy, xh%dz, &
134 coef%drdx, coef%dsdx, coef%dtdx, &
135 coef%drdy, coef%dsdy, coef%dtdy, &
136 coef%drdz, coef%dsdz, coef%dtdz, &
137 xh%w3, coef%B, msh%nelv)
138 case (5)
139 call sx_lambda2_lx5(lambda2%x, u%x, v%x, w%x, &
140 xh%dx, xh%dy, xh%dz, &
141 coef%drdx, coef%dsdx, coef%dtdx, &
142 coef%drdy, coef%dsdy, coef%dtdy, &
143 coef%drdz, coef%dsdz, coef%dtdz, &
144 xh%w3, coef%B, msh%nelv)
145 case (4)
146 call sx_lambda2_lx4(lambda2%x, u%x, v%x, w%x, &
147 xh%dx, xh%dy, xh%dz, &
148 coef%drdx, coef%dsdx, coef%dtdx, &
149 coef%drdy, coef%dsdy, coef%dtdy, &
150 coef%drdz, coef%dsdz, coef%dtdz, &
151 xh%w3, coef%B, msh%nelv)
152 case (3)
153 call sx_lambda2_lx3(lambda2%x, u%x, v%x, w%x, &
154 xh%dx, xh%dy, xh%dz, &
155 coef%drdx, coef%dsdx, coef%dtdx, &
156 coef%drdy, coef%dsdy, coef%dtdy, &
157 coef%drdz, coef%dsdz, coef%dtdz, &
158 xh%w3, coef%B, msh%nelv)
159 case (2)
160 call sx_lambda2_lx2(lambda2%x, u%x, v%x, w%x, &
161 xh%dx, xh%dy, xh%dz, &
162 coef%drdx, coef%dsdx, coef%dtdx, &
163 coef%drdy, coef%dsdy, coef%dtdy, &
164 coef%drdz, coef%dsdz, coef%dtdz, &
165 xh%w3, coef%B, msh%nelv)
166 case default
167 call sx_lambda2_lx(lambda2%x, u%x, v%x, w%x, &
168 xh%dx, xh%dy, xh%dz, &
169 coef%drdx, coef%dsdx, coef%dtdx, &
170 coef%drdy, coef%dsdy, coef%dtdy, &
171 coef%drdz, coef%dsdz, coef%dtdz, &
172 xh%w3, coef%B, msh%nelv, xh%lx)
173 end select
174 end associate
175
176 end subroutine opr_sx_lambda2
177
178 subroutine sx_lambda2_lx(lambda2, u, v, w, dx, dy, dz, &
179 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n, lx)
180 integer, intent(in) :: n, lx
181 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
182 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
183 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
184 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
185 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
186 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
187 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
188 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
189 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
190 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
191 real(kind=rp) :: grad(lx*lx*lx,3,3)
192 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
193 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
194 real(kind=rp) :: a11, a22, a33, a12, a13, a23
195 real(kind=rp) :: msk1, msk2, msk3
196 real(kind=rp) :: ur(lx, lx, lx)
197 real(kind=rp) :: vr(lx, lx, lx)
198 real(kind=rp) :: wr(lx, lx, lx)
199 real(kind=rp) :: us(lx, lx, lx)
200 real(kind=rp) :: vs(lx, lx, lx)
201 real(kind=rp) :: ws(lx, lx, lx)
202 real(kind=rp) :: ut(lx, lx, lx)
203 real(kind=rp) :: vt(lx, lx, lx)
204 real(kind=rp) :: wt(lx, lx, lx)
205 real(kind=rp) :: tmp1, tmp2, tmp3
206 integer :: e, i, j, k, l
207
208 do e = 1, n
209 do j = 1, lx * lx
210 do i = 1, lx
211 tmp1 = 0.0_rp
212 tmp2 = 0.0_rp
213 tmp3 = 0.0_rp
214 do k = 1, lx
215 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
216 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
217 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
218 end do
219 ur(i,j,1) = tmp1
220 vr(i,j,1) = tmp2
221 wr(i,j,1) = tmp3
222 end do
223 end do
224
225 do k = 1, lx
226 do j = 1, lx
227 do i = 1, lx
228 tmp1 = 0.0_rp
229 tmp2 = 0.0_rp
230 tmp3 = 0.0_rp
231 do l = 1, lx
232 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
233 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
234 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
235 end do
236 us(i,j,k) = tmp1
237 vs(i,j,k) = tmp2
238 ws(i,j,k) = tmp3
239 end do
240 end do
241 end do
242
243 do k = 1, lx
244 do i = 1, lx*lx
245 tmp1 = 0.0_rp
246 tmp2 = 0.0_rp
247 tmp3 = 0.0_rp
248 do l = 1, lx
249 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
250 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
251 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
252 end do
253 ut(i,1,k) = tmp1
254 vt(i,1,k) = tmp2
255 wt(i,1,k) = tmp3
256 end do
257 end do
258
259 do i = 1, lx * lx * lx
260 grad(1,1,1) = w3(i,1,1) &
261 * ( drdx(i,1,1,e) * ur(i,1,1) &
262 + dsdx(i,1,1,e) * us(i,1,1) &
263 + dtdx(i,1,1,e) * ut(i,1,1) )
264 grad(1,1,2) = w3(i,1,1) &
265 * ( dsdy(i,1,1,e) * us(i,1,1) &
266 + drdy(i,1,1,e) * ur(i,1,1) &
267 + dtdy(i,1,1,e) * ut(i,1,1) )
268 grad(1,1,3) = w3(i,1,1) &
269 * ( dtdz(i,1,1,e) * ut(i,1,1) &
270 + drdz(i,1,1,e) * ur(i,1,1) &
271 + dsdz(i,1,1,e) * us(i,1,1) )
272
273 grad(1,2,1) = w3(i,1,1) &
274 * ( drdx(i,1,1,e) * vr(i,1,1) &
275 + dsdx(i,1,1,e) * vs(i,1,1) &
276 + dtdx(i,1,1,e) * vt(i,1,1) )
277 grad(1,2,2) = w3(i,1,1) &
278 * ( dsdy(i,1,1,e) * vs(i,1,1) &
279 + drdy(i,1,1,e) * vr(i,1,1) &
280 + dtdy(i,1,1,e) * vt(i,1,1) )
281 grad(1,2,3) = w3(i,1,1) &
282 * ( dtdz(i,1,1,e) * vt(i,1,1) &
283 + drdz(i,1,1,e) * vr(i,1,1) &
284 + dsdz(i,1,1,e) * vs(i,1,1) )
285
286 grad(1,3,1) = w3(i,1,1) &
287 * ( drdx(i,1,1,e) * wr(i,1,1) &
288 + dsdx(i,1,1,e) * ws(i,1,1) &
289 + dtdx(i,1,1,e) * wt(i,1,1) )
290 grad(1,3,2) = w3(i,1,1) &
291 * ( dsdy(i,1,1,e) * ws(i,1,1) &
292 + drdy(i,1,1,e) * wr(i,1,1) &
293 + dtdy(i,1,1,e) * wt(i,1,1) )
294 grad(1,3,3) = w3(i,1,1) &
295 * ( dtdz(i,1,1,e) * wt(i,1,1) &
296 + drdz(i,1,1,e) * wr(i,1,1) &
297 + dsdz(i,1,1,e) * ws(i,1,1) )
298 end do
299
300
301 do i = 1, lx * lx * lx
302 s11 = grad(i,1,1)
303 s22 = grad(i,2,2)
304 s33 = grad(i,3,3)
305
306
307 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
308 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
309 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
310
311 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
312 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
313 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
314
315 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
316 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
317 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
318
319 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
320 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
321 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
322
323
324 b = -(a11 + a22 + a33)
325 c = -(a12*a12 + a13*a13 + a23*a23 &
326 - a11 * a22 - a11 * a33 - a22 * a33)
327 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
328 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
329
330
331 q = (3.0 * c - b*b) / 9.0
332 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
333 theta = acos( r / sqrt(-q*q*q) )
334
335 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
336 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
337 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
338
339 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
340 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
341 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
342 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
343 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
344 .le. eigen(2) .and. eigen(2) .le. eigen(1))
345 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
346 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
347 .le. eigen(3) .and. eigen(3) .le. eigen(1))
348
349 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
350
351 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
352 end do
353 end do
354 end subroutine sx_lambda2_lx
355
356 subroutine sx_lambda2_lx18(lambda2, u, v, w, dx, dy, dz, &
357 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
358 integer, parameter :: lx = 18
359 integer, intent(in) :: n
360 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
361 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
362 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
363 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
364 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
365 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
366 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
367 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
368 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
369 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
370 real(kind=rp) :: grad(lx*lx*lx,3,3)
371 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
372 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
373 real(kind=rp) :: a11, a22, a33, a12, a13, a23
374 real(kind=rp) :: msk1, msk2, msk3
375 real(kind=rp) :: ur(lx, lx, lx)
376 real(kind=rp) :: vr(lx, lx, lx)
377 real(kind=rp) :: wr(lx, lx, lx)
378 real(kind=rp) :: us(lx, lx, lx)
379 real(kind=rp) :: vs(lx, lx, lx)
380 real(kind=rp) :: ws(lx, lx, lx)
381 real(kind=rp) :: ut(lx, lx, lx)
382 real(kind=rp) :: vt(lx, lx, lx)
383 real(kind=rp) :: wt(lx, lx, lx)
384 real(kind=rp) :: tmp1, tmp2, tmp3
385 integer :: e, i, j, k, l
386
387 do e = 1, n
388 do j = 1, lx * lx
389 do i = 1, lx
390 tmp1 = 0.0_rp
391 tmp2 = 0.0_rp
392 tmp3 = 0.0_rp
393 do k = 1, lx
394 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
395 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
396 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
397 end do
398 ur(i,j,1) = tmp1
399 vr(i,j,1) = tmp2
400 wr(i,j,1) = tmp3
401 end do
402 end do
403
404 do k = 1, lx
405 do j = 1, lx
406 do i = 1, lx
407 tmp1 = 0.0_rp
408 tmp2 = 0.0_rp
409 tmp3 = 0.0_rp
410 do l = 1, lx
411 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
412 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
413 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
414 end do
415 us(i,j,k) = tmp1
416 vs(i,j,k) = tmp2
417 ws(i,j,k) = tmp3
418 end do
419 end do
420 end do
421
422 do k = 1, lx
423 do i = 1, lx*lx
424 tmp1 = 0.0_rp
425 tmp2 = 0.0_rp
426 tmp3 = 0.0_rp
427 do l = 1, lx
428 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
429 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
430 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
431 end do
432 ut(i,1,k) = tmp1
433 vt(i,1,k) = tmp2
434 wt(i,1,k) = tmp3
435 end do
436 end do
437
438 do i = 1, lx * lx * lx
439 grad(1,1,1) = w3(i,1,1) &
440 * ( drdx(i,1,1,e) * ur(i,1,1) &
441 + dsdx(i,1,1,e) * us(i,1,1) &
442 + dtdx(i,1,1,e) * ut(i,1,1) )
443 grad(1,1,2) = w3(i,1,1) &
444 * ( dsdy(i,1,1,e) * us(i,1,1) &
445 + drdy(i,1,1,e) * ur(i,1,1) &
446 + dtdy(i,1,1,e) * ut(i,1,1) )
447 grad(1,1,3) = w3(i,1,1) &
448 * ( dtdz(i,1,1,e) * ut(i,1,1) &
449 + drdz(i,1,1,e) * ur(i,1,1) &
450 + dsdz(i,1,1,e) * us(i,1,1) )
451
452 grad(1,2,1) = w3(i,1,1) &
453 * ( drdx(i,1,1,e) * vr(i,1,1) &
454 + dsdx(i,1,1,e) * vs(i,1,1) &
455 + dtdx(i,1,1,e) * vt(i,1,1) )
456 grad(1,2,2) = w3(i,1,1) &
457 * ( dsdy(i,1,1,e) * vs(i,1,1) &
458 + drdy(i,1,1,e) * vr(i,1,1) &
459 + dtdy(i,1,1,e) * vt(i,1,1) )
460 grad(1,2,3) = w3(i,1,1) &
461 * ( dtdz(i,1,1,e) * vt(i,1,1) &
462 + drdz(i,1,1,e) * vr(i,1,1) &
463 + dsdz(i,1,1,e) * vs(i,1,1) )
464
465 grad(1,3,1) = w3(i,1,1) &
466 * ( drdx(i,1,1,e) * wr(i,1,1) &
467 + dsdx(i,1,1,e) * ws(i,1,1) &
468 + dtdx(i,1,1,e) * wt(i,1,1) )
469 grad(1,3,2) = w3(i,1,1) &
470 * ( dsdy(i,1,1,e) * ws(i,1,1) &
471 + drdy(i,1,1,e) * wr(i,1,1) &
472 + dtdy(i,1,1,e) * wt(i,1,1) )
473 grad(1,3,3) = w3(i,1,1) &
474 * ( dtdz(i,1,1,e) * wt(i,1,1) &
475 + drdz(i,1,1,e) * wr(i,1,1) &
476 + dsdz(i,1,1,e) * ws(i,1,1) )
477 end do
478
479
480 do i = 1, lx * lx * lx
481 s11 = grad(i,1,1)
482 s22 = grad(i,2,2)
483 s33 = grad(i,3,3)
484
485
486 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
487 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
488 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
489
490 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
491 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
492 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
493
494 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
495 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
496 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
497
498 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
499 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
500 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
501
502
503 b = -(a11 + a22 + a33)
504 c = -(a12*a12 + a13*a13 + a23*a23 &
505 - a11 * a22 - a11 * a33 - a22 * a33)
506 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
507 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
508
509
510 q = (3.0 * c - b*b) / 9.0
511 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
512 theta = acos( r / sqrt(-q*q*q) )
513
514 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
515 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
516 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
517
518 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
519 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
520 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
521 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
522 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
523 .le. eigen(2) .and. eigen(2) .le. eigen(1))
524 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
525 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
526 .le. eigen(3) .and. eigen(3) .le. eigen(1))
527
528 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
529
530 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
531 end do
532 end do
533 end subroutine sx_lambda2_lx18
534
535 subroutine sx_lambda2_lx17(lambda2, u, v, w, dx, dy, dz, &
536 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
537 integer, parameter :: lx = 17
538 integer, intent(in) :: n
539 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
540 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
541 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
542 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
543 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
544 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
545 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
546 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
547 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
548 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
549 real(kind=rp) :: grad(lx*lx*lx,3,3)
550 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
551 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
552 real(kind=rp) :: a11, a22, a33, a12, a13, a23
553 real(kind=rp) :: msk1, msk2, msk3
554 real(kind=rp) :: ur(lx, lx, lx)
555 real(kind=rp) :: vr(lx, lx, lx)
556 real(kind=rp) :: wr(lx, lx, lx)
557 real(kind=rp) :: us(lx, lx, lx)
558 real(kind=rp) :: vs(lx, lx, lx)
559 real(kind=rp) :: ws(lx, lx, lx)
560 real(kind=rp) :: ut(lx, lx, lx)
561 real(kind=rp) :: vt(lx, lx, lx)
562 real(kind=rp) :: wt(lx, lx, lx)
563 real(kind=rp) :: tmp1, tmp2, tmp3
564 integer :: e, i, j, k, l
565
566 do e = 1, n
567 do j = 1, lx * lx
568 do i = 1, lx
569 tmp1 = 0.0_rp
570 tmp2 = 0.0_rp
571 tmp3 = 0.0_rp
572 do k = 1, lx
573 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
574 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
575 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
576 end do
577 ur(i,j,1) = tmp1
578 vr(i,j,1) = tmp2
579 wr(i,j,1) = tmp3
580 end do
581 end do
582
583 do k = 1, lx
584 do j = 1, lx
585 do i = 1, lx
586 tmp1 = 0.0_rp
587 tmp2 = 0.0_rp
588 tmp3 = 0.0_rp
589 do l = 1, lx
590 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
591 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
592 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
593 end do
594 us(i,j,k) = tmp1
595 vs(i,j,k) = tmp2
596 ws(i,j,k) = tmp3
597 end do
598 end do
599 end do
600
601 do k = 1, lx
602 do i = 1, lx*lx
603 tmp1 = 0.0_rp
604 tmp2 = 0.0_rp
605 tmp3 = 0.0_rp
606 do l = 1, lx
607 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
608 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
609 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
610 end do
611 ut(i,1,k) = tmp1
612 vt(i,1,k) = tmp2
613 wt(i,1,k) = tmp3
614 end do
615 end do
616
617 do i = 1, lx * lx * lx
618 grad(1,1,1) = w3(i,1,1) &
619 * ( drdx(i,1,1,e) * ur(i,1,1) &
620 + dsdx(i,1,1,e) * us(i,1,1) &
621 + dtdx(i,1,1,e) * ut(i,1,1) )
622 grad(1,1,2) = w3(i,1,1) &
623 * ( dsdy(i,1,1,e) * us(i,1,1) &
624 + drdy(i,1,1,e) * ur(i,1,1) &
625 + dtdy(i,1,1,e) * ut(i,1,1) )
626 grad(1,1,3) = w3(i,1,1) &
627 * ( dtdz(i,1,1,e) * ut(i,1,1) &
628 + drdz(i,1,1,e) * ur(i,1,1) &
629 + dsdz(i,1,1,e) * us(i,1,1) )
630
631 grad(1,2,1) = w3(i,1,1) &
632 * ( drdx(i,1,1,e) * vr(i,1,1) &
633 + dsdx(i,1,1,e) * vs(i,1,1) &
634 + dtdx(i,1,1,e) * vt(i,1,1) )
635 grad(1,2,2) = w3(i,1,1) &
636 * ( dsdy(i,1,1,e) * vs(i,1,1) &
637 + drdy(i,1,1,e) * vr(i,1,1) &
638 + dtdy(i,1,1,e) * vt(i,1,1) )
639 grad(1,2,3) = w3(i,1,1) &
640 * ( dtdz(i,1,1,e) * vt(i,1,1) &
641 + drdz(i,1,1,e) * vr(i,1,1) &
642 + dsdz(i,1,1,e) * vs(i,1,1) )
643
644 grad(1,3,1) = w3(i,1,1) &
645 * ( drdx(i,1,1,e) * wr(i,1,1) &
646 + dsdx(i,1,1,e) * ws(i,1,1) &
647 + dtdx(i,1,1,e) * wt(i,1,1) )
648 grad(1,3,2) = w3(i,1,1) &
649 * ( dsdy(i,1,1,e) * ws(i,1,1) &
650 + drdy(i,1,1,e) * wr(i,1,1) &
651 + dtdy(i,1,1,e) * wt(i,1,1) )
652 grad(1,3,3) = w3(i,1,1) &
653 * ( dtdz(i,1,1,e) * wt(i,1,1) &
654 + drdz(i,1,1,e) * wr(i,1,1) &
655 + dsdz(i,1,1,e) * ws(i,1,1) )
656 end do
657
658
659 do i = 1, lx * lx * lx
660 s11 = grad(i,1,1)
661 s22 = grad(i,2,2)
662 s33 = grad(i,3,3)
663
664
665 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
666 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
667 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
668
669 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
670 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
671 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
672
673 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
674 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
675 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
676
677 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
678 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
679 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
680
681
682 b = -(a11 + a22 + a33)
683 c = -(a12*a12 + a13*a13 + a23*a23 &
684 - a11 * a22 - a11 * a33 - a22 * a33)
685 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
686 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
687
688
689 q = (3.0 * c - b*b) / 9.0
690 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
691 theta = acos( r / sqrt(-q*q*q) )
692
693 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
694 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
695 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
696
697 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
698 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
699 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
700 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
701 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
702 .le. eigen(2) .and. eigen(2) .le. eigen(1))
703 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
704 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
705 .le. eigen(3) .and. eigen(3) .le. eigen(1))
706
707 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
708
709 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
710 end do
711 end do
712 end subroutine sx_lambda2_lx17
713
714 subroutine sx_lambda2_lx16(lambda2, u, v, w, dx, dy, dz, &
715 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
716 integer, parameter :: lx = 16
717 integer, intent(in) :: n
718 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
719 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
720 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
721 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
722 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
723 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
724 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
725 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
726 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
727 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
728 real(kind=rp) :: grad(lx*lx*lx,3,3)
729 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
730 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
731 real(kind=rp) :: a11, a22, a33, a12, a13, a23
732 real(kind=rp) :: msk1, msk2, msk3
733 real(kind=rp) :: ur(lx, lx, lx)
734 real(kind=rp) :: vr(lx, lx, lx)
735 real(kind=rp) :: wr(lx, lx, lx)
736 real(kind=rp) :: us(lx, lx, lx)
737 real(kind=rp) :: vs(lx, lx, lx)
738 real(kind=rp) :: ws(lx, lx, lx)
739 real(kind=rp) :: ut(lx, lx, lx)
740 real(kind=rp) :: vt(lx, lx, lx)
741 real(kind=rp) :: wt(lx, lx, lx)
742 real(kind=rp) :: tmp1, tmp2, tmp3
743 integer :: e, i, j, k, l
744
745 do e = 1, n
746 do j = 1, lx * lx
747 do i = 1, lx
748 tmp1 = 0.0_rp
749 tmp2 = 0.0_rp
750 tmp3 = 0.0_rp
751 do k = 1, lx
752 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
753 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
754 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
755 end do
756 ur(i,j,1) = tmp1
757 vr(i,j,1) = tmp2
758 wr(i,j,1) = tmp3
759 end do
760 end do
761
762 do k = 1, lx
763 do j = 1, lx
764 do i = 1, lx
765 tmp1 = 0.0_rp
766 tmp2 = 0.0_rp
767 tmp3 = 0.0_rp
768 do l = 1, lx
769 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
770 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
771 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
772 end do
773 us(i,j,k) = tmp1
774 vs(i,j,k) = tmp2
775 ws(i,j,k) = tmp3
776 end do
777 end do
778 end do
779
780 do k = 1, lx
781 do i = 1, lx*lx
782 tmp1 = 0.0_rp
783 tmp2 = 0.0_rp
784 tmp3 = 0.0_rp
785 do l = 1, lx
786 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
787 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
788 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
789 end do
790 ut(i,1,k) = tmp1
791 vt(i,1,k) = tmp2
792 wt(i,1,k) = tmp3
793 end do
794 end do
795
796 do i = 1, lx * lx * lx
797 grad(1,1,1) = w3(i,1,1) &
798 * ( drdx(i,1,1,e) * ur(i,1,1) &
799 + dsdx(i,1,1,e) * us(i,1,1) &
800 + dtdx(i,1,1,e) * ut(i,1,1) )
801 grad(1,1,2) = w3(i,1,1) &
802 * ( dsdy(i,1,1,e) * us(i,1,1) &
803 + drdy(i,1,1,e) * ur(i,1,1) &
804 + dtdy(i,1,1,e) * ut(i,1,1) )
805 grad(1,1,3) = w3(i,1,1) &
806 * ( dtdz(i,1,1,e) * ut(i,1,1) &
807 + drdz(i,1,1,e) * ur(i,1,1) &
808 + dsdz(i,1,1,e) * us(i,1,1) )
809
810 grad(1,2,1) = w3(i,1,1) &
811 * ( drdx(i,1,1,e) * vr(i,1,1) &
812 + dsdx(i,1,1,e) * vs(i,1,1) &
813 + dtdx(i,1,1,e) * vt(i,1,1) )
814 grad(1,2,2) = w3(i,1,1) &
815 * ( dsdy(i,1,1,e) * vs(i,1,1) &
816 + drdy(i,1,1,e) * vr(i,1,1) &
817 + dtdy(i,1,1,e) * vt(i,1,1) )
818 grad(1,2,3) = w3(i,1,1) &
819 * ( dtdz(i,1,1,e) * vt(i,1,1) &
820 + drdz(i,1,1,e) * vr(i,1,1) &
821 + dsdz(i,1,1,e) * vs(i,1,1) )
822
823 grad(1,3,1) = w3(i,1,1) &
824 * ( drdx(i,1,1,e) * wr(i,1,1) &
825 + dsdx(i,1,1,e) * ws(i,1,1) &
826 + dtdx(i,1,1,e) * wt(i,1,1) )
827 grad(1,3,2) = w3(i,1,1) &
828 * ( dsdy(i,1,1,e) * ws(i,1,1) &
829 + drdy(i,1,1,e) * wr(i,1,1) &
830 + dtdy(i,1,1,e) * wt(i,1,1) )
831 grad(1,3,3) = w3(i,1,1) &
832 * ( dtdz(i,1,1,e) * wt(i,1,1) &
833 + drdz(i,1,1,e) * wr(i,1,1) &
834 + dsdz(i,1,1,e) * ws(i,1,1) )
835 end do
836
837
838 do i = 1, lx * lx * lx
839 s11 = grad(i,1,1)
840 s22 = grad(i,2,2)
841 s33 = grad(i,3,3)
842
843
844 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
845 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
846 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
847
848 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
849 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
850 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
851
852 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
853 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
854 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
855
856 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
857 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
858 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
859
860
861 b = -(a11 + a22 + a33)
862 c = -(a12*a12 + a13*a13 + a23*a23 &
863 - a11 * a22 - a11 * a33 - a22 * a33)
864 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
865 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
866
867
868 q = (3.0 * c - b*b) / 9.0
869 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
870 theta = acos( r / sqrt(-q*q*q) )
871
872 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
873 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
874 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
875
876 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
877 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
878 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
879 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
880 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
881 .le. eigen(2) .and. eigen(2) .le. eigen(1))
882 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
883 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
884 .le. eigen(3) .and. eigen(3) .le. eigen(1))
885
886 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
887
888 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
889 end do
890 end do
891 end subroutine sx_lambda2_lx16
892
893 subroutine sx_lambda2_lx15(lambda2, u, v, w, dx, dy, dz, &
894 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
895 integer, parameter :: lx = 15
896 integer, intent(in) :: n
897 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
898 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
899 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
900 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
901 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
902 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
903 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
904 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
905 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
906 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
907 real(kind=rp) :: grad(lx*lx*lx,3,3)
908 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
909 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
910 real(kind=rp) :: a11, a22, a33, a12, a13, a23
911 real(kind=rp) :: msk1, msk2, msk3
912 real(kind=rp) :: ur(lx, lx, lx)
913 real(kind=rp) :: vr(lx, lx, lx)
914 real(kind=rp) :: wr(lx, lx, lx)
915 real(kind=rp) :: us(lx, lx, lx)
916 real(kind=rp) :: vs(lx, lx, lx)
917 real(kind=rp) :: ws(lx, lx, lx)
918 real(kind=rp) :: ut(lx, lx, lx)
919 real(kind=rp) :: vt(lx, lx, lx)
920 real(kind=rp) :: wt(lx, lx, lx)
921 real(kind=rp) :: tmp1, tmp2, tmp3
922 integer :: e, i, j, k, l
923
924 do e = 1, n
925 do j = 1, lx * lx
926 do i = 1, lx
927 tmp1 = 0.0_rp
928 tmp2 = 0.0_rp
929 tmp3 = 0.0_rp
930 do k = 1, lx
931 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
932 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
933 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
934 end do
935 ur(i,j,1) = tmp1
936 vr(i,j,1) = tmp2
937 wr(i,j,1) = tmp3
938 end do
939 end do
940
941 do k = 1, lx
942 do j = 1, lx
943 do i = 1, lx
944 tmp1 = 0.0_rp
945 tmp2 = 0.0_rp
946 tmp3 = 0.0_rp
947 do l = 1, lx
948 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
949 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
950 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
951 end do
952 us(i,j,k) = tmp1
953 vs(i,j,k) = tmp2
954 ws(i,j,k) = tmp3
955 end do
956 end do
957 end do
958
959 do k = 1, lx
960 do i = 1, lx*lx
961 tmp1 = 0.0_rp
962 tmp2 = 0.0_rp
963 tmp3 = 0.0_rp
964 do l = 1, lx
965 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
966 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
967 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
968 end do
969 ut(i,1,k) = tmp1
970 vt(i,1,k) = tmp2
971 wt(i,1,k) = tmp3
972 end do
973 end do
974
975 do i = 1, lx * lx * lx
976 grad(1,1,1) = w3(i,1,1) &
977 * ( drdx(i,1,1,e) * ur(i,1,1) &
978 + dsdx(i,1,1,e) * us(i,1,1) &
979 + dtdx(i,1,1,e) * ut(i,1,1) )
980 grad(1,1,2) = w3(i,1,1) &
981 * ( dsdy(i,1,1,e) * us(i,1,1) &
982 + drdy(i,1,1,e) * ur(i,1,1) &
983 + dtdy(i,1,1,e) * ut(i,1,1) )
984 grad(1,1,3) = w3(i,1,1) &
985 * ( dtdz(i,1,1,e) * ut(i,1,1) &
986 + drdz(i,1,1,e) * ur(i,1,1) &
987 + dsdz(i,1,1,e) * us(i,1,1) )
988
989 grad(1,2,1) = w3(i,1,1) &
990 * ( drdx(i,1,1,e) * vr(i,1,1) &
991 + dsdx(i,1,1,e) * vs(i,1,1) &
992 + dtdx(i,1,1,e) * vt(i,1,1) )
993 grad(1,2,2) = w3(i,1,1) &
994 * ( dsdy(i,1,1,e) * vs(i,1,1) &
995 + drdy(i,1,1,e) * vr(i,1,1) &
996 + dtdy(i,1,1,e) * vt(i,1,1) )
997 grad(1,2,3) = w3(i,1,1) &
998 * ( dtdz(i,1,1,e) * vt(i,1,1) &
999 + drdz(i,1,1,e) * vr(i,1,1) &
1000 + dsdz(i,1,1,e) * vs(i,1,1) )
1001
1002 grad(1,3,1) = w3(i,1,1) &
1003 * ( drdx(i,1,1,e) * wr(i,1,1) &
1004 + dsdx(i,1,1,e) * ws(i,1,1) &
1005 + dtdx(i,1,1,e) * wt(i,1,1) )
1006 grad(1,3,2) = w3(i,1,1) &
1007 * ( dsdy(i,1,1,e) * ws(i,1,1) &
1008 + drdy(i,1,1,e) * wr(i,1,1) &
1009 + dtdy(i,1,1,e) * wt(i,1,1) )
1010 grad(1,3,3) = w3(i,1,1) &
1011 * ( dtdz(i,1,1,e) * wt(i,1,1) &
1012 + drdz(i,1,1,e) * wr(i,1,1) &
1013 + dsdz(i,1,1,e) * ws(i,1,1) )
1014 end do
1015
1016
1017 do i = 1, lx * lx * lx
1018 s11 = grad(i,1,1)
1019 s22 = grad(i,2,2)
1020 s33 = grad(i,3,3)
1021
1022
1023 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
1024 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
1025 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
1026
1027 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
1028 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
1029 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
1030
1031 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
1032 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
1033 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
1034
1035 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
1036 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
1037 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
1038
1039
1040 b = -(a11 + a22 + a33)
1041 c = -(a12*a12 + a13*a13 + a23*a23 &
1042 - a11 * a22 - a11 * a33 - a22 * a33)
1043 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
1044 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
1045
1046
1047 q = (3.0 * c - b*b) / 9.0
1048 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
1049 theta = acos( r / sqrt(-q*q*q) )
1050
1051 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
1052 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
1053 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
1054
1055 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
1056 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
1057 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
1058 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
1059 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
1060 .le. eigen(2) .and. eigen(2) .le. eigen(1))
1061 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
1062 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
1063 .le. eigen(3) .and. eigen(3) .le. eigen(1))
1064
1065 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1066
1067 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1068 end do
1069 end do
1070 end subroutine sx_lambda2_lx15
1071
1072 subroutine sx_lambda2_lx14(lambda2, u, v, w, dx, dy, dz, &
1073 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
1074 integer, parameter :: lx = 14
1075 integer, intent(in) :: n
1076 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
1077 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1078 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
1079 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
1080 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1081 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1082 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1083 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1084 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
1085 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1086 real(kind=rp) :: grad(lx*lx*lx,3,3)
1087 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
1088 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
1089 real(kind=rp) :: a11, a22, a33, a12, a13, a23
1090 real(kind=rp) :: msk1, msk2, msk3
1091 real(kind=rp) :: ur(lx, lx, lx)
1092 real(kind=rp) :: vr(lx, lx, lx)
1093 real(kind=rp) :: wr(lx, lx, lx)
1094 real(kind=rp) :: us(lx, lx, lx)
1095 real(kind=rp) :: vs(lx, lx, lx)
1096 real(kind=rp) :: ws(lx, lx, lx)
1097 real(kind=rp) :: ut(lx, lx, lx)
1098 real(kind=rp) :: vt(lx, lx, lx)
1099 real(kind=rp) :: wt(lx, lx, lx)
1100 real(kind=rp) :: tmp1, tmp2, tmp3
1101 integer :: e, i, j, k, l
1102
1103 do e = 1, n
1104 do j = 1, lx * lx
1105 do i = 1, lx
1106 tmp1 = 0.0_rp
1107 tmp2 = 0.0_rp
1108 tmp3 = 0.0_rp
1109 do k = 1, lx
1110 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
1111 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
1112 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
1113 end do
1114 ur(i,j,1) = tmp1
1115 vr(i,j,1) = tmp2
1116 wr(i,j,1) = tmp3
1117 end do
1118 end do
1119
1120 do k = 1, lx
1121 do j = 1, lx
1122 do i = 1, lx
1123 tmp1 = 0.0_rp
1124 tmp2 = 0.0_rp
1125 tmp3 = 0.0_rp
1126 do l = 1, lx
1127 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
1128 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
1129 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
1130 end do
1131 us(i,j,k) = tmp1
1132 vs(i,j,k) = tmp2
1133 ws(i,j,k) = tmp3
1134 end do
1135 end do
1136 end do
1137
1138 do k = 1, lx
1139 do i = 1, lx*lx
1140 tmp1 = 0.0_rp
1141 tmp2 = 0.0_rp
1142 tmp3 = 0.0_rp
1143 do l = 1, lx
1144 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
1145 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
1146 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
1147 end do
1148 ut(i,1,k) = tmp1
1149 vt(i,1,k) = tmp2
1150 wt(i,1,k) = tmp3
1151 end do
1152 end do
1153
1154 do i = 1, lx * lx * lx
1155 grad(1,1,1) = w3(i,1,1) &
1156 * ( drdx(i,1,1,e) * ur(i,1,1) &
1157 + dsdx(i,1,1,e) * us(i,1,1) &
1158 + dtdx(i,1,1,e) * ut(i,1,1) )
1159 grad(1,1,2) = w3(i,1,1) &
1160 * ( dsdy(i,1,1,e) * us(i,1,1) &
1161 + drdy(i,1,1,e) * ur(i,1,1) &
1162 + dtdy(i,1,1,e) * ut(i,1,1) )
1163 grad(1,1,3) = w3(i,1,1) &
1164 * ( dtdz(i,1,1,e) * ut(i,1,1) &
1165 + drdz(i,1,1,e) * ur(i,1,1) &
1166 + dsdz(i,1,1,e) * us(i,1,1) )
1167
1168 grad(1,2,1) = w3(i,1,1) &
1169 * ( drdx(i,1,1,e) * vr(i,1,1) &
1170 + dsdx(i,1,1,e) * vs(i,1,1) &
1171 + dtdx(i,1,1,e) * vt(i,1,1) )
1172 grad(1,2,2) = w3(i,1,1) &
1173 * ( dsdy(i,1,1,e) * vs(i,1,1) &
1174 + drdy(i,1,1,e) * vr(i,1,1) &
1175 + dtdy(i,1,1,e) * vt(i,1,1) )
1176 grad(1,2,3) = w3(i,1,1) &
1177 * ( dtdz(i,1,1,e) * vt(i,1,1) &
1178 + drdz(i,1,1,e) * vr(i,1,1) &
1179 + dsdz(i,1,1,e) * vs(i,1,1) )
1180
1181 grad(1,3,1) = w3(i,1,1) &
1182 * ( drdx(i,1,1,e) * wr(i,1,1) &
1183 + dsdx(i,1,1,e) * ws(i,1,1) &
1184 + dtdx(i,1,1,e) * wt(i,1,1) )
1185 grad(1,3,2) = w3(i,1,1) &
1186 * ( dsdy(i,1,1,e) * ws(i,1,1) &
1187 + drdy(i,1,1,e) * wr(i,1,1) &
1188 + dtdy(i,1,1,e) * wt(i,1,1) )
1189 grad(1,3,3) = w3(i,1,1) &
1190 * ( dtdz(i,1,1,e) * wt(i,1,1) &
1191 + drdz(i,1,1,e) * wr(i,1,1) &
1192 + dsdz(i,1,1,e) * ws(i,1,1) )
1193 end do
1194
1195
1196 do i = 1, lx * lx * lx
1197 s11 = grad(i,1,1)
1198 s22 = grad(i,2,2)
1199 s33 = grad(i,3,3)
1200
1201
1202 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
1203 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
1204 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
1205
1206 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
1207 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
1208 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
1209
1210 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
1211 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
1212 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
1213
1214 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
1215 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
1216 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
1217
1218
1219 b = -(a11 + a22 + a33)
1220 c = -(a12*a12 + a13*a13 + a23*a23 &
1221 - a11 * a22 - a11 * a33 - a22 * a33)
1222 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
1223 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
1224
1225
1226 q = (3.0 * c - b*b) / 9.0
1227 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
1228 theta = acos( r / sqrt(-q*q*q) )
1229
1230 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
1231 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
1232 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
1233
1234 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
1235 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
1236 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
1237 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
1238 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
1239 .le. eigen(2) .and. eigen(2) .le. eigen(1))
1240 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
1241 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
1242 .le. eigen(3) .and. eigen(3) .le. eigen(1))
1243
1244 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1245
1246 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1247 end do
1248 end do
1249 end subroutine sx_lambda2_lx14
1250
1251 subroutine sx_lambda2_lx13(lambda2, u, v, w, dx, dy, dz, &
1252 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
1253 integer, parameter :: lx = 13
1254 integer, intent(in) :: n
1255 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
1256 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1257 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
1258 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
1259 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1260 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1261 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1262 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1263 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
1264 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1265 real(kind=rp) :: grad(lx*lx*lx,3,3)
1266 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
1267 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
1268 real(kind=rp) :: a11, a22, a33, a12, a13, a23
1269 real(kind=rp) :: msk1, msk2, msk3
1270 real(kind=rp) :: ur(lx, lx, lx)
1271 real(kind=rp) :: vr(lx, lx, lx)
1272 real(kind=rp) :: wr(lx, lx, lx)
1273 real(kind=rp) :: us(lx, lx, lx)
1274 real(kind=rp) :: vs(lx, lx, lx)
1275 real(kind=rp) :: ws(lx, lx, lx)
1276 real(kind=rp) :: ut(lx, lx, lx)
1277 real(kind=rp) :: vt(lx, lx, lx)
1278 real(kind=rp) :: wt(lx, lx, lx)
1279 real(kind=rp) :: tmp1, tmp2, tmp3
1280 integer :: e, i, j, k, l
1281
1282 do e = 1, n
1283 do j = 1, lx * lx
1284 do i = 1, lx
1285 tmp1 = 0.0_rp
1286 tmp2 = 0.0_rp
1287 tmp3 = 0.0_rp
1288 do k = 1, lx
1289 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
1290 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
1291 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
1292 end do
1293 ur(i,j,1) = tmp1
1294 vr(i,j,1) = tmp2
1295 wr(i,j,1) = tmp3
1296 end do
1297 end do
1298
1299 do k = 1, lx
1300 do j = 1, lx
1301 do i = 1, lx
1302 tmp1 = 0.0_rp
1303 tmp2 = 0.0_rp
1304 tmp3 = 0.0_rp
1305 do l = 1, lx
1306 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
1307 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
1308 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
1309 end do
1310 us(i,j,k) = tmp1
1311 vs(i,j,k) = tmp2
1312 ws(i,j,k) = tmp3
1313 end do
1314 end do
1315 end do
1316
1317 do k = 1, lx
1318 do i = 1, lx*lx
1319 tmp1 = 0.0_rp
1320 tmp2 = 0.0_rp
1321 tmp3 = 0.0_rp
1322 do l = 1, lx
1323 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
1324 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
1325 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
1326 end do
1327 ut(i,1,k) = tmp1
1328 vt(i,1,k) = tmp2
1329 wt(i,1,k) = tmp3
1330 end do
1331 end do
1332
1333 do i = 1, lx * lx * lx
1334 grad(1,1,1) = w3(i,1,1) &
1335 * ( drdx(i,1,1,e) * ur(i,1,1) &
1336 + dsdx(i,1,1,e) * us(i,1,1) &
1337 + dtdx(i,1,1,e) * ut(i,1,1) )
1338 grad(1,1,2) = w3(i,1,1) &
1339 * ( dsdy(i,1,1,e) * us(i,1,1) &
1340 + drdy(i,1,1,e) * ur(i,1,1) &
1341 + dtdy(i,1,1,e) * ut(i,1,1) )
1342 grad(1,1,3) = w3(i,1,1) &
1343 * ( dtdz(i,1,1,e) * ut(i,1,1) &
1344 + drdz(i,1,1,e) * ur(i,1,1) &
1345 + dsdz(i,1,1,e) * us(i,1,1) )
1346
1347 grad(1,2,1) = w3(i,1,1) &
1348 * ( drdx(i,1,1,e) * vr(i,1,1) &
1349 + dsdx(i,1,1,e) * vs(i,1,1) &
1350 + dtdx(i,1,1,e) * vt(i,1,1) )
1351 grad(1,2,2) = w3(i,1,1) &
1352 * ( dsdy(i,1,1,e) * vs(i,1,1) &
1353 + drdy(i,1,1,e) * vr(i,1,1) &
1354 + dtdy(i,1,1,e) * vt(i,1,1) )
1355 grad(1,2,3) = w3(i,1,1) &
1356 * ( dtdz(i,1,1,e) * vt(i,1,1) &
1357 + drdz(i,1,1,e) * vr(i,1,1) &
1358 + dsdz(i,1,1,e) * vs(i,1,1) )
1359
1360 grad(1,3,1) = w3(i,1,1) &
1361 * ( drdx(i,1,1,e) * wr(i,1,1) &
1362 + dsdx(i,1,1,e) * ws(i,1,1) &
1363 + dtdx(i,1,1,e) * wt(i,1,1) )
1364 grad(1,3,2) = w3(i,1,1) &
1365 * ( dsdy(i,1,1,e) * ws(i,1,1) &
1366 + drdy(i,1,1,e) * wr(i,1,1) &
1367 + dtdy(i,1,1,e) * wt(i,1,1) )
1368 grad(1,3,3) = w3(i,1,1) &
1369 * ( dtdz(i,1,1,e) * wt(i,1,1) &
1370 + drdz(i,1,1,e) * wr(i,1,1) &
1371 + dsdz(i,1,1,e) * ws(i,1,1) )
1372 end do
1373
1374
1375 do i = 1, lx * lx * lx
1376 s11 = grad(i,1,1)
1377 s22 = grad(i,2,2)
1378 s33 = grad(i,3,3)
1379
1380
1381 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
1382 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
1383 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
1384
1385 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
1386 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
1387 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
1388
1389 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
1390 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
1391 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
1392
1393 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
1394 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
1395 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
1396
1397
1398 b = -(a11 + a22 + a33)
1399 c = -(a12*a12 + a13*a13 + a23*a23 &
1400 - a11 * a22 - a11 * a33 - a22 * a33)
1401 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
1402 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
1403
1404
1405 q = (3.0 * c - b*b) / 9.0
1406 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
1407 theta = acos( r / sqrt(-q*q*q) )
1408
1409 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
1410 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
1411 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
1412
1413 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
1414 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
1415 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
1416 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
1417 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
1418 .le. eigen(2) .and. eigen(2) .le. eigen(1))
1419 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
1420 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
1421 .le. eigen(3) .and. eigen(3) .le. eigen(1))
1422
1423 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1424
1425 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1426 end do
1427 end do
1428 end subroutine sx_lambda2_lx13
1429
1430 subroutine sx_lambda2_lx12(lambda2, u, v, w, dx, dy, dz, &
1431 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
1432 integer, parameter :: lx = 12
1433 integer, intent(in) :: n
1434 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
1435 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1436 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
1437 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
1438 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1439 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1440 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1441 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1442 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
1443 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1444 real(kind=rp) :: grad(lx*lx*lx,3,3)
1445 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
1446 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
1447 real(kind=rp) :: a11, a22, a33, a12, a13, a23
1448 real(kind=rp) :: msk1, msk2, msk3
1449 real(kind=rp) :: ur(lx, lx, lx)
1450 real(kind=rp) :: vr(lx, lx, lx)
1451 real(kind=rp) :: wr(lx, lx, lx)
1452 real(kind=rp) :: us(lx, lx, lx)
1453 real(kind=rp) :: vs(lx, lx, lx)
1454 real(kind=rp) :: ws(lx, lx, lx)
1455 real(kind=rp) :: ut(lx, lx, lx)
1456 real(kind=rp) :: vt(lx, lx, lx)
1457 real(kind=rp) :: wt(lx, lx, lx)
1458 real(kind=rp) :: tmp1, tmp2, tmp3
1459 integer :: e, i, j, k, l
1460
1461 do e = 1, n
1462 do j = 1, lx * lx
1463 do i = 1, lx
1464 tmp1 = 0.0_rp
1465 tmp2 = 0.0_rp
1466 tmp3 = 0.0_rp
1467 do k = 1, lx
1468 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
1469 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
1470 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
1471 end do
1472 ur(i,j,1) = tmp1
1473 vr(i,j,1) = tmp2
1474 wr(i,j,1) = tmp3
1475 end do
1476 end do
1477
1478 do k = 1, lx
1479 do j = 1, lx
1480 do i = 1, lx
1481 tmp1 = 0.0_rp
1482 tmp2 = 0.0_rp
1483 tmp3 = 0.0_rp
1484 do l = 1, lx
1485 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
1486 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
1487 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
1488 end do
1489 us(i,j,k) = tmp1
1490 vs(i,j,k) = tmp2
1491 ws(i,j,k) = tmp3
1492 end do
1493 end do
1494 end do
1495
1496 do k = 1, lx
1497 do i = 1, lx*lx
1498 tmp1 = 0.0_rp
1499 tmp2 = 0.0_rp
1500 tmp3 = 0.0_rp
1501 do l = 1, lx
1502 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
1503 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
1504 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
1505 end do
1506 ut(i,1,k) = tmp1
1507 vt(i,1,k) = tmp2
1508 wt(i,1,k) = tmp3
1509 end do
1510 end do
1511
1512 do i = 1, lx * lx * lx
1513 grad(1,1,1) = w3(i,1,1) &
1514 * ( drdx(i,1,1,e) * ur(i,1,1) &
1515 + dsdx(i,1,1,e) * us(i,1,1) &
1516 + dtdx(i,1,1,e) * ut(i,1,1) )
1517 grad(1,1,2) = w3(i,1,1) &
1518 * ( dsdy(i,1,1,e) * us(i,1,1) &
1519 + drdy(i,1,1,e) * ur(i,1,1) &
1520 + dtdy(i,1,1,e) * ut(i,1,1) )
1521 grad(1,1,3) = w3(i,1,1) &
1522 * ( dtdz(i,1,1,e) * ut(i,1,1) &
1523 + drdz(i,1,1,e) * ur(i,1,1) &
1524 + dsdz(i,1,1,e) * us(i,1,1) )
1525
1526 grad(1,2,1) = w3(i,1,1) &
1527 * ( drdx(i,1,1,e) * vr(i,1,1) &
1528 + dsdx(i,1,1,e) * vs(i,1,1) &
1529 + dtdx(i,1,1,e) * vt(i,1,1) )
1530 grad(1,2,2) = w3(i,1,1) &
1531 * ( dsdy(i,1,1,e) * vs(i,1,1) &
1532 + drdy(i,1,1,e) * vr(i,1,1) &
1533 + dtdy(i,1,1,e) * vt(i,1,1) )
1534 grad(1,2,3) = w3(i,1,1) &
1535 * ( dtdz(i,1,1,e) * vt(i,1,1) &
1536 + drdz(i,1,1,e) * vr(i,1,1) &
1537 + dsdz(i,1,1,e) * vs(i,1,1) )
1538
1539 grad(1,3,1) = w3(i,1,1) &
1540 * ( drdx(i,1,1,e) * wr(i,1,1) &
1541 + dsdx(i,1,1,e) * ws(i,1,1) &
1542 + dtdx(i,1,1,e) * wt(i,1,1) )
1543 grad(1,3,2) = w3(i,1,1) &
1544 * ( dsdy(i,1,1,e) * ws(i,1,1) &
1545 + drdy(i,1,1,e) * wr(i,1,1) &
1546 + dtdy(i,1,1,e) * wt(i,1,1) )
1547 grad(1,3,3) = w3(i,1,1) &
1548 * ( dtdz(i,1,1,e) * wt(i,1,1) &
1549 + drdz(i,1,1,e) * wr(i,1,1) &
1550 + dsdz(i,1,1,e) * ws(i,1,1) )
1551 end do
1552
1553
1554 do i = 1, lx * lx * lx
1555 s11 = grad(i,1,1)
1556 s22 = grad(i,2,2)
1557 s33 = grad(i,3,3)
1558
1559
1560 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
1561 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
1562 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
1563
1564 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
1565 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
1566 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
1567
1568 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
1569 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
1570 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
1571
1572 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
1573 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
1574 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
1575
1576
1577 b = -(a11 + a22 + a33)
1578 c = -(a12*a12 + a13*a13 + a23*a23 &
1579 - a11 * a22 - a11 * a33 - a22 * a33)
1580 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
1581 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
1582
1583
1584 q = (3.0 * c - b*b) / 9.0
1585 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
1586 theta = acos( r / sqrt(-q*q*q) )
1587
1588 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
1589 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
1590 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
1591
1592 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
1593 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
1594 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
1595 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
1596 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
1597 .le. eigen(2) .and. eigen(2) .le. eigen(1))
1598 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
1599 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
1600 .le. eigen(3) .and. eigen(3) .le. eigen(1))
1601
1602 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1603
1604 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1605 end do
1606 end do
1607 end subroutine sx_lambda2_lx12
1608
1609 subroutine sx_lambda2_lx11(lambda2, u, v, w, dx, dy, dz, &
1610 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
1611 integer, parameter :: lx = 11
1612 integer, intent(in) :: n
1613 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
1614 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1615 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
1616 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
1617 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1618 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1619 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1620 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1621 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
1622 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1623 real(kind=rp) :: grad(lx*lx*lx,3,3)
1624 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
1625 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
1626 real(kind=rp) :: a11, a22, a33, a12, a13, a23
1627 real(kind=rp) :: msk1, msk2, msk3
1628 real(kind=rp) :: ur(lx, lx, lx)
1629 real(kind=rp) :: vr(lx, lx, lx)
1630 real(kind=rp) :: wr(lx, lx, lx)
1631 real(kind=rp) :: us(lx, lx, lx)
1632 real(kind=rp) :: vs(lx, lx, lx)
1633 real(kind=rp) :: ws(lx, lx, lx)
1634 real(kind=rp) :: ut(lx, lx, lx)
1635 real(kind=rp) :: vt(lx, lx, lx)
1636 real(kind=rp) :: wt(lx, lx, lx)
1637 real(kind=rp) :: tmp1, tmp2, tmp3
1638 integer :: e, i, j, k, l
1639
1640 do e = 1, n
1641 do j = 1, lx * lx
1642 do i = 1, lx
1643 tmp1 = 0.0_rp
1644 tmp2 = 0.0_rp
1645 tmp3 = 0.0_rp
1646 do k = 1, lx
1647 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
1648 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
1649 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
1650 end do
1651 ur(i,j,1) = tmp1
1652 vr(i,j,1) = tmp2
1653 wr(i,j,1) = tmp3
1654 end do
1655 end do
1656
1657 do k = 1, lx
1658 do j = 1, lx
1659 do i = 1, lx
1660 tmp1 = 0.0_rp
1661 tmp2 = 0.0_rp
1662 tmp3 = 0.0_rp
1663 do l = 1, lx
1664 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
1665 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
1666 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
1667 end do
1668 us(i,j,k) = tmp1
1669 vs(i,j,k) = tmp2
1670 ws(i,j,k) = tmp3
1671 end do
1672 end do
1673 end do
1674
1675 do k = 1, lx
1676 do i = 1, lx*lx
1677 tmp1 = 0.0_rp
1678 tmp2 = 0.0_rp
1679 tmp3 = 0.0_rp
1680 do l = 1, lx
1681 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
1682 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
1683 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
1684 end do
1685 ut(i,1,k) = tmp1
1686 vt(i,1,k) = tmp2
1687 wt(i,1,k) = tmp3
1688 end do
1689 end do
1690
1691 do i = 1, lx * lx * lx
1692 grad(1,1,1) = w3(i,1,1) &
1693 * ( drdx(i,1,1,e) * ur(i,1,1) &
1694 + dsdx(i,1,1,e) * us(i,1,1) &
1695 + dtdx(i,1,1,e) * ut(i,1,1) )
1696 grad(1,1,2) = w3(i,1,1) &
1697 * ( dsdy(i,1,1,e) * us(i,1,1) &
1698 + drdy(i,1,1,e) * ur(i,1,1) &
1699 + dtdy(i,1,1,e) * ut(i,1,1) )
1700 grad(1,1,3) = w3(i,1,1) &
1701 * ( dtdz(i,1,1,e) * ut(i,1,1) &
1702 + drdz(i,1,1,e) * ur(i,1,1) &
1703 + dsdz(i,1,1,e) * us(i,1,1) )
1704
1705 grad(1,2,1) = w3(i,1,1) &
1706 * ( drdx(i,1,1,e) * vr(i,1,1) &
1707 + dsdx(i,1,1,e) * vs(i,1,1) &
1708 + dtdx(i,1,1,e) * vt(i,1,1) )
1709 grad(1,2,2) = w3(i,1,1) &
1710 * ( dsdy(i,1,1,e) * vs(i,1,1) &
1711 + drdy(i,1,1,e) * vr(i,1,1) &
1712 + dtdy(i,1,1,e) * vt(i,1,1) )
1713 grad(1,2,3) = w3(i,1,1) &
1714 * ( dtdz(i,1,1,e) * vt(i,1,1) &
1715 + drdz(i,1,1,e) * vr(i,1,1) &
1716 + dsdz(i,1,1,e) * vs(i,1,1) )
1717
1718 grad(1,3,1) = w3(i,1,1) &
1719 * ( drdx(i,1,1,e) * wr(i,1,1) &
1720 + dsdx(i,1,1,e) * ws(i,1,1) &
1721 + dtdx(i,1,1,e) * wt(i,1,1) )
1722 grad(1,3,2) = w3(i,1,1) &
1723 * ( dsdy(i,1,1,e) * ws(i,1,1) &
1724 + drdy(i,1,1,e) * wr(i,1,1) &
1725 + dtdy(i,1,1,e) * wt(i,1,1) )
1726 grad(1,3,3) = w3(i,1,1) &
1727 * ( dtdz(i,1,1,e) * wt(i,1,1) &
1728 + drdz(i,1,1,e) * wr(i,1,1) &
1729 + dsdz(i,1,1,e) * ws(i,1,1) )
1730 end do
1731
1732
1733 do i = 1, lx * lx * lx
1734 s11 = grad(i,1,1)
1735 s22 = grad(i,2,2)
1736 s33 = grad(i,3,3)
1737
1738
1739 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
1740 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
1741 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
1742
1743 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
1744 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
1745 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
1746
1747 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
1748 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
1749 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
1750
1751 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
1752 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
1753 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
1754
1755
1756 b = -(a11 + a22 + a33)
1757 c = -(a12*a12 + a13*a13 + a23*a23 &
1758 - a11 * a22 - a11 * a33 - a22 * a33)
1759 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
1760 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
1761
1762
1763 q = (3.0 * c - b*b) / 9.0
1764 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
1765 theta = acos( r / sqrt(-q*q*q) )
1766
1767 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
1768 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
1769 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
1770
1771 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
1772 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
1773 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
1774 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
1775 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
1776 .le. eigen(2) .and. eigen(2) .le. eigen(1))
1777 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
1778 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
1779 .le. eigen(3) .and. eigen(3) .le. eigen(1))
1780
1781 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1782
1783 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1784 end do
1785 end do
1786 end subroutine sx_lambda2_lx11
1787
1788 subroutine sx_lambda2_lx10(lambda2, u, v, w, dx, dy, dz, &
1789 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
1790 integer, parameter :: lx = 10
1791 integer, intent(in) :: n
1792 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
1793 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1794 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
1795 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
1796 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1797 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1798 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1799 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1800 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
1801 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1802 real(kind=rp) :: grad(lx*lx*lx,3,3)
1803 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
1804 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
1805 real(kind=rp) :: a11, a22, a33, a12, a13, a23
1806 real(kind=rp) :: msk1, msk2, msk3
1807 real(kind=rp) :: ur(lx, lx, lx)
1808 real(kind=rp) :: vr(lx, lx, lx)
1809 real(kind=rp) :: wr(lx, lx, lx)
1810 real(kind=rp) :: us(lx, lx, lx)
1811 real(kind=rp) :: vs(lx, lx, lx)
1812 real(kind=rp) :: ws(lx, lx, lx)
1813 real(kind=rp) :: ut(lx, lx, lx)
1814 real(kind=rp) :: vt(lx, lx, lx)
1815 real(kind=rp) :: wt(lx, lx, lx)
1816 real(kind=rp) :: tmp1, tmp2, tmp3
1817 integer :: e, i, j, k, l
1818
1819 do e = 1, n
1820 do j = 1, lx * lx
1821 do i = 1, lx
1822 tmp1 = 0.0_rp
1823 tmp2 = 0.0_rp
1824 tmp3 = 0.0_rp
1825 do k = 1, lx
1826 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
1827 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
1828 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
1829 end do
1830 ur(i,j,1) = tmp1
1831 vr(i,j,1) = tmp2
1832 wr(i,j,1) = tmp3
1833 end do
1834 end do
1835
1836 do k = 1, lx
1837 do j = 1, lx
1838 do i = 1, lx
1839 tmp1 = 0.0_rp
1840 tmp2 = 0.0_rp
1841 tmp3 = 0.0_rp
1842 do l = 1, lx
1843 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
1844 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
1845 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
1846 end do
1847 us(i,j,k) = tmp1
1848 vs(i,j,k) = tmp2
1849 ws(i,j,k) = tmp3
1850 end do
1851 end do
1852 end do
1853
1854 do k = 1, lx
1855 do i = 1, lx*lx
1856 tmp1 = 0.0_rp
1857 tmp2 = 0.0_rp
1858 tmp3 = 0.0_rp
1859 do l = 1, lx
1860 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
1861 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
1862 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
1863 end do
1864 ut(i,1,k) = tmp1
1865 vt(i,1,k) = tmp2
1866 wt(i,1,k) = tmp3
1867 end do
1868 end do
1869
1870 do i = 1, lx * lx * lx
1871 grad(1,1,1) = w3(i,1,1) &
1872 * ( drdx(i,1,1,e) * ur(i,1,1) &
1873 + dsdx(i,1,1,e) * us(i,1,1) &
1874 + dtdx(i,1,1,e) * ut(i,1,1) )
1875 grad(1,1,2) = w3(i,1,1) &
1876 * ( dsdy(i,1,1,e) * us(i,1,1) &
1877 + drdy(i,1,1,e) * ur(i,1,1) &
1878 + dtdy(i,1,1,e) * ut(i,1,1) )
1879 grad(1,1,3) = w3(i,1,1) &
1880 * ( dtdz(i,1,1,e) * ut(i,1,1) &
1881 + drdz(i,1,1,e) * ur(i,1,1) &
1882 + dsdz(i,1,1,e) * us(i,1,1) )
1883
1884 grad(1,2,1) = w3(i,1,1) &
1885 * ( drdx(i,1,1,e) * vr(i,1,1) &
1886 + dsdx(i,1,1,e) * vs(i,1,1) &
1887 + dtdx(i,1,1,e) * vt(i,1,1) )
1888 grad(1,2,2) = w3(i,1,1) &
1889 * ( dsdy(i,1,1,e) * vs(i,1,1) &
1890 + drdy(i,1,1,e) * vr(i,1,1) &
1891 + dtdy(i,1,1,e) * vt(i,1,1) )
1892 grad(1,2,3) = w3(i,1,1) &
1893 * ( dtdz(i,1,1,e) * vt(i,1,1) &
1894 + drdz(i,1,1,e) * vr(i,1,1) &
1895 + dsdz(i,1,1,e) * vs(i,1,1) )
1896
1897 grad(1,3,1) = w3(i,1,1) &
1898 * ( drdx(i,1,1,e) * wr(i,1,1) &
1899 + dsdx(i,1,1,e) * ws(i,1,1) &
1900 + dtdx(i,1,1,e) * wt(i,1,1) )
1901 grad(1,3,2) = w3(i,1,1) &
1902 * ( dsdy(i,1,1,e) * ws(i,1,1) &
1903 + drdy(i,1,1,e) * wr(i,1,1) &
1904 + dtdy(i,1,1,e) * wt(i,1,1) )
1905 grad(1,3,3) = w3(i,1,1) &
1906 * ( dtdz(i,1,1,e) * wt(i,1,1) &
1907 + drdz(i,1,1,e) * wr(i,1,1) &
1908 + dsdz(i,1,1,e) * ws(i,1,1) )
1909 end do
1910
1911
1912 do i = 1, lx * lx * lx
1913 s11 = grad(i,1,1)
1914 s22 = grad(i,2,2)
1915 s33 = grad(i,3,3)
1916
1917
1918 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
1919 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
1920 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
1921
1922 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
1923 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
1924 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
1925
1926 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
1927 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
1928 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
1929
1930 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
1931 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
1932 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
1933
1934
1935 b = -(a11 + a22 + a33)
1936 c = -(a12*a12 + a13*a13 + a23*a23 &
1937 - a11 * a22 - a11 * a33 - a22 * a33)
1938 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
1939 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
1940
1941
1942 q = (3.0 * c - b*b) / 9.0
1943 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
1944 theta = acos( r / sqrt(-q*q*q) )
1945
1946 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
1947 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
1948 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
1949
1950 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
1951 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
1952 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
1953 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
1954 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
1955 .le. eigen(2) .and. eigen(2) .le. eigen(1))
1956 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
1957 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
1958 .le. eigen(3) .and. eigen(3) .le. eigen(1))
1959
1960 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1961
1962 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1963 end do
1964 end do
1965 end subroutine sx_lambda2_lx10
1966
1967 subroutine sx_lambda2_lx9(lambda2, u, v, w, dx, dy, dz, &
1968 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
1969 integer, parameter :: lx = 9
1970 integer, intent(in) :: n
1971 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
1972 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1973 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
1974 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
1975 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1976 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1977 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1978 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1979 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
1980 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1981 real(kind=rp) :: grad(lx*lx*lx,3,3)
1982 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
1983 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
1984 real(kind=rp) :: a11, a22, a33, a12, a13, a23
1985 real(kind=rp) :: msk1, msk2, msk3
1986 real(kind=rp) :: ur(lx, lx, lx)
1987 real(kind=rp) :: vr(lx, lx, lx)
1988 real(kind=rp) :: wr(lx, lx, lx)
1989 real(kind=rp) :: us(lx, lx, lx)
1990 real(kind=rp) :: vs(lx, lx, lx)
1991 real(kind=rp) :: ws(lx, lx, lx)
1992 real(kind=rp) :: ut(lx, lx, lx)
1993 real(kind=rp) :: vt(lx, lx, lx)
1994 real(kind=rp) :: wt(lx, lx, lx)
1995 real(kind=rp) :: tmp1, tmp2, tmp3
1996 integer :: e, i, j, k, l
1997
1998 do e = 1, n
1999 do j = 1, lx * lx
2000 do i = 1, lx
2001 tmp1 = 0.0_rp
2002 tmp2 = 0.0_rp
2003 tmp3 = 0.0_rp
2004 do k = 1, lx
2005 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
2006 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
2007 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
2008 end do
2009 ur(i,j,1) = tmp1
2010 vr(i,j,1) = tmp2
2011 wr(i,j,1) = tmp3
2012 end do
2013 end do
2014
2015 do k = 1, lx
2016 do j = 1, lx
2017 do i = 1, lx
2018 tmp1 = 0.0_rp
2019 tmp2 = 0.0_rp
2020 tmp3 = 0.0_rp
2021 do l = 1, lx
2022 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
2023 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
2024 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
2025 end do
2026 us(i,j,k) = tmp1
2027 vs(i,j,k) = tmp2
2028 ws(i,j,k) = tmp3
2029 end do
2030 end do
2031 end do
2032
2033 do k = 1, lx
2034 do i = 1, lx*lx
2035 tmp1 = 0.0_rp
2036 tmp2 = 0.0_rp
2037 tmp3 = 0.0_rp
2038 do l = 1, lx
2039 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
2040 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
2041 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
2042 end do
2043 ut(i,1,k) = tmp1
2044 vt(i,1,k) = tmp2
2045 wt(i,1,k) = tmp3
2046 end do
2047 end do
2048
2049 do i = 1, lx * lx * lx
2050 grad(1,1,1) = w3(i,1,1) &
2051 * ( drdx(i,1,1,e) * ur(i,1,1) &
2052 + dsdx(i,1,1,e) * us(i,1,1) &
2053 + dtdx(i,1,1,e) * ut(i,1,1) )
2054 grad(1,1,2) = w3(i,1,1) &
2055 * ( dsdy(i,1,1,e) * us(i,1,1) &
2056 + drdy(i,1,1,e) * ur(i,1,1) &
2057 + dtdy(i,1,1,e) * ut(i,1,1) )
2058 grad(1,1,3) = w3(i,1,1) &
2059 * ( dtdz(i,1,1,e) * ut(i,1,1) &
2060 + drdz(i,1,1,e) * ur(i,1,1) &
2061 + dsdz(i,1,1,e) * us(i,1,1) )
2062
2063 grad(1,2,1) = w3(i,1,1) &
2064 * ( drdx(i,1,1,e) * vr(i,1,1) &
2065 + dsdx(i,1,1,e) * vs(i,1,1) &
2066 + dtdx(i,1,1,e) * vt(i,1,1) )
2067 grad(1,2,2) = w3(i,1,1) &
2068 * ( dsdy(i,1,1,e) * vs(i,1,1) &
2069 + drdy(i,1,1,e) * vr(i,1,1) &
2070 + dtdy(i,1,1,e) * vt(i,1,1) )
2071 grad(1,2,3) = w3(i,1,1) &
2072 * ( dtdz(i,1,1,e) * vt(i,1,1) &
2073 + drdz(i,1,1,e) * vr(i,1,1) &
2074 + dsdz(i,1,1,e) * vs(i,1,1) )
2075
2076 grad(1,3,1) = w3(i,1,1) &
2077 * ( drdx(i,1,1,e) * wr(i,1,1) &
2078 + dsdx(i,1,1,e) * ws(i,1,1) &
2079 + dtdx(i,1,1,e) * wt(i,1,1) )
2080 grad(1,3,2) = w3(i,1,1) &
2081 * ( dsdy(i,1,1,e) * ws(i,1,1) &
2082 + drdy(i,1,1,e) * wr(i,1,1) &
2083 + dtdy(i,1,1,e) * wt(i,1,1) )
2084 grad(1,3,3) = w3(i,1,1) &
2085 * ( dtdz(i,1,1,e) * wt(i,1,1) &
2086 + drdz(i,1,1,e) * wr(i,1,1) &
2087 + dsdz(i,1,1,e) * ws(i,1,1) )
2088 end do
2089
2090
2091 do i = 1, lx * lx * lx
2092 s11 = grad(i,1,1)
2093 s22 = grad(i,2,2)
2094 s33 = grad(i,3,3)
2095
2096
2097 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
2098 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
2099 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
2100
2101 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
2102 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
2103 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
2104
2105 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
2106 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
2107 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
2108
2109 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
2110 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
2111 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
2112
2113
2114 b = -(a11 + a22 + a33)
2115 c = -(a12*a12 + a13*a13 + a23*a23 &
2116 - a11 * a22 - a11 * a33 - a22 * a33)
2117 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
2118 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
2119
2120
2121 q = (3.0 * c - b*b) / 9.0
2122 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
2123 theta = acos( r / sqrt(-q*q*q) )
2124
2125 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
2126 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
2127 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
2128
2129 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
2130 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
2131 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
2132 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
2133 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
2134 .le. eigen(2) .and. eigen(2) .le. eigen(1))
2135 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
2136 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
2137 .le. eigen(3) .and. eigen(3) .le. eigen(1))
2138
2139 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2140
2141 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2142 end do
2143 end do
2144 end subroutine sx_lambda2_lx9
2145
2146 subroutine sx_lambda2_lx8(lambda2, u, v, w, dx, dy, dz, &
2147 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
2148 integer, parameter :: lx = 8
2149 integer, intent(in) :: n
2150 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
2151 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
2152 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
2153 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
2154 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2155 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
2156 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
2157 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
2158 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
2159 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2160 real(kind=rp) :: grad(lx*lx*lx,3,3)
2161 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
2162 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
2163 real(kind=rp) :: a11, a22, a33, a12, a13, a23
2164 real(kind=rp) :: msk1, msk2, msk3
2165 real(kind=rp) :: ur(lx, lx, lx)
2166 real(kind=rp) :: vr(lx, lx, lx)
2167 real(kind=rp) :: wr(lx, lx, lx)
2168 real(kind=rp) :: us(lx, lx, lx)
2169 real(kind=rp) :: vs(lx, lx, lx)
2170 real(kind=rp) :: ws(lx, lx, lx)
2171 real(kind=rp) :: ut(lx, lx, lx)
2172 real(kind=rp) :: vt(lx, lx, lx)
2173 real(kind=rp) :: wt(lx, lx, lx)
2174 real(kind=rp) :: tmp1, tmp2, tmp3
2175 integer :: e, i, j, k, l
2176
2177 do e = 1, n
2178 do j = 1, lx * lx
2179 do i = 1, lx
2180 tmp1 = 0.0_rp
2181 tmp2 = 0.0_rp
2182 tmp3 = 0.0_rp
2183 do k = 1, lx
2184 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
2185 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
2186 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
2187 end do
2188 ur(i,j,1) = tmp1
2189 vr(i,j,1) = tmp2
2190 wr(i,j,1) = tmp3
2191 end do
2192 end do
2193
2194 do k = 1, lx
2195 do j = 1, lx
2196 do i = 1, lx
2197 tmp1 = 0.0_rp
2198 tmp2 = 0.0_rp
2199 tmp3 = 0.0_rp
2200 do l = 1, lx
2201 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
2202 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
2203 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
2204 end do
2205 us(i,j,k) = tmp1
2206 vs(i,j,k) = tmp2
2207 ws(i,j,k) = tmp3
2208 end do
2209 end do
2210 end do
2211
2212 do k = 1, lx
2213 do i = 1, lx*lx
2214 tmp1 = 0.0_rp
2215 tmp2 = 0.0_rp
2216 tmp3 = 0.0_rp
2217 do l = 1, lx
2218 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
2219 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
2220 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
2221 end do
2222 ut(i,1,k) = tmp1
2223 vt(i,1,k) = tmp2
2224 wt(i,1,k) = tmp3
2225 end do
2226 end do
2227
2228 do i = 1, lx * lx * lx
2229 grad(1,1,1) = w3(i,1,1) &
2230 * ( drdx(i,1,1,e) * ur(i,1,1) &
2231 + dsdx(i,1,1,e) * us(i,1,1) &
2232 + dtdx(i,1,1,e) * ut(i,1,1) )
2233 grad(1,1,2) = w3(i,1,1) &
2234 * ( dsdy(i,1,1,e) * us(i,1,1) &
2235 + drdy(i,1,1,e) * ur(i,1,1) &
2236 + dtdy(i,1,1,e) * ut(i,1,1) )
2237 grad(1,1,3) = w3(i,1,1) &
2238 * ( dtdz(i,1,1,e) * ut(i,1,1) &
2239 + drdz(i,1,1,e) * ur(i,1,1) &
2240 + dsdz(i,1,1,e) * us(i,1,1) )
2241
2242 grad(1,2,1) = w3(i,1,1) &
2243 * ( drdx(i,1,1,e) * vr(i,1,1) &
2244 + dsdx(i,1,1,e) * vs(i,1,1) &
2245 + dtdx(i,1,1,e) * vt(i,1,1) )
2246 grad(1,2,2) = w3(i,1,1) &
2247 * ( dsdy(i,1,1,e) * vs(i,1,1) &
2248 + drdy(i,1,1,e) * vr(i,1,1) &
2249 + dtdy(i,1,1,e) * vt(i,1,1) )
2250 grad(1,2,3) = w3(i,1,1) &
2251 * ( dtdz(i,1,1,e) * vt(i,1,1) &
2252 + drdz(i,1,1,e) * vr(i,1,1) &
2253 + dsdz(i,1,1,e) * vs(i,1,1) )
2254
2255 grad(1,3,1) = w3(i,1,1) &
2256 * ( drdx(i,1,1,e) * wr(i,1,1) &
2257 + dsdx(i,1,1,e) * ws(i,1,1) &
2258 + dtdx(i,1,1,e) * wt(i,1,1) )
2259 grad(1,3,2) = w3(i,1,1) &
2260 * ( dsdy(i,1,1,e) * ws(i,1,1) &
2261 + drdy(i,1,1,e) * wr(i,1,1) &
2262 + dtdy(i,1,1,e) * wt(i,1,1) )
2263 grad(1,3,3) = w3(i,1,1) &
2264 * ( dtdz(i,1,1,e) * wt(i,1,1) &
2265 + drdz(i,1,1,e) * wr(i,1,1) &
2266 + dsdz(i,1,1,e) * ws(i,1,1) )
2267 end do
2268
2269
2270 do i = 1, lx * lx * lx
2271 s11 = grad(i,1,1)
2272 s22 = grad(i,2,2)
2273 s33 = grad(i,3,3)
2274
2275
2276 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
2277 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
2278 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
2279
2280 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
2281 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
2282 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
2283
2284 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
2285 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
2286 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
2287
2288 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
2289 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
2290 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
2291
2292
2293 b = -(a11 + a22 + a33)
2294 c = -(a12*a12 + a13*a13 + a23*a23 &
2295 - a11 * a22 - a11 * a33 - a22 * a33)
2296 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
2297 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
2298
2299
2300 q = (3.0 * c - b*b) / 9.0
2301 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
2302 theta = acos( r / sqrt(-q*q*q) )
2303
2304 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
2305 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
2306 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
2307
2308 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
2309 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
2310 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
2311 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
2312 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
2313 .le. eigen(2) .and. eigen(2) .le. eigen(1))
2314 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
2315 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
2316 .le. eigen(3) .and. eigen(3) .le. eigen(1))
2317
2318 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2319
2320 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2321 end do
2322 end do
2323 end subroutine sx_lambda2_lx8
2324
2325 subroutine sx_lambda2_lx7(lambda2, u, v, w, dx, dy, dz, &
2326 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
2327 integer, parameter :: lx = 7
2328 integer, intent(in) :: n
2329 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
2330 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
2331 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
2332 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
2333 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2334 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
2335 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
2336 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
2337 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
2338 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2339 real(kind=rp) :: grad(lx*lx*lx,3,3)
2340 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
2341 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
2342 real(kind=rp) :: a11, a22, a33, a12, a13, a23
2343 real(kind=rp) :: msk1, msk2, msk3
2344 real(kind=rp) :: ur(lx, lx, lx)
2345 real(kind=rp) :: vr(lx, lx, lx)
2346 real(kind=rp) :: wr(lx, lx, lx)
2347 real(kind=rp) :: us(lx, lx, lx)
2348 real(kind=rp) :: vs(lx, lx, lx)
2349 real(kind=rp) :: ws(lx, lx, lx)
2350 real(kind=rp) :: ut(lx, lx, lx)
2351 real(kind=rp) :: vt(lx, lx, lx)
2352 real(kind=rp) :: wt(lx, lx, lx)
2353 real(kind=rp) :: tmp1, tmp2, tmp3
2354 integer :: e, i, j, k, l
2355
2356 do e = 1, n
2357 do j = 1, lx * lx
2358 do i = 1, lx
2359 tmp1 = 0.0_rp
2360 tmp2 = 0.0_rp
2361 tmp3 = 0.0_rp
2362 do k = 1, lx
2363 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
2364 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
2365 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
2366 end do
2367 ur(i,j,1) = tmp1
2368 vr(i,j,1) = tmp2
2369 wr(i,j,1) = tmp3
2370 end do
2371 end do
2372
2373 do k = 1, lx
2374 do j = 1, lx
2375 do i = 1, lx
2376 tmp1 = 0.0_rp
2377 tmp2 = 0.0_rp
2378 tmp3 = 0.0_rp
2379 do l = 1, lx
2380 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
2381 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
2382 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
2383 end do
2384 us(i,j,k) = tmp1
2385 vs(i,j,k) = tmp2
2386 ws(i,j,k) = tmp3
2387 end do
2388 end do
2389 end do
2390
2391 do k = 1, lx
2392 do i = 1, lx*lx
2393 tmp1 = 0.0_rp
2394 tmp2 = 0.0_rp
2395 tmp3 = 0.0_rp
2396 do l = 1, lx
2397 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
2398 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
2399 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
2400 end do
2401 ut(i,1,k) = tmp1
2402 vt(i,1,k) = tmp2
2403 wt(i,1,k) = tmp3
2404 end do
2405 end do
2406
2407 do i = 1, lx * lx * lx
2408 grad(1,1,1) = w3(i,1,1) &
2409 * ( drdx(i,1,1,e) * ur(i,1,1) &
2410 + dsdx(i,1,1,e) * us(i,1,1) &
2411 + dtdx(i,1,1,e) * ut(i,1,1) )
2412 grad(1,1,2) = w3(i,1,1) &
2413 * ( dsdy(i,1,1,e) * us(i,1,1) &
2414 + drdy(i,1,1,e) * ur(i,1,1) &
2415 + dtdy(i,1,1,e) * ut(i,1,1) )
2416 grad(1,1,3) = w3(i,1,1) &
2417 * ( dtdz(i,1,1,e) * ut(i,1,1) &
2418 + drdz(i,1,1,e) * ur(i,1,1) &
2419 + dsdz(i,1,1,e) * us(i,1,1) )
2420
2421 grad(1,2,1) = w3(i,1,1) &
2422 * ( drdx(i,1,1,e) * vr(i,1,1) &
2423 + dsdx(i,1,1,e) * vs(i,1,1) &
2424 + dtdx(i,1,1,e) * vt(i,1,1) )
2425 grad(1,2,2) = w3(i,1,1) &
2426 * ( dsdy(i,1,1,e) * vs(i,1,1) &
2427 + drdy(i,1,1,e) * vr(i,1,1) &
2428 + dtdy(i,1,1,e) * vt(i,1,1) )
2429 grad(1,2,3) = w3(i,1,1) &
2430 * ( dtdz(i,1,1,e) * vt(i,1,1) &
2431 + drdz(i,1,1,e) * vr(i,1,1) &
2432 + dsdz(i,1,1,e) * vs(i,1,1) )
2433
2434 grad(1,3,1) = w3(i,1,1) &
2435 * ( drdx(i,1,1,e) * wr(i,1,1) &
2436 + dsdx(i,1,1,e) * ws(i,1,1) &
2437 + dtdx(i,1,1,e) * wt(i,1,1) )
2438 grad(1,3,2) = w3(i,1,1) &
2439 * ( dsdy(i,1,1,e) * ws(i,1,1) &
2440 + drdy(i,1,1,e) * wr(i,1,1) &
2441 + dtdy(i,1,1,e) * wt(i,1,1) )
2442 grad(1,3,3) = w3(i,1,1) &
2443 * ( dtdz(i,1,1,e) * wt(i,1,1) &
2444 + drdz(i,1,1,e) * wr(i,1,1) &
2445 + dsdz(i,1,1,e) * ws(i,1,1) )
2446 end do
2447
2448
2449 do i = 1, lx * lx * lx
2450 s11 = grad(i,1,1)
2451 s22 = grad(i,2,2)
2452 s33 = grad(i,3,3)
2453
2454
2455 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
2456 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
2457 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
2458
2459 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
2460 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
2461 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
2462
2463 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
2464 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
2465 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
2466
2467 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
2468 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
2469 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
2470
2471
2472 b = -(a11 + a22 + a33)
2473 c = -(a12*a12 + a13*a13 + a23*a23 &
2474 - a11 * a22 - a11 * a33 - a22 * a33)
2475 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
2476 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
2477
2478
2479 q = (3.0 * c - b*b) / 9.0
2480 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
2481 theta = acos( r / sqrt(-q*q*q) )
2482
2483 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
2484 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
2485 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
2486
2487 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
2488 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
2489 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
2490 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
2491 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
2492 .le. eigen(2) .and. eigen(2) .le. eigen(1))
2493 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
2494 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
2495 .le. eigen(3) .and. eigen(3) .le. eigen(1))
2496
2497 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2498
2499 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2500 end do
2501 end do
2502 end subroutine sx_lambda2_lx7
2503
2504 subroutine sx_lambda2_lx6(lambda2, u, v, w, dx, dy, dz, &
2505 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
2506 integer, parameter :: lx = 6
2507 integer, intent(in) :: n
2508 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
2509 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
2510 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
2511 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
2512 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2513 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
2514 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
2515 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
2516 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
2517 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2518 real(kind=rp) :: grad(lx*lx*lx,3,3)
2519 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
2520 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
2521 real(kind=rp) :: a11, a22, a33, a12, a13, a23
2522 real(kind=rp) :: msk1, msk2, msk3
2523 real(kind=rp) :: ur(lx, lx, lx)
2524 real(kind=rp) :: vr(lx, lx, lx)
2525 real(kind=rp) :: wr(lx, lx, lx)
2526 real(kind=rp) :: us(lx, lx, lx)
2527 real(kind=rp) :: vs(lx, lx, lx)
2528 real(kind=rp) :: ws(lx, lx, lx)
2529 real(kind=rp) :: ut(lx, lx, lx)
2530 real(kind=rp) :: vt(lx, lx, lx)
2531 real(kind=rp) :: wt(lx, lx, lx)
2532 real(kind=rp) :: tmp1, tmp2, tmp3
2533 integer :: e, i, j, k, l
2534
2535 do e = 1, n
2536 do j = 1, lx * lx
2537 do i = 1, lx
2538 tmp1 = 0.0_rp
2539 tmp2 = 0.0_rp
2540 tmp3 = 0.0_rp
2541 do k = 1, lx
2542 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
2543 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
2544 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
2545 end do
2546 ur(i,j,1) = tmp1
2547 vr(i,j,1) = tmp2
2548 wr(i,j,1) = tmp3
2549 end do
2550 end do
2551
2552 do k = 1, lx
2553 do j = 1, lx
2554 do i = 1, lx
2555 tmp1 = 0.0_rp
2556 tmp2 = 0.0_rp
2557 tmp3 = 0.0_rp
2558 do l = 1, lx
2559 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
2560 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
2561 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
2562 end do
2563 us(i,j,k) = tmp1
2564 vs(i,j,k) = tmp2
2565 ws(i,j,k) = tmp3
2566 end do
2567 end do
2568 end do
2569
2570 do k = 1, lx
2571 do i = 1, lx*lx
2572 tmp1 = 0.0_rp
2573 tmp2 = 0.0_rp
2574 tmp3 = 0.0_rp
2575 do l = 1, lx
2576 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
2577 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
2578 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
2579 end do
2580 ut(i,1,k) = tmp1
2581 vt(i,1,k) = tmp2
2582 wt(i,1,k) = tmp3
2583 end do
2584 end do
2585
2586 do i = 1, lx * lx * lx
2587 grad(1,1,1) = w3(i,1,1) &
2588 * ( drdx(i,1,1,e) * ur(i,1,1) &
2589 + dsdx(i,1,1,e) * us(i,1,1) &
2590 + dtdx(i,1,1,e) * ut(i,1,1) )
2591 grad(1,1,2) = w3(i,1,1) &
2592 * ( dsdy(i,1,1,e) * us(i,1,1) &
2593 + drdy(i,1,1,e) * ur(i,1,1) &
2594 + dtdy(i,1,1,e) * ut(i,1,1) )
2595 grad(1,1,3) = w3(i,1,1) &
2596 * ( dtdz(i,1,1,e) * ut(i,1,1) &
2597 + drdz(i,1,1,e) * ur(i,1,1) &
2598 + dsdz(i,1,1,e) * us(i,1,1) )
2599
2600 grad(1,2,1) = w3(i,1,1) &
2601 * ( drdx(i,1,1,e) * vr(i,1,1) &
2602 + dsdx(i,1,1,e) * vs(i,1,1) &
2603 + dtdx(i,1,1,e) * vt(i,1,1) )
2604 grad(1,2,2) = w3(i,1,1) &
2605 * ( dsdy(i,1,1,e) * vs(i,1,1) &
2606 + drdy(i,1,1,e) * vr(i,1,1) &
2607 + dtdy(i,1,1,e) * vt(i,1,1) )
2608 grad(1,2,3) = w3(i,1,1) &
2609 * ( dtdz(i,1,1,e) * vt(i,1,1) &
2610 + drdz(i,1,1,e) * vr(i,1,1) &
2611 + dsdz(i,1,1,e) * vs(i,1,1) )
2612
2613 grad(1,3,1) = w3(i,1,1) &
2614 * ( drdx(i,1,1,e) * wr(i,1,1) &
2615 + dsdx(i,1,1,e) * ws(i,1,1) &
2616 + dtdx(i,1,1,e) * wt(i,1,1) )
2617 grad(1,3,2) = w3(i,1,1) &
2618 * ( dsdy(i,1,1,e) * ws(i,1,1) &
2619 + drdy(i,1,1,e) * wr(i,1,1) &
2620 + dtdy(i,1,1,e) * wt(i,1,1) )
2621 grad(1,3,3) = w3(i,1,1) &
2622 * ( dtdz(i,1,1,e) * wt(i,1,1) &
2623 + drdz(i,1,1,e) * wr(i,1,1) &
2624 + dsdz(i,1,1,e) * ws(i,1,1) )
2625 end do
2626
2627
2628 do i = 1, lx * lx * lx
2629 s11 = grad(i,1,1)
2630 s22 = grad(i,2,2)
2631 s33 = grad(i,3,3)
2632
2633
2634 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
2635 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
2636 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
2637
2638 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
2639 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
2640 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
2641
2642 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
2643 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
2644 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
2645
2646 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
2647 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
2648 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
2649
2650
2651 b = -(a11 + a22 + a33)
2652 c = -(a12*a12 + a13*a13 + a23*a23 &
2653 - a11 * a22 - a11 * a33 - a22 * a33)
2654 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
2655 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
2656
2657
2658 q = (3.0 * c - b*b) / 9.0
2659 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
2660 theta = acos( r / sqrt(-q*q*q) )
2661
2662 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
2663 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
2664 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
2665
2666 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
2667 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
2668 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
2669 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
2670 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
2671 .le. eigen(2) .and. eigen(2) .le. eigen(1))
2672 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
2673 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
2674 .le. eigen(3) .and. eigen(3) .le. eigen(1))
2675
2676 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2677
2678 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2679 end do
2680 end do
2681 end subroutine sx_lambda2_lx6
2682
2683 subroutine sx_lambda2_lx5(lambda2, u, v, w, dx, dy, dz, &
2684 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
2685 integer, parameter :: lx = 5
2686 integer, intent(in) :: n
2687 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
2688 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
2689 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
2690 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
2691 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2692 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
2693 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
2694 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
2695 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
2696 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2697 real(kind=rp) :: grad(lx*lx*lx,3,3)
2698 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
2699 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
2700 real(kind=rp) :: a11, a22, a33, a12, a13, a23
2701 real(kind=rp) :: msk1, msk2, msk3
2702 real(kind=rp) :: ur(lx, lx, lx)
2703 real(kind=rp) :: vr(lx, lx, lx)
2704 real(kind=rp) :: wr(lx, lx, lx)
2705 real(kind=rp) :: us(lx, lx, lx)
2706 real(kind=rp) :: vs(lx, lx, lx)
2707 real(kind=rp) :: ws(lx, lx, lx)
2708 real(kind=rp) :: ut(lx, lx, lx)
2709 real(kind=rp) :: vt(lx, lx, lx)
2710 real(kind=rp) :: wt(lx, lx, lx)
2711 real(kind=rp) :: tmp1, tmp2, tmp3
2712 integer :: e, i, j, k, l
2713
2714 do e = 1, n
2715 do j = 1, lx * lx
2716 do i = 1, lx
2717 tmp1 = 0.0_rp
2718 tmp2 = 0.0_rp
2719 tmp3 = 0.0_rp
2720 do k = 1, lx
2721 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
2722 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
2723 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
2724 end do
2725 ur(i,j,1) = tmp1
2726 vr(i,j,1) = tmp2
2727 wr(i,j,1) = tmp3
2728 end do
2729 end do
2730
2731 do k = 1, lx
2732 do j = 1, lx
2733 do i = 1, lx
2734 tmp1 = 0.0_rp
2735 tmp2 = 0.0_rp
2736 tmp3 = 0.0_rp
2737 do l = 1, lx
2738 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
2739 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
2740 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
2741 end do
2742 us(i,j,k) = tmp1
2743 vs(i,j,k) = tmp2
2744 ws(i,j,k) = tmp3
2745 end do
2746 end do
2747 end do
2748
2749 do k = 1, lx
2750 do i = 1, lx*lx
2751 tmp1 = 0.0_rp
2752 tmp2 = 0.0_rp
2753 tmp3 = 0.0_rp
2754 do l = 1, lx
2755 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
2756 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
2757 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
2758 end do
2759 ut(i,1,k) = tmp1
2760 vt(i,1,k) = tmp2
2761 wt(i,1,k) = tmp3
2762 end do
2763 end do
2764
2765 do i = 1, lx * lx * lx
2766 grad(1,1,1) = w3(i,1,1) &
2767 * ( drdx(i,1,1,e) * ur(i,1,1) &
2768 + dsdx(i,1,1,e) * us(i,1,1) &
2769 + dtdx(i,1,1,e) * ut(i,1,1) )
2770 grad(1,1,2) = w3(i,1,1) &
2771 * ( dsdy(i,1,1,e) * us(i,1,1) &
2772 + drdy(i,1,1,e) * ur(i,1,1) &
2773 + dtdy(i,1,1,e) * ut(i,1,1) )
2774 grad(1,1,3) = w3(i,1,1) &
2775 * ( dtdz(i,1,1,e) * ut(i,1,1) &
2776 + drdz(i,1,1,e) * ur(i,1,1) &
2777 + dsdz(i,1,1,e) * us(i,1,1) )
2778
2779 grad(1,2,1) = w3(i,1,1) &
2780 * ( drdx(i,1,1,e) * vr(i,1,1) &
2781 + dsdx(i,1,1,e) * vs(i,1,1) &
2782 + dtdx(i,1,1,e) * vt(i,1,1) )
2783 grad(1,2,2) = w3(i,1,1) &
2784 * ( dsdy(i,1,1,e) * vs(i,1,1) &
2785 + drdy(i,1,1,e) * vr(i,1,1) &
2786 + dtdy(i,1,1,e) * vt(i,1,1) )
2787 grad(1,2,3) = w3(i,1,1) &
2788 * ( dtdz(i,1,1,e) * vt(i,1,1) &
2789 + drdz(i,1,1,e) * vr(i,1,1) &
2790 + dsdz(i,1,1,e) * vs(i,1,1) )
2791
2792 grad(1,3,1) = w3(i,1,1) &
2793 * ( drdx(i,1,1,e) * wr(i,1,1) &
2794 + dsdx(i,1,1,e) * ws(i,1,1) &
2795 + dtdx(i,1,1,e) * wt(i,1,1) )
2796 grad(1,3,2) = w3(i,1,1) &
2797 * ( dsdy(i,1,1,e) * ws(i,1,1) &
2798 + drdy(i,1,1,e) * wr(i,1,1) &
2799 + dtdy(i,1,1,e) * wt(i,1,1) )
2800 grad(1,3,3) = w3(i,1,1) &
2801 * ( dtdz(i,1,1,e) * wt(i,1,1) &
2802 + drdz(i,1,1,e) * wr(i,1,1) &
2803 + dsdz(i,1,1,e) * ws(i,1,1) )
2804 end do
2805
2806
2807 do i = 1, lx * lx * lx
2808 s11 = grad(i,1,1)
2809 s22 = grad(i,2,2)
2810 s33 = grad(i,3,3)
2811
2812
2813 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
2814 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
2815 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
2816
2817 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
2818 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
2819 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
2820
2821 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
2822 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
2823 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
2824
2825 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
2826 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
2827 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
2828
2829
2830 b = -(a11 + a22 + a33)
2831 c = -(a12*a12 + a13*a13 + a23*a23 &
2832 - a11 * a22 - a11 * a33 - a22 * a33)
2833 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
2834 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
2835
2836
2837 q = (3.0 * c - b*b) / 9.0
2838 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
2839 theta = acos( r / sqrt(-q*q*q) )
2840
2841 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
2842 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
2843 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
2844
2845 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
2846 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
2847 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
2848 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
2849 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
2850 .le. eigen(2) .and. eigen(2) .le. eigen(1))
2851 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
2852 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
2853 .le. eigen(3) .and. eigen(3) .le. eigen(1))
2854
2855 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2856
2857 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2858 end do
2859 end do
2860 end subroutine sx_lambda2_lx5
2861
2862 subroutine sx_lambda2_lx4(lambda2, u, v, w, dx, dy, dz, &
2863 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
2864 integer, parameter :: lx = 4
2865 integer, intent(in) :: n
2866 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
2867 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
2868 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
2869 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
2870 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2871 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
2872 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
2873 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
2874 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
2875 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2876 real(kind=rp) :: grad(lx*lx*lx,3,3)
2877 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
2878 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
2879 real(kind=rp) :: a11, a22, a33, a12, a13, a23
2880 real(kind=rp) :: msk1, msk2, msk3
2881 real(kind=rp) :: ur(lx, lx, lx)
2882 real(kind=rp) :: vr(lx, lx, lx)
2883 real(kind=rp) :: wr(lx, lx, lx)
2884 real(kind=rp) :: us(lx, lx, lx)
2885 real(kind=rp) :: vs(lx, lx, lx)
2886 real(kind=rp) :: ws(lx, lx, lx)
2887 real(kind=rp) :: ut(lx, lx, lx)
2888 real(kind=rp) :: vt(lx, lx, lx)
2889 real(kind=rp) :: wt(lx, lx, lx)
2890 real(kind=rp) :: tmp1, tmp2, tmp3
2891 integer :: e, i, j, k, l
2892
2893 do e = 1, n
2894 do j = 1, lx * lx
2895 do i = 1, lx
2896 tmp1 = 0.0_rp
2897 tmp2 = 0.0_rp
2898 tmp3 = 0.0_rp
2899 do k = 1, lx
2900 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
2901 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
2902 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
2903 end do
2904 ur(i,j,1) = tmp1
2905 vr(i,j,1) = tmp2
2906 wr(i,j,1) = tmp3
2907 end do
2908 end do
2909
2910 do k = 1, lx
2911 do j = 1, lx
2912 do i = 1, lx
2913 tmp1 = 0.0_rp
2914 tmp2 = 0.0_rp
2915 tmp3 = 0.0_rp
2916 do l = 1, lx
2917 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
2918 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
2919 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
2920 end do
2921 us(i,j,k) = tmp1
2922 vs(i,j,k) = tmp2
2923 ws(i,j,k) = tmp3
2924 end do
2925 end do
2926 end do
2927
2928 do k = 1, lx
2929 do i = 1, lx*lx
2930 tmp1 = 0.0_rp
2931 tmp2 = 0.0_rp
2932 tmp3 = 0.0_rp
2933 do l = 1, lx
2934 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
2935 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
2936 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
2937 end do
2938 ut(i,1,k) = tmp1
2939 vt(i,1,k) = tmp2
2940 wt(i,1,k) = tmp3
2941 end do
2942 end do
2943
2944 do i = 1, lx * lx * lx
2945 grad(1,1,1) = w3(i,1,1) &
2946 * ( drdx(i,1,1,e) * ur(i,1,1) &
2947 + dsdx(i,1,1,e) * us(i,1,1) &
2948 + dtdx(i,1,1,e) * ut(i,1,1) )
2949 grad(1,1,2) = w3(i,1,1) &
2950 * ( dsdy(i,1,1,e) * us(i,1,1) &
2951 + drdy(i,1,1,e) * ur(i,1,1) &
2952 + dtdy(i,1,1,e) * ut(i,1,1) )
2953 grad(1,1,3) = w3(i,1,1) &
2954 * ( dtdz(i,1,1,e) * ut(i,1,1) &
2955 + drdz(i,1,1,e) * ur(i,1,1) &
2956 + dsdz(i,1,1,e) * us(i,1,1) )
2957
2958 grad(1,2,1) = w3(i,1,1) &
2959 * ( drdx(i,1,1,e) * vr(i,1,1) &
2960 + dsdx(i,1,1,e) * vs(i,1,1) &
2961 + dtdx(i,1,1,e) * vt(i,1,1) )
2962 grad(1,2,2) = w3(i,1,1) &
2963 * ( dsdy(i,1,1,e) * vs(i,1,1) &
2964 + drdy(i,1,1,e) * vr(i,1,1) &
2965 + dtdy(i,1,1,e) * vt(i,1,1) )
2966 grad(1,2,3) = w3(i,1,1) &
2967 * ( dtdz(i,1,1,e) * vt(i,1,1) &
2968 + drdz(i,1,1,e) * vr(i,1,1) &
2969 + dsdz(i,1,1,e) * vs(i,1,1) )
2970
2971 grad(1,3,1) = w3(i,1,1) &
2972 * ( drdx(i,1,1,e) * wr(i,1,1) &
2973 + dsdx(i,1,1,e) * ws(i,1,1) &
2974 + dtdx(i,1,1,e) * wt(i,1,1) )
2975 grad(1,3,2) = w3(i,1,1) &
2976 * ( dsdy(i,1,1,e) * ws(i,1,1) &
2977 + drdy(i,1,1,e) * wr(i,1,1) &
2978 + dtdy(i,1,1,e) * wt(i,1,1) )
2979 grad(1,3,3) = w3(i,1,1) &
2980 * ( dtdz(i,1,1,e) * wt(i,1,1) &
2981 + drdz(i,1,1,e) * wr(i,1,1) &
2982 + dsdz(i,1,1,e) * ws(i,1,1) )
2983 end do
2984
2985
2986 do i = 1, lx * lx * lx
2987 s11 = grad(i,1,1)
2988 s22 = grad(i,2,2)
2989 s33 = grad(i,3,3)
2990
2991
2992 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
2993 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
2994 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
2995
2996 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
2997 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
2998 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
2999
3000 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
3001 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
3002 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
3003
3004 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
3005 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
3006 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
3007
3008
3009 b = -(a11 + a22 + a33)
3010 c = -(a12*a12 + a13*a13 + a23*a23 &
3011 - a11 * a22 - a11 * a33 - a22 * a33)
3012 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
3013 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
3014
3015
3016 q = (3.0 * c - b*b) / 9.0
3017 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
3018 theta = acos( r / sqrt(-q*q*q) )
3019
3020 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
3021 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
3022 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
3023
3024 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
3025 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
3026 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
3027 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
3028 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
3029 .le. eigen(2) .and. eigen(2) .le. eigen(1))
3030 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
3031 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
3032 .le. eigen(3) .and. eigen(3) .le. eigen(1))
3033
3034 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
3035
3036 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
3037 end do
3038 end do
3039 end subroutine sx_lambda2_lx4
3040
3041 subroutine sx_lambda2_lx3(lambda2, u, v, w, dx, dy, dz, &
3042 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
3043 integer, parameter :: lx = 3
3044 integer, intent(in) :: n
3045 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
3046 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
3047 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
3048 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
3049 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
3050 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
3051 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
3052 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
3053 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
3054 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
3055 real(kind=rp) :: grad(lx*lx*lx,3,3)
3056 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
3057 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
3058 real(kind=rp) :: a11, a22, a33, a12, a13, a23
3059 real(kind=rp) :: msk1, msk2, msk3
3060 real(kind=rp) :: ur(lx, lx, lx)
3061 real(kind=rp) :: vr(lx, lx, lx)
3062 real(kind=rp) :: wr(lx, lx, lx)
3063 real(kind=rp) :: us(lx, lx, lx)
3064 real(kind=rp) :: vs(lx, lx, lx)
3065 real(kind=rp) :: ws(lx, lx, lx)
3066 real(kind=rp) :: ut(lx, lx, lx)
3067 real(kind=rp) :: vt(lx, lx, lx)
3068 real(kind=rp) :: wt(lx, lx, lx)
3069 real(kind=rp) :: tmp1, tmp2, tmp3
3070 integer :: e, i, j, k, l
3071
3072 do e = 1, n
3073 do j = 1, lx * lx
3074 do i = 1, lx
3075 tmp1 = 0.0_rp
3076 tmp2 = 0.0_rp
3077 tmp3 = 0.0_rp
3078 do k = 1, lx
3079 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
3080 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
3081 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
3082 end do
3083 ur(i,j,1) = tmp1
3084 vr(i,j,1) = tmp2
3085 wr(i,j,1) = tmp3
3086 end do
3087 end do
3088
3089 do k = 1, lx
3090 do j = 1, lx
3091 do i = 1, lx
3092 tmp1 = 0.0_rp
3093 tmp2 = 0.0_rp
3094 tmp3 = 0.0_rp
3095 do l = 1, lx
3096 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
3097 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
3098 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
3099 end do
3100 us(i,j,k) = tmp1
3101 vs(i,j,k) = tmp2
3102 ws(i,j,k) = tmp3
3103 end do
3104 end do
3105 end do
3106
3107 do k = 1, lx
3108 do i = 1, lx*lx
3109 tmp1 = 0.0_rp
3110 tmp2 = 0.0_rp
3111 tmp3 = 0.0_rp
3112 do l = 1, lx
3113 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
3114 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
3115 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
3116 end do
3117 ut(i,1,k) = tmp1
3118 vt(i,1,k) = tmp2
3119 wt(i,1,k) = tmp3
3120 end do
3121 end do
3122
3123 do i = 1, lx * lx * lx
3124 grad(1,1,1) = w3(i,1,1) &
3125 * ( drdx(i,1,1,e) * ur(i,1,1) &
3126 + dsdx(i,1,1,e) * us(i,1,1) &
3127 + dtdx(i,1,1,e) * ut(i,1,1) )
3128 grad(1,1,2) = w3(i,1,1) &
3129 * ( dsdy(i,1,1,e) * us(i,1,1) &
3130 + drdy(i,1,1,e) * ur(i,1,1) &
3131 + dtdy(i,1,1,e) * ut(i,1,1) )
3132 grad(1,1,3) = w3(i,1,1) &
3133 * ( dtdz(i,1,1,e) * ut(i,1,1) &
3134 + drdz(i,1,1,e) * ur(i,1,1) &
3135 + dsdz(i,1,1,e) * us(i,1,1) )
3136
3137 grad(1,2,1) = w3(i,1,1) &
3138 * ( drdx(i,1,1,e) * vr(i,1,1) &
3139 + dsdx(i,1,1,e) * vs(i,1,1) &
3140 + dtdx(i,1,1,e) * vt(i,1,1) )
3141 grad(1,2,2) = w3(i,1,1) &
3142 * ( dsdy(i,1,1,e) * vs(i,1,1) &
3143 + drdy(i,1,1,e) * vr(i,1,1) &
3144 + dtdy(i,1,1,e) * vt(i,1,1) )
3145 grad(1,2,3) = w3(i,1,1) &
3146 * ( dtdz(i,1,1,e) * vt(i,1,1) &
3147 + drdz(i,1,1,e) * vr(i,1,1) &
3148 + dsdz(i,1,1,e) * vs(i,1,1) )
3149
3150 grad(1,3,1) = w3(i,1,1) &
3151 * ( drdx(i,1,1,e) * wr(i,1,1) &
3152 + dsdx(i,1,1,e) * ws(i,1,1) &
3153 + dtdx(i,1,1,e) * wt(i,1,1) )
3154 grad(1,3,2) = w3(i,1,1) &
3155 * ( dsdy(i,1,1,e) * ws(i,1,1) &
3156 + drdy(i,1,1,e) * wr(i,1,1) &
3157 + dtdy(i,1,1,e) * wt(i,1,1) )
3158 grad(1,3,3) = w3(i,1,1) &
3159 * ( dtdz(i,1,1,e) * wt(i,1,1) &
3160 + drdz(i,1,1,e) * wr(i,1,1) &
3161 + dsdz(i,1,1,e) * ws(i,1,1) )
3162 end do
3163
3164
3165 do i = 1, lx * lx * lx
3166 s11 = grad(i,1,1)
3167 s22 = grad(i,2,2)
3168 s33 = grad(i,3,3)
3169
3170
3171 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
3172 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
3173 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
3174
3175 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
3176 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
3177 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
3178
3179 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
3180 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
3181 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
3182
3183 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
3184 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
3185 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
3186
3187
3188 b = -(a11 + a22 + a33)
3189 c = -(a12*a12 + a13*a13 + a23*a23 &
3190 - a11 * a22 - a11 * a33 - a22 * a33)
3191 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
3192 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
3193
3194
3195 q = (3.0 * c - b*b) / 9.0
3196 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
3197 theta = acos( r / sqrt(-q*q*q) )
3198
3199 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
3200 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
3201 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
3202
3203 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
3204 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
3205 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
3206 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
3207 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
3208 .le. eigen(2) .and. eigen(2) .le. eigen(1))
3209 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
3210 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
3211 .le. eigen(3) .and. eigen(3) .le. eigen(1))
3212
3213 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
3214
3215 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
3216 end do
3217 end do
3218 end subroutine sx_lambda2_lx3
3219
3220 subroutine sx_lambda2_lx2(lambda2, u, v, w, dx, dy, dz, &
3221 drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, cB, n)
3222 integer, parameter :: lx = 2
3223 integer, intent(in) :: n
3224 real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: lambda2
3225 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
3226 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: v
3227 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: w
3228 real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
3229 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
3230 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
3231 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
3232 real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: cb
3233 real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
3234 real(kind=rp) :: grad(lx*lx*lx,3,3)
3235 real(kind=rp) :: eigen(3), b, c, d, q, r, theta, l2
3236 real(kind=rp) :: s11, s22, s33, s12, s13, s23, o12, o13, o23
3237 real(kind=rp) :: a11, a22, a33, a12, a13, a23
3238 real(kind=rp) :: msk1, msk2, msk3
3239 real(kind=rp) :: ur(lx, lx, lx)
3240 real(kind=rp) :: vr(lx, lx, lx)
3241 real(kind=rp) :: wr(lx, lx, lx)
3242 real(kind=rp) :: us(lx, lx, lx)
3243 real(kind=rp) :: vs(lx, lx, lx)
3244 real(kind=rp) :: ws(lx, lx, lx)
3245 real(kind=rp) :: ut(lx, lx, lx)
3246 real(kind=rp) :: vt(lx, lx, lx)
3247 real(kind=rp) :: wt(lx, lx, lx)
3248 real(kind=rp) :: tmp1, tmp2, tmp3
3249 integer :: e, i, j, k, l
3250
3251 do e = 1, n
3252 do j = 1, lx * lx
3253 do i = 1, lx
3254 tmp1 = 0.0_rp
3255 tmp2 = 0.0_rp
3256 tmp3 = 0.0_rp
3257 do k = 1, lx
3258 tmp1 = tmp1 + dx(i,k) * u(k,j,1,e)
3259 tmp2 = tmp2 + dx(i,k) * v(k,j,1,e)
3260 tmp3 = tmp3 + dx(i,k) * w(k,j,1,e)
3261 end do
3262 ur(i,j,1) = tmp1
3263 vr(i,j,1) = tmp2
3264 wr(i,j,1) = tmp3
3265 end do
3266 end do
3267
3268 do k = 1, lx
3269 do j = 1, lx
3270 do i = 1, lx
3271 tmp1 = 0.0_rp
3272 tmp2 = 0.0_rp
3273 tmp3 = 0.0_rp
3274 do l = 1, lx
3275 tmp1 = tmp1 + dy(j,l) * u(i,l,k,e)
3276 tmp2 = tmp2 + dy(j,l) * v(i,l,k,e)
3277 tmp3 = tmp3 + dy(j,l) * w(i,l,k,e)
3278 end do
3279 us(i,j,k) = tmp1
3280 vs(i,j,k) = tmp2
3281 ws(i,j,k) = tmp3
3282 end do
3283 end do
3284 end do
3285
3286 do k = 1, lx
3287 do i = 1, lx*lx
3288 tmp1 = 0.0_rp
3289 tmp2 = 0.0_rp
3290 tmp3 = 0.0_rp
3291 do l = 1, lx
3292 tmp1 = tmp1 + dz(k,l) * u(i,1,l,e)
3293 tmp2 = tmp2 + dz(k,l) * v(i,1,l,e)
3294 tmp3 = tmp3 + dz(k,l) * w(i,1,l,e)
3295 end do
3296 ut(i,1,k) = tmp1
3297 vt(i,1,k) = tmp2
3298 wt(i,1,k) = tmp3
3299 end do
3300 end do
3301
3302 do i = 1, lx * lx * lx
3303 grad(1,1,1) = w3(i,1,1) &
3304 * ( drdx(i,1,1,e) * ur(i,1,1) &
3305 + dsdx(i,1,1,e) * us(i,1,1) &
3306 + dtdx(i,1,1,e) * ut(i,1,1) )
3307 grad(1,1,2) = w3(i,1,1) &
3308 * ( dsdy(i,1,1,e) * us(i,1,1) &
3309 + drdy(i,1,1,e) * ur(i,1,1) &
3310 + dtdy(i,1,1,e) * ut(i,1,1) )
3311 grad(1,1,3) = w3(i,1,1) &
3312 * ( dtdz(i,1,1,e) * ut(i,1,1) &
3313 + drdz(i,1,1,e) * ur(i,1,1) &
3314 + dsdz(i,1,1,e) * us(i,1,1) )
3315
3316 grad(1,2,1) = w3(i,1,1) &
3317 * ( drdx(i,1,1,e) * vr(i,1,1) &
3318 + dsdx(i,1,1,e) * vs(i,1,1) &
3319 + dtdx(i,1,1,e) * vt(i,1,1) )
3320 grad(1,2,2) = w3(i,1,1) &
3321 * ( dsdy(i,1,1,e) * vs(i,1,1) &
3322 + drdy(i,1,1,e) * vr(i,1,1) &
3323 + dtdy(i,1,1,e) * vt(i,1,1) )
3324 grad(1,2,3) = w3(i,1,1) &
3325 * ( dtdz(i,1,1,e) * vt(i,1,1) &
3326 + drdz(i,1,1,e) * vr(i,1,1) &
3327 + dsdz(i,1,1,e) * vs(i,1,1) )
3328
3329 grad(1,3,1) = w3(i,1,1) &
3330 * ( drdx(i,1,1,e) * wr(i,1,1) &
3331 + dsdx(i,1,1,e) * ws(i,1,1) &
3332 + dtdx(i,1,1,e) * wt(i,1,1) )
3333 grad(1,3,2) = w3(i,1,1) &
3334 * ( dsdy(i,1,1,e) * ws(i,1,1) &
3335 + drdy(i,1,1,e) * wr(i,1,1) &
3336 + dtdy(i,1,1,e) * wt(i,1,1) )
3337 grad(1,3,3) = w3(i,1,1) &
3338 * ( dtdz(i,1,1,e) * wt(i,1,1) &
3339 + drdz(i,1,1,e) * wr(i,1,1) &
3340 + dsdz(i,1,1,e) * ws(i,1,1) )
3341 end do
3342
3343
3344 do i = 1, lx * lx * lx
3345 s11 = grad(i,1,1)
3346 s22 = grad(i,2,2)
3347 s33 = grad(i,3,3)
3348
3349
3350 s12 = 0.5*(grad(i,1,2) + grad(i,2,1))
3351 s13 = 0.5*(grad(i,1,3) + grad(i,3,1))
3352 s23 = 0.5*(grad(i,2,3) + grad(i,3,2))
3353
3354 o12 = 0.5*(grad(i,1,2) - grad(i,2,1))
3355 o13 = 0.5*(grad(i,1,3) - grad(i,3,1))
3356 o23 = 0.5*(grad(i,2,3) - grad(i,3,2))
3357
3358 a11 = s11*s11 + s12*s12 + s13*s13 - o12*o12 - o13*o13
3359 a12 = s11 * s12 + s12 * s22 + s13 * s23 - o13 * o23
3360 a13 = s11 * s13 + s12 * s23 + s13 * s33 + o12 * o23
3361
3362 a22 = s12*s12 + s22*s22 + s23*s23 - o12*o12 - o23*o23
3363 a23 = s12 * s13 + s22 * s23 + s23 * s33 - o12 * o13
3364 a33 = s13*s13 + s23*s23 + s33*s33 - o13*o13 - o23*o23
3365
3366
3367 b = -(a11 + a22 + a33)
3368 c = -(a12*a12 + a13*a13 + a23*a23 &
3369 - a11 * a22 - a11 * a33 - a22 * a33)
3370 d = -(2.0 * a12 * a13 * a23 - a11 * a23*a23 &
3371 - a22 * a13*a13 - a33 * a12*a12 + a11 * a22 * a33)
3372
3373
3374 q = (3.0 * c - b*b) / 9.0
3375 r = (9.0 * c * b - 27.0 * d - 2.0 * b*b*b) / 54.0
3376 theta = acos( r / sqrt(-q*q*q) )
3377
3378 eigen(1) = 2.0 * sqrt(-q) * cos(theta / 3.0) - b / 3.0
3379 eigen(2) = 2.0 * sqrt(-q) * cos((theta + 2.0 * pi) / 3.0) - b / 3.0
3380 eigen(3) = 2.0 * sqrt(-q) * cos((theta + 4.0 * pi) / 3.0) - b / 3.0
3381
3382 msk1 = merge(1.0_rp, 0.0_rp, eigen(2) .le. eigen(1) &
3383 .and. eigen(1) .le. eigen(3) .or. eigen(3) &
3384 .le. eigen(1) .and. eigen(1) .le. eigen(2) )
3385 msk2 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(2) &
3386 .and. eigen(2) .le. eigen(3) .or. eigen(3) &
3387 .le. eigen(2) .and. eigen(2) .le. eigen(1))
3388 msk3 = merge(1.0_rp, 0.0_rp, eigen(1) .le. eigen(3) &
3389 .and. eigen(3) .le. eigen(2) .or. eigen(2) &
3390 .le. eigen(3) .and. eigen(3) .le. eigen(1))
3391
3392 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
3393
3394 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
3395 end do
3396 end do
3397 end subroutine sx_lambda2_lx2
3398
3399end submodule sx_lambda2
3400
3401
lambda2
A simulation component that computes lambda2 The values are stored in the field registry under the na...
Definition lambda2.f90:37
math
Definition math.f90:60
math::pi
real(kind=rp), parameter, public pi
Definition math.f90:75
opr_sx
Operators SX-Aurora backend.
Definition opr_sx.f90:2