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