34submodule(
opr_sx) sx_lambda2
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
45 associate(xh => coef%Xh, msh => coef%msh)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
176 end subroutine opr_sx_lambda2
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
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)
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)
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)
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) )
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) )
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) )
301 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
349 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
351 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
354 end subroutine sx_lambda2_lx
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
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)
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)
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)
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) )
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) )
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) )
480 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
528 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
530 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
533 end subroutine sx_lambda2_lx18
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
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)
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)
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)
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) )
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) )
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) )
659 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
707 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
709 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
712 end subroutine sx_lambda2_lx17
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
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)
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)
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)
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) )
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) )
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) )
838 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
886 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
888 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
891 end subroutine sx_lambda2_lx16
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
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)
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)
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)
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) )
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) )
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) )
1017 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
1065 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1067 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1070 end subroutine sx_lambda2_lx15
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
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)
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)
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)
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) )
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) )
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) )
1196 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
1244 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1246 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1249 end subroutine sx_lambda2_lx14
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
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)
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)
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)
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) )
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) )
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) )
1375 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
1423 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1425 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1428 end subroutine sx_lambda2_lx13
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
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)
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)
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)
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) )
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) )
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) )
1554 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
1602 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1604 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1607 end subroutine sx_lambda2_lx12
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
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)
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)
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)
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) )
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) )
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) )
1733 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
1781 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1783 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1786 end subroutine sx_lambda2_lx11
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
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)
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)
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)
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) )
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) )
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) )
1912 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
1960 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
1962 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
1965 end subroutine sx_lambda2_lx10
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
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)
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)
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)
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) )
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) )
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) )
2091 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
2139 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2141 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2144 end subroutine sx_lambda2_lx9
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
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)
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)
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)
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) )
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) )
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) )
2270 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
2318 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2320 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2323 end subroutine sx_lambda2_lx8
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
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)
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)
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)
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) )
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) )
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) )
2449 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
2497 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2499 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2502 end subroutine sx_lambda2_lx7
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
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)
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)
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)
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) )
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) )
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) )
2628 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
2676 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2678 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2681 end subroutine sx_lambda2_lx6
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
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)
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)
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)
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) )
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) )
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) )
2807 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
2855 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
2857 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
2860 end subroutine sx_lambda2_lx5
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
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)
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)
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)
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) )
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) )
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) )
2986 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
3034 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
3036 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
3039 end subroutine sx_lambda2_lx4
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
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)
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)
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)
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) )
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) )
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) )
3165 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
3213 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
3215 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
3218 end subroutine sx_lambda2_lx3
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
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)
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)
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)
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) )
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) )
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) )
3344 do i = 1, lx * lx * lx
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))
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))
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
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
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)
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) )
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
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))
3392 l2 = msk1 * eigen(1) + msk2 * eigen(2) + msk3 * eigen(3)
3394 lambda2(i,1,1,e) = l2/(cb(i,1,1,e)**2)
3397 end subroutine sx_lambda2_lx2
3399end submodule sx_lambda2
A simulation component that computes lambda2 The values are stored in the field registry under the na...
real(kind=rp), parameter, public pi
Operators SX-Aurora backend.