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