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