Neko 1.99.3
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
ax_helm_full_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 2023-2026, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
33module ax_helm_full_cpu
34 use ax_helm_full, only : ax_helm_full_t
35 use num_types, only : rp
36 use coefs, only : coef_t
37 use space, only : space_t
38 use mesh, only : mesh_t
39 implicit none
40 private
41
43 type, public, extends(ax_helm_full_t) :: ax_helm_full_cpu_t
44 contains
46 procedure, pass(this) :: compute_vector => ax_helm_full_compute_vector
47 end type ax_helm_full_cpu_t
48
49contains
50
62 subroutine ax_helm_full_compute_vector(this, au, av, aw, u, v, w, coef, msh, &
63 Xh)
64 class(ax_helm_full_cpu_t), intent(in) :: this
65 type(mesh_t), intent(in) :: msh
66 type(space_t), intent(in) :: Xh
67 type(coef_t), intent(in) :: coef
68 real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, msh%nelv)
69 real(kind=rp), intent(in) :: v(xh%lx, xh%ly, xh%lz, msh%nelv)
70 real(kind=rp), intent(in) :: w(xh%lx, xh%ly, xh%lz, msh%nelv)
71 real(kind=rp), intent(inout) :: au(xh%lx, xh%ly, xh%lz, msh%nelv)
72 real(kind=rp), intent(inout) :: av(xh%lx, xh%ly, xh%lz, msh%nelv)
73 real(kind=rp), intent(inout) :: aw(xh%lx, xh%ly, xh%lz, msh%nelv)
74 integer :: i
75
76 !$omp parallel
77 select case (xh%lx)
78 case (14)
79 call ax_helm_stress_lx14(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
80 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
81 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
82 coef%dtdx, coef%dtdy, coef%dtdz, &
83 coef%jacinv, xh%w3, msh%nelv)
84 case (13)
85 call ax_helm_stress_lx13(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
86 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
87 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
88 coef%dtdx, coef%dtdy, coef%dtdz, &
89 coef%jacinv, xh%w3, msh%nelv)
90 case (12)
91 call ax_helm_stress_lx12(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
92 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
93 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
94 coef%dtdx, coef%dtdy, coef%dtdz, &
95 coef%jacinv, xh%w3, msh%nelv)
96 case (11)
97 call ax_helm_stress_lx11(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
98 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
99 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
100 coef%dtdx, coef%dtdy, coef%dtdz, &
101 coef%jacinv, xh%w3, msh%nelv)
102 case (10)
103 call ax_helm_stress_lx10(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
104 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
105 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
106 coef%dtdx, coef%dtdy, coef%dtdz, &
107 coef%jacinv, xh%w3, msh%nelv)
108 case (9)
109 call ax_helm_stress_lx9(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
110 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
111 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
112 coef%dtdx, coef%dtdy, coef%dtdz, &
113 coef%jacinv, xh%w3, msh%nelv)
114 case (8)
115 call ax_helm_stress_lx8(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
116 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
117 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
118 coef%dtdx, coef%dtdy, coef%dtdz, &
119 coef%jacinv, xh%w3, msh%nelv)
120 case (7)
121 call ax_helm_stress_lx7(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
122 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
123 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
124 coef%dtdx, coef%dtdy, coef%dtdz, &
125 coef%jacinv, xh%w3, msh%nelv)
126 case (6)
127 call ax_helm_stress_lx6(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
128 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
129 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
130 coef%dtdx, coef%dtdy, coef%dtdz, &
131 coef%jacinv, xh%w3, msh%nelv)
132 case (5)
133 call ax_helm_stress_lx5(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
134 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
135 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
136 coef%dtdx, coef%dtdy, coef%dtdz, &
137 coef%jacinv, xh%w3, msh%nelv)
138 case (4)
139 call ax_helm_stress_lx4(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
140 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
141 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
142 coef%dtdx, coef%dtdy, coef%dtdz, &
143 coef%jacinv, xh%w3, msh%nelv)
144 case (3)
145 call ax_helm_stress_lx3(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
146 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
147 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
148 coef%dtdx, coef%dtdy, coef%dtdz, &
149 coef%jacinv, xh%w3, msh%nelv)
150 case (2)
151 call ax_helm_stress_lx2(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
152 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
153 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
154 coef%dtdx, coef%dtdy, coef%dtdz, &
155 coef%jacinv, xh%w3, msh%nelv)
156 case default
157 call ax_helm_stress_lx(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
158 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%h2, &
159 coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, coef%dsdz, &
160 coef%dtdx, coef%dtdy, coef%dtdz, &
161 coef%jacinv, xh%w3, msh%nelv, xh%lx)
162 end select
163
164 if (coef%ifh2) then
165 !$omp do private(i)
166 do i = 1, coef%dof%size()
167 au(i,1,1,1) = au(i,1,1,1) + &
168 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * u(i,1,1,1)
169 av(i,1,1,1) = av(i,1,1,1) + &
170 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * v(i,1,1,1)
171 aw(i,1,1,1) = aw(i,1,1,1) + &
172 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * w(i,1,1,1)
173 end do
174 !$omp end do
175
176 end if
177 !$omp end parallel
178
179 end subroutine ax_helm_full_compute_vector
180
181 subroutine ax_helm_stress_lx(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
182 h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
183 jacinv, weights3, n, lx)
184
185 integer, intent(in) :: n, lx
186 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
187 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
188 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
189 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
190 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
191 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
192 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
193 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
194 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
195 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
196 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
197 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
198 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
199 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
200 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
201 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
202 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
203 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
204 real(kind=rp), intent(in) :: dx(lx, lx)
205 real(kind=rp), intent(in) :: dy(lx, lx)
206 real(kind=rp), intent(in) :: dz(lx, lx)
207 real(kind=rp), intent(in) :: dxt(lx, lx)
208 real(kind=rp), intent(in) :: dyt(lx, lx)
209 real(kind=rp), intent(in) :: dzt(lx, lx)
210 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
211
212 real(kind=rp) :: wur(lx, lx, lx)
213 real(kind=rp) :: wus(lx, lx, lx)
214 real(kind=rp) :: wut(lx, lx, lx)
215 real(kind=rp) :: wvr(lx, lx, lx)
216 real(kind=rp) :: wvs(lx, lx, lx)
217 real(kind=rp) :: wvt(lx, lx, lx)
218 real(kind=rp) :: wwr(lx, lx, lx)
219 real(kind=rp) :: wws(lx, lx, lx)
220 real(kind=rp) :: wwt(lx, lx, lx)
221
222 integer :: e, i, j, k, l
223 real(kind=rp) :: dj, t1, t2, t3
224 real(kind=rp) :: s11, s22, s33, s12, s13, s23
225 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
226
227 !$omp do private(e,i,j,k,l)
228 do e = 1, n
229
230 do j = 1, lx * lx
231 do i = 1, lx
232 t1 = 0.0_rp
233 t2 = 0.0_rp
234 t3 = 0.0_rp
235 do k = 1, lx
236 t1 = t1 + dx(i,k) * u(k,j,1,e)
237 t2 = t2 + dx(i,k) * v(k,j,1,e)
238 t3 = t3 + dx(i,k) * w(k,j,1,e)
239 end do
240 wur(i,j,1) = t1
241 wvr(i,j,1) = t2
242 wwr(i,j,1) = t3
243 end do
244 end do
245
246 do k = 1, lx
247 do j = 1, lx
248 do i = 1, lx
249 t1 = 0.0_rp
250 t2 = 0.0_rp
251 t3 = 0.0_rp
252 do l = 1, lx
253 t1 = t1 + dy(j,l) * u(i,l,k,e)
254 t2 = t2 + dy(j,l) * v(i,l,k,e)
255 t3 = t3 + dy(j,l) * w(i,l,k,e)
256 end do
257 wus(i,j,k) = t1
258 wvs(i,j,k) = t2
259 wws(i,j,k) = t3
260 end do
261 end do
262 end do
263
264 do k = 1, lx
265 do i = 1, lx*lx
266 t1 = 0.0_rp
267 t2 = 0.0_rp
268 t3 = 0.0_rp
269 do l = 1, lx
270 t1 = t1 + dz(k,l) * u(i,1,l,e)
271 t2 = t2 + dz(k,l) * v(i,1,l,e)
272 t3 = t3 + dz(k,l) * w(i,1,l,e)
273 end do
274 wut(i,1,k) = t1
275 wvt(i,1,k) = t2
276 wwt(i,1,k) = t3
277 end do
278 end do
279
280 do i = 1, lx*lx*lx
281
282 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
283 + wut(i,1,1) * dtdx(i,1,1,e)
284 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
285 + wut(i,1,1) * dtdy(i,1,1,e)
286 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
287 + wut(i,1,1) * dtdz(i,1,1,e)
288
289 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
290 + wvt(i,1,1) * dtdx(i,1,1,e)
291 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
292 + wvt(i,1,1) * dtdy(i,1,1,e)
293 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
294 + wvt(i,1,1) * dtdz(i,1,1,e)
295
296 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
297 + wwt(i,1,1) * dtdx(i,1,1,e)
298 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
299 + wwt(i,1,1) * dtdy(i,1,1,e)
300 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
301 + wwt(i,1,1) * dtdz(i,1,1,e)
302
303 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
304 s11 = dj*(u1 + u1)
305 s22 = dj*(v2 + v2)
306 s33 = dj*(w3 + w3)
307 s12 = dj*(u2 + v1)
308 s13 = dj*(u3 + w1)
309 s23 = dj*(v3 + w2)
310
311 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
312 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
313 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
314
315 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
316 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
317 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
318
319 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
320 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
321 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
322 end do
323
324 do j = 1, lx*lx
325 do i = 1, lx
326 t1 = 0.0_rp
327 t2 = 0.0_rp
328 t3 = 0.0_rp
329 do k = 1, lx
330 t1 = t1 + dxt(i,k) * wur(k,j,1)
331 t2 = t2 + dxt(i,k) * wvr(k,j,1)
332 t3 = t3 + dxt(i,k) * wwr(k,j,1)
333 end do
334 au(i,j,1,e) = t1
335 av(i,j,1,e) = t2
336 aw(i,j,1,e) = t3
337 end do
338 end do
339
340 do k = 1, lx
341 do j = 1, lx
342 do i = 1, lx
343 t1 = 0.0_rp
344 t2 = 0.0_rp
345 t3 = 0.0_rp
346 do l = 1, lx
347 t1 = t1 + dyt(j,l) * wus(i,l,k)
348 t2 = t2 + dyt(j,l) * wvs(i,l,k)
349 t3 = t3 + dyt(j,l) * wws(i,l,k)
350 end do
351 au(i,j,k,e) = au(i,j,k,e) + t1
352 av(i,j,k,e) = av(i,j,k,e) + t2
353 aw(i,j,k,e) = aw(i,j,k,e) + t3
354 end do
355 end do
356 end do
357
358 do k = 1, lx
359 do i = 1, lx*lx
360 t1 = 0.0_rp
361 t2 = 0.0_rp
362 t3 = 0.0_rp
363 do l = 1, lx
364 t1 = t1 + dzt(k,l) * wut(i,1,l)
365 t2 = t2 + dzt(k,l) * wvt(i,1,l)
366 t3 = t3 + dzt(k,l) * wwt(i,1,l)
367 end do
368 au(i,1,k,e) = au(i,1,k,e) + t1
369 av(i,1,k,e) = av(i,1,k,e) + t2
370 aw(i,1,k,e) = aw(i,1,k,e) + t3
371 end do
372 end do
373 end do
374 !$omp end do
375 end subroutine ax_helm_stress_lx
376
377 subroutine ax_helm_stress_lx14(au, av, aw, u, v, w, Dx, Dy, Dz, &
378 Dxt, Dyt, Dzt, &
379 h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
380 jacinv, weights3, n)
381 integer, parameter :: lx = 14
382 integer, intent(in) :: n
383 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
384 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
385 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
386 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
387 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
388 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
389 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
390 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
391 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
392 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
393 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
394 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
395 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
396 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
397 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
398 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
399 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
400 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
401 real(kind=rp), intent(in) :: dx(lx, lx)
402 real(kind=rp), intent(in) :: dy(lx, lx)
403 real(kind=rp), intent(in) :: dz(lx, lx)
404 real(kind=rp), intent(in) :: dxt(lx, lx)
405 real(kind=rp), intent(in) :: dyt(lx, lx)
406 real(kind=rp), intent(in) :: dzt(lx, lx)
407 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
408
409 real(kind=rp) :: wur(lx, lx, lx)
410 real(kind=rp) :: wus(lx, lx, lx)
411 real(kind=rp) :: wut(lx, lx, lx)
412 real(kind=rp) :: wvr(lx, lx, lx)
413 real(kind=rp) :: wvs(lx, lx, lx)
414 real(kind=rp) :: wvt(lx, lx, lx)
415 real(kind=rp) :: wwr(lx, lx, lx)
416 real(kind=rp) :: wws(lx, lx, lx)
417 real(kind=rp) :: wwt(lx, lx, lx)
418
419 integer :: e, i, j, k, l
420 real(kind=rp) :: dj
421 real(kind=rp) :: s11, s22, s33, s12, s13, s23
422 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
423
424 !$omp do
425 do e = 1, n
426
427 do j = 1, lx * lx
428 do i = 1, lx
429 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
430 + dx(i,2) * u(2,j,1,e) &
431 + dx(i,3) * u(3,j,1,e) &
432 + dx(i,4) * u(4,j,1,e) &
433 + dx(i,5) * u(5,j,1,e) &
434 + dx(i,6) * u(6,j,1,e) &
435 + dx(i,7) * u(7,j,1,e) &
436 + dx(i,8) * u(8,j,1,e) &
437 + dx(i,9) * u(9,j,1,e) &
438 + dx(i,10) * u(10,j,1,e) &
439 + dx(i,11) * u(11,j,1,e) &
440 + dx(i,12) * u(12,j,1,e) &
441 + dx(i,13) * u(13,j,1,e) &
442 + dx(i,14) * u(14,j,1,e)
443
444 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
445 + dx(i,2) * v(2,j,1,e) &
446 + dx(i,3) * v(3,j,1,e) &
447 + dx(i,4) * v(4,j,1,e) &
448 + dx(i,5) * v(5,j,1,e) &
449 + dx(i,6) * v(6,j,1,e) &
450 + dx(i,7) * v(7,j,1,e) &
451 + dx(i,8) * v(8,j,1,e) &
452 + dx(i,9) * v(9,j,1,e) &
453 + dx(i,10) * v(10,j,1,e) &
454 + dx(i,11) * v(11,j,1,e) &
455 + dx(i,12) * v(12,j,1,e) &
456 + dx(i,13) * v(13,j,1,e) &
457 + dx(i,14) * v(14,j,1,e)
458
459 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
460 + dx(i,2) * w(2,j,1,e) &
461 + dx(i,3) * w(3,j,1,e) &
462 + dx(i,4) * w(4,j,1,e) &
463 + dx(i,5) * w(5,j,1,e) &
464 + dx(i,6) * w(6,j,1,e) &
465 + dx(i,7) * w(7,j,1,e) &
466 + dx(i,8) * w(8,j,1,e) &
467 + dx(i,9) * w(9,j,1,e) &
468 + dx(i,10) * w(10,j,1,e) &
469 + dx(i,11) * w(11,j,1,e) &
470 + dx(i,12) * w(12,j,1,e) &
471 + dx(i,13) * w(13,j,1,e) &
472 + dx(i,14) * w(14,j,1,e)
473 end do
474 end do
475
476 do k = 1, lx
477 do j = 1, lx
478 do i = 1, lx
479 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
480 + dy(j,2) * u(i,2,k,e) &
481 + dy(j,3) * u(i,3,k,e) &
482 + dy(j,4) * u(i,4,k,e) &
483 + dy(j,5) * u(i,5,k,e) &
484 + dy(j,6) * u(i,6,k,e) &
485 + dy(j,7) * u(i,7,k,e) &
486 + dy(j,8) * u(i,8,k,e) &
487 + dy(j,9) * u(i,9,k,e) &
488 + dy(j,10) * u(i,10,k,e) &
489 + dy(j,11) * u(i,11,k,e) &
490 + dy(j,12) * u(i,12,k,e) &
491 + dy(j,13) * u(i,13,k,e) &
492 + dy(j,14) * u(i,14,k,e)
493
494
495 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
496 + dy(j,2) * v(i,2,k,e) &
497 + dy(j,3) * v(i,3,k,e) &
498 + dy(j,4) * v(i,4,k,e) &
499 + dy(j,5) * v(i,5,k,e) &
500 + dy(j,6) * v(i,6,k,e) &
501 + dy(j,7) * v(i,7,k,e) &
502 + dy(j,8) * v(i,8,k,e) &
503 + dy(j,9) * v(i,9,k,e) &
504 + dy(j,10) * v(i,10,k,e) &
505 + dy(j,11) * v(i,11,k,e) &
506 + dy(j,12) * v(i,12,k,e) &
507 + dy(j,13) * v(i,13,k,e) &
508 + dy(j,14) * v(i,14,k,e)
509
510 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
511 + dy(j,2) * w(i,2,k,e) &
512 + dy(j,3) * w(i,3,k,e) &
513 + dy(j,4) * w(i,4,k,e) &
514 + dy(j,5) * w(i,5,k,e) &
515 + dy(j,6) * w(i,6,k,e) &
516 + dy(j,7) * w(i,7,k,e) &
517 + dy(j,8) * w(i,8,k,e) &
518 + dy(j,9) * w(i,9,k,e) &
519 + dy(j,10) * w(i,10,k,e) &
520 + dy(j,11) * w(i,11,k,e) &
521 + dy(j,12) * w(i,12,k,e) &
522 + dy(j,13) * w(i,13,k,e) &
523 + dy(j,14) * w(i,14,k,e)
524 end do
525 end do
526 end do
527
528 do k = 1, lx
529 do i = 1, lx*lx
530 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
531 + dz(k,2) * u(i,1,2,e) &
532 + dz(k,3) * u(i,1,3,e) &
533 + dz(k,4) * u(i,1,4,e) &
534 + dz(k,5) * u(i,1,5,e) &
535 + dz(k,6) * u(i,1,6,e) &
536 + dz(k,7) * u(i,1,7,e) &
537 + dz(k,8) * u(i,1,8,e) &
538 + dz(k,9) * u(i,1,9,e) &
539 + dz(k,10) * u(i,1,10,e) &
540 + dz(k,11) * u(i,1,11,e) &
541 + dz(k,12) * u(i,1,12,e) &
542 + dz(k,13) * u(i,1,13,e) &
543 + dz(k,14) * u(i,1,14,e)
544
545 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
546 + dz(k,2) * v(i,1,2,e) &
547 + dz(k,3) * v(i,1,3,e) &
548 + dz(k,4) * v(i,1,4,e) &
549 + dz(k,5) * v(i,1,5,e) &
550 + dz(k,6) * v(i,1,6,e) &
551 + dz(k,7) * v(i,1,7,e) &
552 + dz(k,8) * v(i,1,8,e) &
553 + dz(k,9) * v(i,1,9,e) &
554 + dz(k,10) * v(i,1,10,e) &
555 + dz(k,11) * v(i,1,11,e) &
556 + dz(k,12) * v(i,1,12,e) &
557 + dz(k,13) * v(i,1,13,e) &
558 + dz(k,14) * v(i,1,14,e)
559
560 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
561 + dz(k,2) * w(i,1,2,e) &
562 + dz(k,3) * w(i,1,3,e) &
563 + dz(k,4) * w(i,1,4,e) &
564 + dz(k,5) * w(i,1,5,e) &
565 + dz(k,6) * w(i,1,6,e) &
566 + dz(k,7) * w(i,1,7,e) &
567 + dz(k,8) * w(i,1,8,e) &
568 + dz(k,9) * w(i,1,9,e) &
569 + dz(k,10) * w(i,1,10,e) &
570 + dz(k,11) * w(i,1,11,e) &
571 + dz(k,12) * w(i,1,12,e) &
572 + dz(k,13) * w(i,1,13,e) &
573 + dz(k,14) * w(i,1,14,e)
574 end do
575 end do
576
577 do i = 1, lx*lx*lx
578
579 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
580 + wut(i,1,1) * dtdx(i,1,1,e)
581 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
582 + wut(i,1,1) * dtdy(i,1,1,e)
583 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
584 + wut(i,1,1) * dtdz(i,1,1,e)
585
586 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
587 + wvt(i,1,1) * dtdx(i,1,1,e)
588 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
589 + wvt(i,1,1) * dtdy(i,1,1,e)
590 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
591 + wvt(i,1,1) * dtdz(i,1,1,e)
592
593 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
594 + wwt(i,1,1) * dtdx(i,1,1,e)
595 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
596 + wwt(i,1,1) * dtdy(i,1,1,e)
597 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
598 + wwt(i,1,1) * dtdz(i,1,1,e)
599
600 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
601 s11 = dj*(u1 + u1)
602 s22 = dj*(v2 + v2)
603 s33 = dj*(w3 + w3)
604 s12 = dj*(u2 + v1)
605 s13 = dj*(u3 + w1)
606 s23 = dj*(v3 + w2)
607
608 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
609 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
610 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
611
612 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
613 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
614 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
615
616 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
617 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
618 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
619 end do
620
621 do j = 1, lx*lx
622 do i = 1, lx
623 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
624 + dxt(i,2) * wur(2,j,1) &
625 + dxt(i,3) * wur(3,j,1) &
626 + dxt(i,4) * wur(4,j,1) &
627 + dxt(i,5) * wur(5,j,1) &
628 + dxt(i,6) * wur(6,j,1) &
629 + dxt(i,7) * wur(7,j,1) &
630 + dxt(i,8) * wur(8,j,1) &
631 + dxt(i,9) * wur(9,j,1) &
632 + dxt(i,10) * wur(10,j,1) &
633 + dxt(i,11) * wur(11,j,1) &
634 + dxt(i,12) * wur(12,j,1) &
635 + dxt(i,13) * wur(13,j,1) &
636 + dxt(i,14) * wur(14,j,1)
637
638 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
639 + dxt(i,2) * wvr(2,j,1) &
640 + dxt(i,3) * wvr(3,j,1) &
641 + dxt(i,4) * wvr(4,j,1) &
642 + dxt(i,5) * wvr(5,j,1) &
643 + dxt(i,6) * wvr(6,j,1) &
644 + dxt(i,7) * wvr(7,j,1) &
645 + dxt(i,8) * wvr(8,j,1) &
646 + dxt(i,9) * wvr(9,j,1) &
647 + dxt(i,10) * wvr(10,j,1) &
648 + dxt(i,11) * wvr(11,j,1) &
649 + dxt(i,12) * wvr(12,j,1) &
650 + dxt(i,13) * wvr(13,j,1) &
651 + dxt(i,14) * wvr(14,j,1)
652
653 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
654 + dxt(i,2) * wwr(2,j,1) &
655 + dxt(i,3) * wwr(3,j,1) &
656 + dxt(i,4) * wwr(4,j,1) &
657 + dxt(i,5) * wwr(5,j,1) &
658 + dxt(i,6) * wwr(6,j,1) &
659 + dxt(i,7) * wwr(7,j,1) &
660 + dxt(i,8) * wwr(8,j,1) &
661 + dxt(i,9) * wwr(9,j,1) &
662 + dxt(i,10) * wwr(10,j,1) &
663 + dxt(i,11) * wwr(11,j,1) &
664 + dxt(i,12) * wwr(12,j,1) &
665 + dxt(i,13) * wwr(13,j,1) &
666 + dxt(i,14) * wwr(14,j,1)
667
668 end do
669 end do
670
671 do k = 1, lx
672 do j = 1, lx
673 do i = 1, lx
674 au(i,j,k,e) = au(i,j,k,e) &
675 + dyt(j,1) * wus(i,1,k) &
676 + dyt(j,2) * wus(i,2,k) &
677 + dyt(j,3) * wus(i,3,k) &
678 + dyt(j,4) * wus(i,4,k) &
679 + dyt(j,5) * wus(i,5,k) &
680 + dyt(j,6) * wus(i,6,k) &
681 + dyt(j,7) * wus(i,7,k) &
682 + dyt(j,8) * wus(i,8,k) &
683 + dyt(j,9) * wus(i,9,k) &
684 + dyt(j,10) * wus(i,10,k) &
685 + dyt(j,11) * wus(i,11,k) &
686 + dyt(j,12) * wus(i,12,k) &
687 + dyt(j,13) * wus(i,13,k) &
688 + dyt(j,14) * wus(i,14,k)
689
690 av(i,j,k,e) = av(i,j,k,e) &
691 + dyt(j,1) * wvs(i,1,k) &
692 + dyt(j,2) * wvs(i,2,k) &
693 + dyt(j,3) * wvs(i,3,k) &
694 + dyt(j,4) * wvs(i,4,k) &
695 + dyt(j,5) * wvs(i,5,k) &
696 + dyt(j,6) * wvs(i,6,k) &
697 + dyt(j,7) * wvs(i,7,k) &
698 + dyt(j,8) * wvs(i,8,k) &
699 + dyt(j,9) * wvs(i,9,k) &
700 + dyt(j,10) * wvs(i,10,k) &
701 + dyt(j,11) * wvs(i,11,k) &
702 + dyt(j,12) * wvs(i,12,k) &
703 + dyt(j,13) * wvs(i,13,k) &
704 + dyt(j,14) * wvs(i,14,k)
705
706 aw(i,j,k,e) = aw(i,j,k,e) &
707 + dyt(j,1) * wws(i,1,k) &
708 + dyt(j,2) * wws(i,2,k) &
709 + dyt(j,3) * wws(i,3,k) &
710 + dyt(j,4) * wws(i,4,k) &
711 + dyt(j,5) * wws(i,5,k) &
712 + dyt(j,6) * wws(i,6,k) &
713 + dyt(j,7) * wws(i,7,k) &
714 + dyt(j,8) * wws(i,8,k) &
715 + dyt(j,9) * wws(i,9,k) &
716 + dyt(j,10) * wws(i,10,k) &
717 + dyt(j,11) * wws(i,11,k) &
718 + dyt(j,12) * wws(i,12,k) &
719 + dyt(j,13) * wws(i,13,k) &
720 + dyt(j,14) * wws(i,14,k)
721 end do
722 end do
723 end do
724
725 do k = 1, lx
726 do i = 1, lx*lx
727 au(i,1,k,e) = au(i,1,k,e) &
728 + dzt(k,1) * wut(i,1,1) &
729 + dzt(k,2) * wut(i,1,2) &
730 + dzt(k,3) * wut(i,1,3) &
731 + dzt(k,4) * wut(i,1,4) &
732 + dzt(k,5) * wut(i,1,5) &
733 + dzt(k,6) * wut(i,1,6) &
734 + dzt(k,7) * wut(i,1,7) &
735 + dzt(k,8) * wut(i,1,8) &
736 + dzt(k,9) * wut(i,1,9) &
737 + dzt(k,10) * wut(i,1,10) &
738 + dzt(k,11) * wut(i,1,11) &
739 + dzt(k,12) * wut(i,1,12) &
740 + dzt(k,13) * wut(i,1,13) &
741 + dzt(k,14) * wut(i,1,14)
742
743 av(i,1,k,e) = av(i,1,k,e) &
744 + dzt(k,1) * wvt(i,1,1) &
745 + dzt(k,2) * wvt(i,1,2) &
746 + dzt(k,3) * wvt(i,1,3) &
747 + dzt(k,4) * wvt(i,1,4) &
748 + dzt(k,5) * wvt(i,1,5) &
749 + dzt(k,6) * wvt(i,1,6) &
750 + dzt(k,7) * wvt(i,1,7) &
751 + dzt(k,8) * wvt(i,1,8) &
752 + dzt(k,9) * wvt(i,1,9) &
753 + dzt(k,10) * wvt(i,1,10) &
754 + dzt(k,11) * wvt(i,1,11) &
755 + dzt(k,12) * wvt(i,1,12) &
756 + dzt(k,13) * wvt(i,1,13) &
757 + dzt(k,14) * wvt(i,1,14)
758
759 aw(i,1,k,e) = aw(i,1,k,e) &
760 + dzt(k,1) * wwt(i,1,1) &
761 + dzt(k,2) * wwt(i,1,2) &
762 + dzt(k,3) * wwt(i,1,3) &
763 + dzt(k,4) * wwt(i,1,4) &
764 + dzt(k,5) * wwt(i,1,5) &
765 + dzt(k,6) * wwt(i,1,6) &
766 + dzt(k,7) * wwt(i,1,7) &
767 + dzt(k,8) * wwt(i,1,8) &
768 + dzt(k,9) * wwt(i,1,9) &
769 + dzt(k,10) * wwt(i,1,10) &
770 + dzt(k,11) * wwt(i,1,11) &
771 + dzt(k,12) * wwt(i,1,12) &
772 + dzt(k,13) * wwt(i,1,13) &
773 + dzt(k,14) * wwt(i,1,14)
774 end do
775 end do
776
777 end do
778 !$omp end do
779 end subroutine ax_helm_stress_lx14
780
781 subroutine ax_helm_stress_lx13(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
782 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
783 jacinv, weights3, n)
784 integer, parameter :: lx = 13
785 integer, intent(in) :: n
786 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
787 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
788 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
789 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
790 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
791 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
792 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
793 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
794 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
795 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
796 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
797 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
798 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
799 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
800 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
801 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
802 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
803 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
804 real(kind=rp), intent(in) :: dx(lx, lx)
805 real(kind=rp), intent(in) :: dy(lx, lx)
806 real(kind=rp), intent(in) :: dz(lx, lx)
807 real(kind=rp), intent(in) :: dxt(lx, lx)
808 real(kind=rp), intent(in) :: dyt(lx, lx)
809 real(kind=rp), intent(in) :: dzt(lx, lx)
810 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
811
812 real(kind=rp) :: wur(lx, lx, lx)
813 real(kind=rp) :: wus(lx, lx, lx)
814 real(kind=rp) :: wut(lx, lx, lx)
815 real(kind=rp) :: wvr(lx, lx, lx)
816 real(kind=rp) :: wvs(lx, lx, lx)
817 real(kind=rp) :: wvt(lx, lx, lx)
818 real(kind=rp) :: wwr(lx, lx, lx)
819 real(kind=rp) :: wws(lx, lx, lx)
820 real(kind=rp) :: wwt(lx, lx, lx)
821
822 integer :: e, i, j, k, l
823 real(kind=rp) :: dj
824 real(kind=rp) :: s11, s22, s33, s12, s13, s23
825 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
826
827 !$omp do
828 do e = 1, n
829
830 do j = 1, lx * lx
831 do i = 1, lx
832 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
833 + dx(i,2) * u(2,j,1,e) &
834 + dx(i,3) * u(3,j,1,e) &
835 + dx(i,4) * u(4,j,1,e) &
836 + dx(i,5) * u(5,j,1,e) &
837 + dx(i,6) * u(6,j,1,e) &
838 + dx(i,7) * u(7,j,1,e) &
839 + dx(i,8) * u(8,j,1,e) &
840 + dx(i,9) * u(9,j,1,e) &
841 + dx(i,10) * u(10,j,1,e) &
842 + dx(i,11) * u(11,j,1,e) &
843 + dx(i,12) * u(12,j,1,e) &
844 + dx(i,13) * u(13,j,1,e)
845
846 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
847 + dx(i,2) * v(2,j,1,e) &
848 + dx(i,3) * v(3,j,1,e) &
849 + dx(i,4) * v(4,j,1,e) &
850 + dx(i,5) * v(5,j,1,e) &
851 + dx(i,6) * v(6,j,1,e) &
852 + dx(i,7) * v(7,j,1,e) &
853 + dx(i,8) * v(8,j,1,e) &
854 + dx(i,9) * v(9,j,1,e) &
855 + dx(i,10) * v(10,j,1,e) &
856 + dx(i,11) * v(11,j,1,e) &
857 + dx(i,12) * v(12,j,1,e) &
858 + dx(i,13) * v(13,j,1,e)
859
860 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
861 + dx(i,2) * w(2,j,1,e) &
862 + dx(i,3) * w(3,j,1,e) &
863 + dx(i,4) * w(4,j,1,e) &
864 + dx(i,5) * w(5,j,1,e) &
865 + dx(i,6) * w(6,j,1,e) &
866 + dx(i,7) * w(7,j,1,e) &
867 + dx(i,8) * w(8,j,1,e) &
868 + dx(i,9) * w(9,j,1,e) &
869 + dx(i,10) * w(10,j,1,e) &
870 + dx(i,11) * w(11,j,1,e) &
871 + dx(i,12) * w(12,j,1,e) &
872 + dx(i,13) * w(13,j,1,e)
873 end do
874 end do
875
876 do k = 1, lx
877 do j = 1, lx
878 do i = 1, lx
879 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
880 + dy(j,2) * u(i,2,k,e) &
881 + dy(j,3) * u(i,3,k,e) &
882 + dy(j,4) * u(i,4,k,e) &
883 + dy(j,5) * u(i,5,k,e) &
884 + dy(j,6) * u(i,6,k,e) &
885 + dy(j,7) * u(i,7,k,e) &
886 + dy(j,8) * u(i,8,k,e) &
887 + dy(j,9) * u(i,9,k,e) &
888 + dy(j,10) * u(i,10,k,e) &
889 + dy(j,11) * u(i,11,k,e) &
890 + dy(j,12) * u(i,12,k,e) &
891 + dy(j,13) * u(i,13,k,e)
892
893
894 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
895 + dy(j,2) * v(i,2,k,e) &
896 + dy(j,3) * v(i,3,k,e) &
897 + dy(j,4) * v(i,4,k,e) &
898 + dy(j,5) * v(i,5,k,e) &
899 + dy(j,6) * v(i,6,k,e) &
900 + dy(j,7) * v(i,7,k,e) &
901 + dy(j,8) * v(i,8,k,e) &
902 + dy(j,9) * v(i,9,k,e) &
903 + dy(j,10) * v(i,10,k,e) &
904 + dy(j,11) * v(i,11,k,e) &
905 + dy(j,12) * v(i,12,k,e) &
906 + dy(j,13) * v(i,13,k,e)
907
908 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
909 + dy(j,2) * w(i,2,k,e) &
910 + dy(j,3) * w(i,3,k,e) &
911 + dy(j,4) * w(i,4,k,e) &
912 + dy(j,5) * w(i,5,k,e) &
913 + dy(j,6) * w(i,6,k,e) &
914 + dy(j,7) * w(i,7,k,e) &
915 + dy(j,8) * w(i,8,k,e) &
916 + dy(j,9) * w(i,9,k,e) &
917 + dy(j,10) * w(i,10,k,e) &
918 + dy(j,11) * w(i,11,k,e) &
919 + dy(j,12) * w(i,12,k,e) &
920 + dy(j,13) * w(i,13,k,e)
921 end do
922 end do
923 end do
924
925 do k = 1, lx
926 do i = 1, lx*lx
927 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
928 + dz(k,2) * u(i,1,2,e) &
929 + dz(k,3) * u(i,1,3,e) &
930 + dz(k,4) * u(i,1,4,e) &
931 + dz(k,5) * u(i,1,5,e) &
932 + dz(k,6) * u(i,1,6,e) &
933 + dz(k,7) * u(i,1,7,e) &
934 + dz(k,8) * u(i,1,8,e) &
935 + dz(k,9) * u(i,1,9,e) &
936 + dz(k,10) * u(i,1,10,e) &
937 + dz(k,11) * u(i,1,11,e) &
938 + dz(k,12) * u(i,1,12,e) &
939 + dz(k,13) * u(i,1,13,e)
940
941 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
942 + dz(k,2) * v(i,1,2,e) &
943 + dz(k,3) * v(i,1,3,e) &
944 + dz(k,4) * v(i,1,4,e) &
945 + dz(k,5) * v(i,1,5,e) &
946 + dz(k,6) * v(i,1,6,e) &
947 + dz(k,7) * v(i,1,7,e) &
948 + dz(k,8) * v(i,1,8,e) &
949 + dz(k,9) * v(i,1,9,e) &
950 + dz(k,10) * v(i,1,10,e) &
951 + dz(k,11) * v(i,1,11,e) &
952 + dz(k,12) * v(i,1,12,e) &
953 + dz(k,13) * v(i,1,13,e)
954
955 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
956 + dz(k,2) * w(i,1,2,e) &
957 + dz(k,3) * w(i,1,3,e) &
958 + dz(k,4) * w(i,1,4,e) &
959 + dz(k,5) * w(i,1,5,e) &
960 + dz(k,6) * w(i,1,6,e) &
961 + dz(k,7) * w(i,1,7,e) &
962 + dz(k,8) * w(i,1,8,e) &
963 + dz(k,9) * w(i,1,9,e) &
964 + dz(k,10) * w(i,1,10,e) &
965 + dz(k,11) * w(i,1,11,e) &
966 + dz(k,12) * w(i,1,12,e) &
967 + dz(k,13) * w(i,1,13,e)
968 end do
969 end do
970
971 do i = 1, lx*lx*lx
972
973 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
974 + wut(i,1,1) * dtdx(i,1,1,e)
975 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
976 + wut(i,1,1) * dtdy(i,1,1,e)
977 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
978 + wut(i,1,1) * dtdz(i,1,1,e)
979
980 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
981 + wvt(i,1,1) * dtdx(i,1,1,e)
982 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
983 + wvt(i,1,1) * dtdy(i,1,1,e)
984 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
985 + wvt(i,1,1) * dtdz(i,1,1,e)
986
987 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
988 + wwt(i,1,1) * dtdx(i,1,1,e)
989 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
990 + wwt(i,1,1) * dtdy(i,1,1,e)
991 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
992 + wwt(i,1,1) * dtdz(i,1,1,e)
993
994 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
995 s11 = dj*(u1 + u1)
996 s22 = dj*(v2 + v2)
997 s33 = dj*(w3 + w3)
998 s12 = dj*(u2 + v1)
999 s13 = dj*(u3 + w1)
1000 s23 = dj*(v3 + w2)
1001
1002 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
1003 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
1004 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
1005
1006 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
1007 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
1008 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
1009
1010 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
1011 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
1012 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
1013 end do
1014
1015 do j = 1, lx*lx
1016 do i = 1, lx
1017 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
1018 + dxt(i,2) * wur(2,j,1) &
1019 + dxt(i,3) * wur(3,j,1) &
1020 + dxt(i,4) * wur(4,j,1) &
1021 + dxt(i,5) * wur(5,j,1) &
1022 + dxt(i,6) * wur(6,j,1) &
1023 + dxt(i,7) * wur(7,j,1) &
1024 + dxt(i,8) * wur(8,j,1) &
1025 + dxt(i,9) * wur(9,j,1) &
1026 + dxt(i,10) * wur(10,j,1) &
1027 + dxt(i,11) * wur(11,j,1) &
1028 + dxt(i,12) * wur(12,j,1) &
1029 + dxt(i,13) * wur(13,j,1)
1030
1031 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
1032 + dxt(i,2) * wvr(2,j,1) &
1033 + dxt(i,3) * wvr(3,j,1) &
1034 + dxt(i,4) * wvr(4,j,1) &
1035 + dxt(i,5) * wvr(5,j,1) &
1036 + dxt(i,6) * wvr(6,j,1) &
1037 + dxt(i,7) * wvr(7,j,1) &
1038 + dxt(i,8) * wvr(8,j,1) &
1039 + dxt(i,9) * wvr(9,j,1) &
1040 + dxt(i,10) * wvr(10,j,1) &
1041 + dxt(i,11) * wvr(11,j,1) &
1042 + dxt(i,12) * wvr(12,j,1) &
1043 + dxt(i,13) * wvr(13,j,1)
1044
1045 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
1046 + dxt(i,2) * wwr(2,j,1) &
1047 + dxt(i,3) * wwr(3,j,1) &
1048 + dxt(i,4) * wwr(4,j,1) &
1049 + dxt(i,5) * wwr(5,j,1) &
1050 + dxt(i,6) * wwr(6,j,1) &
1051 + dxt(i,7) * wwr(7,j,1) &
1052 + dxt(i,8) * wwr(8,j,1) &
1053 + dxt(i,9) * wwr(9,j,1) &
1054 + dxt(i,10) * wwr(10,j,1) &
1055 + dxt(i,11) * wwr(11,j,1) &
1056 + dxt(i,12) * wwr(12,j,1) &
1057 + dxt(i,13) * wwr(13,j,1)
1058 end do
1059 end do
1060
1061 do k = 1, lx
1062 do j = 1, lx
1063 do i = 1, lx
1064 au(i,j,k,e) = au(i,j,k,e) &
1065 + dyt(j,1) * wus(i,1,k) &
1066 + dyt(j,2) * wus(i,2,k) &
1067 + dyt(j,3) * wus(i,3,k) &
1068 + dyt(j,4) * wus(i,4,k) &
1069 + dyt(j,5) * wus(i,5,k) &
1070 + dyt(j,6) * wus(i,6,k) &
1071 + dyt(j,7) * wus(i,7,k) &
1072 + dyt(j,8) * wus(i,8,k) &
1073 + dyt(j,9) * wus(i,9,k) &
1074 + dyt(j,10) * wus(i,10,k) &
1075 + dyt(j,11) * wus(i,11,k) &
1076 + dyt(j,12) * wus(i,12,k) &
1077 + dyt(j,13) * wus(i,13,k)
1078
1079 av(i,j,k,e) = av(i,j,k,e) &
1080 + dyt(j,1) * wvs(i,1,k) &
1081 + dyt(j,2) * wvs(i,2,k) &
1082 + dyt(j,3) * wvs(i,3,k) &
1083 + dyt(j,4) * wvs(i,4,k) &
1084 + dyt(j,5) * wvs(i,5,k) &
1085 + dyt(j,6) * wvs(i,6,k) &
1086 + dyt(j,7) * wvs(i,7,k) &
1087 + dyt(j,8) * wvs(i,8,k) &
1088 + dyt(j,9) * wvs(i,9,k) &
1089 + dyt(j,10) * wvs(i,10,k) &
1090 + dyt(j,11) * wvs(i,11,k) &
1091 + dyt(j,12) * wvs(i,12,k) &
1092 + dyt(j,13) * wvs(i,13,k)
1093
1094 aw(i,j,k,e) = aw(i,j,k,e) &
1095 + dyt(j,1) * wws(i,1,k) &
1096 + dyt(j,2) * wws(i,2,k) &
1097 + dyt(j,3) * wws(i,3,k) &
1098 + dyt(j,4) * wws(i,4,k) &
1099 + dyt(j,5) * wws(i,5,k) &
1100 + dyt(j,6) * wws(i,6,k) &
1101 + dyt(j,7) * wws(i,7,k) &
1102 + dyt(j,8) * wws(i,8,k) &
1103 + dyt(j,9) * wws(i,9,k) &
1104 + dyt(j,10) * wws(i,10,k) &
1105 + dyt(j,11) * wws(i,11,k) &
1106 + dyt(j,12) * wws(i,12,k) &
1107 + dyt(j,13) * wws(i,13,k)
1108 end do
1109 end do
1110 end do
1111
1112 do k = 1, lx
1113 do i = 1, lx*lx
1114 au(i,1,k,e) = au(i,1,k,e) &
1115 + dzt(k,1) * wut(i,1,1) &
1116 + dzt(k,2) * wut(i,1,2) &
1117 + dzt(k,3) * wut(i,1,3) &
1118 + dzt(k,4) * wut(i,1,4) &
1119 + dzt(k,5) * wut(i,1,5) &
1120 + dzt(k,6) * wut(i,1,6) &
1121 + dzt(k,7) * wut(i,1,7) &
1122 + dzt(k,8) * wut(i,1,8) &
1123 + dzt(k,9) * wut(i,1,9) &
1124 + dzt(k,10) * wut(i,1,10) &
1125 + dzt(k,11) * wut(i,1,11) &
1126 + dzt(k,12) * wut(i,1,12) &
1127 + dzt(k,13) * wut(i,1,13)
1128
1129 av(i,1,k,e) = av(i,1,k,e) &
1130 + dzt(k,1) * wvt(i,1,1) &
1131 + dzt(k,2) * wvt(i,1,2) &
1132 + dzt(k,3) * wvt(i,1,3) &
1133 + dzt(k,4) * wvt(i,1,4) &
1134 + dzt(k,5) * wvt(i,1,5) &
1135 + dzt(k,6) * wvt(i,1,6) &
1136 + dzt(k,7) * wvt(i,1,7) &
1137 + dzt(k,8) * wvt(i,1,8) &
1138 + dzt(k,9) * wvt(i,1,9) &
1139 + dzt(k,10) * wvt(i,1,10) &
1140 + dzt(k,11) * wvt(i,1,11) &
1141 + dzt(k,12) * wvt(i,1,12) &
1142 + dzt(k,13) * wvt(i,1,13)
1143
1144 aw(i,1,k,e) = aw(i,1,k,e) &
1145 + dzt(k,1) * wwt(i,1,1) &
1146 + dzt(k,2) * wwt(i,1,2) &
1147 + dzt(k,3) * wwt(i,1,3) &
1148 + dzt(k,4) * wwt(i,1,4) &
1149 + dzt(k,5) * wwt(i,1,5) &
1150 + dzt(k,6) * wwt(i,1,6) &
1151 + dzt(k,7) * wwt(i,1,7) &
1152 + dzt(k,8) * wwt(i,1,8) &
1153 + dzt(k,9) * wwt(i,1,9) &
1154 + dzt(k,10) * wwt(i,1,10) &
1155 + dzt(k,11) * wwt(i,1,11) &
1156 + dzt(k,12) * wwt(i,1,12) &
1157 + dzt(k,13) * wwt(i,1,13)
1158 end do
1159 end do
1160
1161 end do
1162 !$omp end do
1163 end subroutine ax_helm_stress_lx13
1164
1165 subroutine ax_helm_stress_lx12(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
1166 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
1167 jacinv, weights3, n)
1168 integer, parameter :: lx = 12
1169 integer, intent(in) :: n
1170 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1171 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
1172 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
1173 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
1174 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
1175 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
1176 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1177 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
1178 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
1179 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
1180 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
1181 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
1182 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
1183 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
1184 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
1185 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
1186 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
1187 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
1188 real(kind=rp), intent(in) :: dx(lx, lx)
1189 real(kind=rp), intent(in) :: dy(lx, lx)
1190 real(kind=rp), intent(in) :: dz(lx, lx)
1191 real(kind=rp), intent(in) :: dxt(lx, lx)
1192 real(kind=rp), intent(in) :: dyt(lx, lx)
1193 real(kind=rp), intent(in) :: dzt(lx, lx)
1194 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
1195
1196 real(kind=rp) :: wur(lx, lx, lx)
1197 real(kind=rp) :: wus(lx, lx, lx)
1198 real(kind=rp) :: wut(lx, lx, lx)
1199 real(kind=rp) :: wvr(lx, lx, lx)
1200 real(kind=rp) :: wvs(lx, lx, lx)
1201 real(kind=rp) :: wvt(lx, lx, lx)
1202 real(kind=rp) :: wwr(lx, lx, lx)
1203 real(kind=rp) :: wws(lx, lx, lx)
1204 real(kind=rp) :: wwt(lx, lx, lx)
1205
1206 integer :: e, i, j, k, l
1207 real(kind=rp) :: dj
1208 real(kind=rp) :: s11, s22, s33, s12, s13, s23
1209 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
1210
1211 !$omp do
1212 do e = 1, n
1213
1214 do j = 1, lx * lx
1215 do i = 1, lx
1216 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1217 + dx(i,2) * u(2,j,1,e) &
1218 + dx(i,3) * u(3,j,1,e) &
1219 + dx(i,4) * u(4,j,1,e) &
1220 + dx(i,5) * u(5,j,1,e) &
1221 + dx(i,6) * u(6,j,1,e) &
1222 + dx(i,7) * u(7,j,1,e) &
1223 + dx(i,8) * u(8,j,1,e) &
1224 + dx(i,9) * u(9,j,1,e) &
1225 + dx(i,10) * u(10,j,1,e) &
1226 + dx(i,11) * u(11,j,1,e) &
1227 + dx(i,12) * u(12,j,1,e)
1228
1229 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
1230 + dx(i,2) * v(2,j,1,e) &
1231 + dx(i,3) * v(3,j,1,e) &
1232 + dx(i,4) * v(4,j,1,e) &
1233 + dx(i,5) * v(5,j,1,e) &
1234 + dx(i,6) * v(6,j,1,e) &
1235 + dx(i,7) * v(7,j,1,e) &
1236 + dx(i,8) * v(8,j,1,e) &
1237 + dx(i,9) * v(9,j,1,e) &
1238 + dx(i,10) * v(10,j,1,e) &
1239 + dx(i,11) * v(11,j,1,e) &
1240 + dx(i,12) * v(12,j,1,e)
1241
1242 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
1243 + dx(i,2) * w(2,j,1,e) &
1244 + dx(i,3) * w(3,j,1,e) &
1245 + dx(i,4) * w(4,j,1,e) &
1246 + dx(i,5) * w(5,j,1,e) &
1247 + dx(i,6) * w(6,j,1,e) &
1248 + dx(i,7) * w(7,j,1,e) &
1249 + dx(i,8) * w(8,j,1,e) &
1250 + dx(i,9) * w(9,j,1,e) &
1251 + dx(i,10) * w(10,j,1,e) &
1252 + dx(i,11) * w(11,j,1,e) &
1253 + dx(i,12) * w(12,j,1,e)
1254 end do
1255 end do
1256
1257 do k = 1, lx
1258 do j = 1, lx
1259 do i = 1, lx
1260 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1261 + dy(j,2) * u(i,2,k,e) &
1262 + dy(j,3) * u(i,3,k,e) &
1263 + dy(j,4) * u(i,4,k,e) &
1264 + dy(j,5) * u(i,5,k,e) &
1265 + dy(j,6) * u(i,6,k,e) &
1266 + dy(j,7) * u(i,7,k,e) &
1267 + dy(j,8) * u(i,8,k,e) &
1268 + dy(j,9) * u(i,9,k,e) &
1269 + dy(j,10) * u(i,10,k,e) &
1270 + dy(j,11) * u(i,11,k,e) &
1271 + dy(j,12) * u(i,12,k,e)
1272
1273
1274 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
1275 + dy(j,2) * v(i,2,k,e) &
1276 + dy(j,3) * v(i,3,k,e) &
1277 + dy(j,4) * v(i,4,k,e) &
1278 + dy(j,5) * v(i,5,k,e) &
1279 + dy(j,6) * v(i,6,k,e) &
1280 + dy(j,7) * v(i,7,k,e) &
1281 + dy(j,8) * v(i,8,k,e) &
1282 + dy(j,9) * v(i,9,k,e) &
1283 + dy(j,10) * v(i,10,k,e) &
1284 + dy(j,11) * v(i,11,k,e) &
1285 + dy(j,12) * v(i,12,k,e)
1286
1287 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
1288 + dy(j,2) * w(i,2,k,e) &
1289 + dy(j,3) * w(i,3,k,e) &
1290 + dy(j,4) * w(i,4,k,e) &
1291 + dy(j,5) * w(i,5,k,e) &
1292 + dy(j,6) * w(i,6,k,e) &
1293 + dy(j,7) * w(i,7,k,e) &
1294 + dy(j,8) * w(i,8,k,e) &
1295 + dy(j,9) * w(i,9,k,e) &
1296 + dy(j,10) * w(i,10,k,e) &
1297 + dy(j,11) * w(i,11,k,e) &
1298 + dy(j,12) * w(i,12,k,e)
1299 end do
1300 end do
1301 end do
1302
1303 do k = 1, lx
1304 do i = 1, lx*lx
1305 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1306 + dz(k,2) * u(i,1,2,e) &
1307 + dz(k,3) * u(i,1,3,e) &
1308 + dz(k,4) * u(i,1,4,e) &
1309 + dz(k,5) * u(i,1,5,e) &
1310 + dz(k,6) * u(i,1,6,e) &
1311 + dz(k,7) * u(i,1,7,e) &
1312 + dz(k,8) * u(i,1,8,e) &
1313 + dz(k,9) * u(i,1,9,e) &
1314 + dz(k,10) * u(i,1,10,e) &
1315 + dz(k,11) * u(i,1,11,e) &
1316 + dz(k,12) * u(i,1,12,e)
1317
1318 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
1319 + dz(k,2) * v(i,1,2,e) &
1320 + dz(k,3) * v(i,1,3,e) &
1321 + dz(k,4) * v(i,1,4,e) &
1322 + dz(k,5) * v(i,1,5,e) &
1323 + dz(k,6) * v(i,1,6,e) &
1324 + dz(k,7) * v(i,1,7,e) &
1325 + dz(k,8) * v(i,1,8,e) &
1326 + dz(k,9) * v(i,1,9,e) &
1327 + dz(k,10) * v(i,1,10,e) &
1328 + dz(k,11) * v(i,1,11,e) &
1329 + dz(k,12) * v(i,1,12,e)
1330
1331 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
1332 + dz(k,2) * w(i,1,2,e) &
1333 + dz(k,3) * w(i,1,3,e) &
1334 + dz(k,4) * w(i,1,4,e) &
1335 + dz(k,5) * w(i,1,5,e) &
1336 + dz(k,6) * w(i,1,6,e) &
1337 + dz(k,7) * w(i,1,7,e) &
1338 + dz(k,8) * w(i,1,8,e) &
1339 + dz(k,9) * w(i,1,9,e) &
1340 + dz(k,10) * w(i,1,10,e) &
1341 + dz(k,11) * w(i,1,11,e) &
1342 + dz(k,12) * w(i,1,12,e)
1343 end do
1344 end do
1345
1346 do i = 1, lx*lx*lx
1347
1348 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
1349 + wut(i,1,1) * dtdx(i,1,1,e)
1350 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
1351 + wut(i,1,1) * dtdy(i,1,1,e)
1352 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
1353 + wut(i,1,1) * dtdz(i,1,1,e)
1354
1355 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
1356 + wvt(i,1,1) * dtdx(i,1,1,e)
1357 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
1358 + wvt(i,1,1) * dtdy(i,1,1,e)
1359 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
1360 + wvt(i,1,1) * dtdz(i,1,1,e)
1361
1362 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
1363 + wwt(i,1,1) * dtdx(i,1,1,e)
1364 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
1365 + wwt(i,1,1) * dtdy(i,1,1,e)
1366 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
1367 + wwt(i,1,1) * dtdz(i,1,1,e)
1368
1369 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
1370 s11 = dj*(u1 + u1)
1371 s22 = dj*(v2 + v2)
1372 s33 = dj*(w3 + w3)
1373 s12 = dj*(u2 + v1)
1374 s13 = dj*(u3 + w1)
1375 s23 = dj*(v3 + w2)
1376
1377 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
1378 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
1379 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
1380
1381 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
1382 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
1383 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
1384
1385 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
1386 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
1387 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
1388 end do
1389
1390 do j = 1, lx*lx
1391 do i = 1, lx
1392 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
1393 + dxt(i,2) * wur(2,j,1) &
1394 + dxt(i,3) * wur(3,j,1) &
1395 + dxt(i,4) * wur(4,j,1) &
1396 + dxt(i,5) * wur(5,j,1) &
1397 + dxt(i,6) * wur(6,j,1) &
1398 + dxt(i,7) * wur(7,j,1) &
1399 + dxt(i,8) * wur(8,j,1) &
1400 + dxt(i,9) * wur(9,j,1) &
1401 + dxt(i,10) * wur(10,j,1) &
1402 + dxt(i,11) * wur(11,j,1) &
1403 + dxt(i,12) * wur(12,j,1)
1404
1405 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
1406 + dxt(i,2) * wvr(2,j,1) &
1407 + dxt(i,3) * wvr(3,j,1) &
1408 + dxt(i,4) * wvr(4,j,1) &
1409 + dxt(i,5) * wvr(5,j,1) &
1410 + dxt(i,6) * wvr(6,j,1) &
1411 + dxt(i,7) * wvr(7,j,1) &
1412 + dxt(i,8) * wvr(8,j,1) &
1413 + dxt(i,9) * wvr(9,j,1) &
1414 + dxt(i,10) * wvr(10,j,1) &
1415 + dxt(i,11) * wvr(11,j,1) &
1416 + dxt(i,12) * wvr(12,j,1)
1417
1418 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
1419 + dxt(i,2) * wwr(2,j,1) &
1420 + dxt(i,3) * wwr(3,j,1) &
1421 + dxt(i,4) * wwr(4,j,1) &
1422 + dxt(i,5) * wwr(5,j,1) &
1423 + dxt(i,6) * wwr(6,j,1) &
1424 + dxt(i,7) * wwr(7,j,1) &
1425 + dxt(i,8) * wwr(8,j,1) &
1426 + dxt(i,9) * wwr(9,j,1) &
1427 + dxt(i,10) * wwr(10,j,1) &
1428 + dxt(i,11) * wwr(11,j,1) &
1429 + dxt(i,12) * wwr(12,j,1)
1430 end do
1431 end do
1432
1433 do k = 1, lx
1434 do j = 1, lx
1435 do i = 1, lx
1436 au(i,j,k,e) = au(i,j,k,e) &
1437 + dyt(j,1) * wus(i,1,k) &
1438 + dyt(j,2) * wus(i,2,k) &
1439 + dyt(j,3) * wus(i,3,k) &
1440 + dyt(j,4) * wus(i,4,k) &
1441 + dyt(j,5) * wus(i,5,k) &
1442 + dyt(j,6) * wus(i,6,k) &
1443 + dyt(j,7) * wus(i,7,k) &
1444 + dyt(j,8) * wus(i,8,k) &
1445 + dyt(j,9) * wus(i,9,k) &
1446 + dyt(j,10) * wus(i,10,k) &
1447 + dyt(j,11) * wus(i,11,k) &
1448 + dyt(j,12) * wus(i,12,k)
1449
1450 av(i,j,k,e) = av(i,j,k,e) &
1451 + dyt(j,1) * wvs(i,1,k) &
1452 + dyt(j,2) * wvs(i,2,k) &
1453 + dyt(j,3) * wvs(i,3,k) &
1454 + dyt(j,4) * wvs(i,4,k) &
1455 + dyt(j,5) * wvs(i,5,k) &
1456 + dyt(j,6) * wvs(i,6,k) &
1457 + dyt(j,7) * wvs(i,7,k) &
1458 + dyt(j,8) * wvs(i,8,k) &
1459 + dyt(j,9) * wvs(i,9,k) &
1460 + dyt(j,10) * wvs(i,10,k) &
1461 + dyt(j,11) * wvs(i,11,k) &
1462 + dyt(j,12) * wvs(i,12,k)
1463
1464 aw(i,j,k,e) = aw(i,j,k,e) &
1465 + dyt(j,1) * wws(i,1,k) &
1466 + dyt(j,2) * wws(i,2,k) &
1467 + dyt(j,3) * wws(i,3,k) &
1468 + dyt(j,4) * wws(i,4,k) &
1469 + dyt(j,5) * wws(i,5,k) &
1470 + dyt(j,6) * wws(i,6,k) &
1471 + dyt(j,7) * wws(i,7,k) &
1472 + dyt(j,8) * wws(i,8,k) &
1473 + dyt(j,9) * wws(i,9,k) &
1474 + dyt(j,10) * wws(i,10,k) &
1475 + dyt(j,11) * wws(i,11,k) &
1476 + dyt(j,12) * wws(i,12,k)
1477 end do
1478 end do
1479 end do
1480
1481 do k = 1, lx
1482 do i = 1, lx*lx
1483 au(i,1,k,e) = au(i,1,k,e) &
1484 + dzt(k,1) * wut(i,1,1) &
1485 + dzt(k,2) * wut(i,1,2) &
1486 + dzt(k,3) * wut(i,1,3) &
1487 + dzt(k,4) * wut(i,1,4) &
1488 + dzt(k,5) * wut(i,1,5) &
1489 + dzt(k,6) * wut(i,1,6) &
1490 + dzt(k,7) * wut(i,1,7) &
1491 + dzt(k,8) * wut(i,1,8) &
1492 + dzt(k,9) * wut(i,1,9) &
1493 + dzt(k,10) * wut(i,1,10) &
1494 + dzt(k,11) * wut(i,1,11) &
1495 + dzt(k,12) * wut(i,1,12)
1496
1497 av(i,1,k,e) = av(i,1,k,e) &
1498 + dzt(k,1) * wvt(i,1,1) &
1499 + dzt(k,2) * wvt(i,1,2) &
1500 + dzt(k,3) * wvt(i,1,3) &
1501 + dzt(k,4) * wvt(i,1,4) &
1502 + dzt(k,5) * wvt(i,1,5) &
1503 + dzt(k,6) * wvt(i,1,6) &
1504 + dzt(k,7) * wvt(i,1,7) &
1505 + dzt(k,8) * wvt(i,1,8) &
1506 + dzt(k,9) * wvt(i,1,9) &
1507 + dzt(k,10) * wvt(i,1,10) &
1508 + dzt(k,11) * wvt(i,1,11) &
1509 + dzt(k,12) * wvt(i,1,12)
1510
1511 aw(i,1,k,e) = aw(i,1,k,e) &
1512 + dzt(k,1) * wwt(i,1,1) &
1513 + dzt(k,2) * wwt(i,1,2) &
1514 + dzt(k,3) * wwt(i,1,3) &
1515 + dzt(k,4) * wwt(i,1,4) &
1516 + dzt(k,5) * wwt(i,1,5) &
1517 + dzt(k,6) * wwt(i,1,6) &
1518 + dzt(k,7) * wwt(i,1,7) &
1519 + dzt(k,8) * wwt(i,1,8) &
1520 + dzt(k,9) * wwt(i,1,9) &
1521 + dzt(k,10) * wwt(i,1,10) &
1522 + dzt(k,11) * wwt(i,1,11) &
1523 + dzt(k,12) * wwt(i,1,12)
1524 end do
1525 end do
1526
1527 end do
1528 !$omp end do
1529 end subroutine ax_helm_stress_lx12
1530
1531 subroutine ax_helm_stress_lx11(au, av, aw, u, v, w, Dx, Dy, Dz, &
1532 Dxt, Dyt, Dzt, &
1533 h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
1534 jacinv, weights3, n)
1535 integer, parameter :: lx = 11
1536 integer, intent(in) :: n
1537 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1538 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
1539 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
1540 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
1541 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
1542 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
1543 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1544 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
1545 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
1546 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
1547 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
1548 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
1549 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
1550 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
1551 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
1552 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
1553 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
1554 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
1555 real(kind=rp), intent(in) :: dx(lx, lx)
1556 real(kind=rp), intent(in) :: dy(lx, lx)
1557 real(kind=rp), intent(in) :: dz(lx, lx)
1558 real(kind=rp), intent(in) :: dxt(lx, lx)
1559 real(kind=rp), intent(in) :: dyt(lx, lx)
1560 real(kind=rp), intent(in) :: dzt(lx, lx)
1561 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
1562
1563 real(kind=rp) :: wur(lx, lx, lx)
1564 real(kind=rp) :: wus(lx, lx, lx)
1565 real(kind=rp) :: wut(lx, lx, lx)
1566 real(kind=rp) :: wvr(lx, lx, lx)
1567 real(kind=rp) :: wvs(lx, lx, lx)
1568 real(kind=rp) :: wvt(lx, lx, lx)
1569 real(kind=rp) :: wwr(lx, lx, lx)
1570 real(kind=rp) :: wws(lx, lx, lx)
1571 real(kind=rp) :: wwt(lx, lx, lx)
1572
1573 integer :: e, i, j, k, l
1574 real(kind=rp) :: dj
1575 real(kind=rp) :: s11, s22, s33, s12, s13, s23
1576 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
1577
1578 !$omp do
1579 do e = 1, n
1580
1581 do j = 1, lx * lx
1582 do i = 1, lx
1583 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1584 + dx(i,2) * u(2,j,1,e) &
1585 + dx(i,3) * u(3,j,1,e) &
1586 + dx(i,4) * u(4,j,1,e) &
1587 + dx(i,5) * u(5,j,1,e) &
1588 + dx(i,6) * u(6,j,1,e) &
1589 + dx(i,7) * u(7,j,1,e) &
1590 + dx(i,8) * u(8,j,1,e) &
1591 + dx(i,9) * u(9,j,1,e) &
1592 + dx(i,10) * u(10,j,1,e) &
1593 + dx(i,11) * u(11,j,1,e)
1594
1595 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
1596 + dx(i,2) * v(2,j,1,e) &
1597 + dx(i,3) * v(3,j,1,e) &
1598 + dx(i,4) * v(4,j,1,e) &
1599 + dx(i,5) * v(5,j,1,e) &
1600 + dx(i,6) * v(6,j,1,e) &
1601 + dx(i,7) * v(7,j,1,e) &
1602 + dx(i,8) * v(8,j,1,e) &
1603 + dx(i,9) * v(9,j,1,e) &
1604 + dx(i,10) * v(10,j,1,e) &
1605 + dx(i,11) * v(11,j,1,e)
1606
1607 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
1608 + dx(i,2) * w(2,j,1,e) &
1609 + dx(i,3) * w(3,j,1,e) &
1610 + dx(i,4) * w(4,j,1,e) &
1611 + dx(i,5) * w(5,j,1,e) &
1612 + dx(i,6) * w(6,j,1,e) &
1613 + dx(i,7) * w(7,j,1,e) &
1614 + dx(i,8) * w(8,j,1,e) &
1615 + dx(i,9) * w(9,j,1,e) &
1616 + dx(i,10) * w(10,j,1,e) &
1617 + dx(i,11) * w(11,j,1,e)
1618 end do
1619 end do
1620
1621 do k = 1, lx
1622 do j = 1, lx
1623 do i = 1, lx
1624 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1625 + dy(j,2) * u(i,2,k,e) &
1626 + dy(j,3) * u(i,3,k,e) &
1627 + dy(j,4) * u(i,4,k,e) &
1628 + dy(j,5) * u(i,5,k,e) &
1629 + dy(j,6) * u(i,6,k,e) &
1630 + dy(j,7) * u(i,7,k,e) &
1631 + dy(j,8) * u(i,8,k,e) &
1632 + dy(j,9) * u(i,9,k,e) &
1633 + dy(j,10) * u(i,10,k,e) &
1634 + dy(j,11) * u(i,11,k,e)
1635
1636 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
1637 + dy(j,2) * v(i,2,k,e) &
1638 + dy(j,3) * v(i,3,k,e) &
1639 + dy(j,4) * v(i,4,k,e) &
1640 + dy(j,5) * v(i,5,k,e) &
1641 + dy(j,6) * v(i,6,k,e) &
1642 + dy(j,7) * v(i,7,k,e) &
1643 + dy(j,8) * v(i,8,k,e) &
1644 + dy(j,9) * v(i,9,k,e) &
1645 + dy(j,10) * v(i,10,k,e) &
1646 + dy(j,11) * v(i,11,k,e)
1647
1648 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
1649 + dy(j,2) * w(i,2,k,e) &
1650 + dy(j,3) * w(i,3,k,e) &
1651 + dy(j,4) * w(i,4,k,e) &
1652 + dy(j,5) * w(i,5,k,e) &
1653 + dy(j,6) * w(i,6,k,e) &
1654 + dy(j,7) * w(i,7,k,e) &
1655 + dy(j,8) * w(i,8,k,e) &
1656 + dy(j,9) * w(i,9,k,e) &
1657 + dy(j,10) * w(i,10,k,e) &
1658 + dy(j,11) * w(i,11,k,e)
1659 end do
1660 end do
1661 end do
1662
1663 do k = 1, lx
1664 do i = 1, lx*lx
1665 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1666 + dz(k,2) * u(i,1,2,e) &
1667 + dz(k,3) * u(i,1,3,e) &
1668 + dz(k,4) * u(i,1,4,e) &
1669 + dz(k,5) * u(i,1,5,e) &
1670 + dz(k,6) * u(i,1,6,e) &
1671 + dz(k,7) * u(i,1,7,e) &
1672 + dz(k,8) * u(i,1,8,e) &
1673 + dz(k,9) * u(i,1,9,e) &
1674 + dz(k,10) * u(i,1,10,e) &
1675 + dz(k,11) * u(i,1,11,e)
1676
1677 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
1678 + dz(k,2) * v(i,1,2,e) &
1679 + dz(k,3) * v(i,1,3,e) &
1680 + dz(k,4) * v(i,1,4,e) &
1681 + dz(k,5) * v(i,1,5,e) &
1682 + dz(k,6) * v(i,1,6,e) &
1683 + dz(k,7) * v(i,1,7,e) &
1684 + dz(k,8) * v(i,1,8,e) &
1685 + dz(k,9) * v(i,1,9,e) &
1686 + dz(k,10) * v(i,1,10,e) &
1687 + dz(k,11) * v(i,1,11,e)
1688
1689 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
1690 + dz(k,2) * w(i,1,2,e) &
1691 + dz(k,3) * w(i,1,3,e) &
1692 + dz(k,4) * w(i,1,4,e) &
1693 + dz(k,5) * w(i,1,5,e) &
1694 + dz(k,6) * w(i,1,6,e) &
1695 + dz(k,7) * w(i,1,7,e) &
1696 + dz(k,8) * w(i,1,8,e) &
1697 + dz(k,9) * w(i,1,9,e) &
1698 + dz(k,10) * w(i,1,10,e) &
1699 + dz(k,11) * w(i,1,11,e)
1700 end do
1701 end do
1702
1703 do i = 1, lx*lx*lx
1704
1705 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
1706 + wut(i,1,1) * dtdx(i,1,1,e)
1707 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
1708 + wut(i,1,1) * dtdy(i,1,1,e)
1709 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
1710 + wut(i,1,1) * dtdz(i,1,1,e)
1711
1712 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
1713 + wvt(i,1,1) * dtdx(i,1,1,e)
1714 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
1715 + wvt(i,1,1) * dtdy(i,1,1,e)
1716 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
1717 + wvt(i,1,1) * dtdz(i,1,1,e)
1718
1719 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
1720 + wwt(i,1,1) * dtdx(i,1,1,e)
1721 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
1722 + wwt(i,1,1) * dtdy(i,1,1,e)
1723 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
1724 + wwt(i,1,1) * dtdz(i,1,1,e)
1725
1726 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
1727 s11 = dj*(u1 + u1)
1728 s22 = dj*(v2 + v2)
1729 s33 = dj*(w3 + w3)
1730 s12 = dj*(u2 + v1)
1731 s13 = dj*(u3 + w1)
1732 s23 = dj*(v3 + w2)
1733
1734 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
1735 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
1736 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
1737
1738 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
1739 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
1740 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
1741
1742 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
1743 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
1744 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
1745 end do
1746
1747 do j = 1, lx*lx
1748 do i = 1, lx
1749 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
1750 + dxt(i,2) * wur(2,j,1) &
1751 + dxt(i,3) * wur(3,j,1) &
1752 + dxt(i,4) * wur(4,j,1) &
1753 + dxt(i,5) * wur(5,j,1) &
1754 + dxt(i,6) * wur(6,j,1) &
1755 + dxt(i,7) * wur(7,j,1) &
1756 + dxt(i,8) * wur(8,j,1) &
1757 + dxt(i,9) * wur(9,j,1) &
1758 + dxt(i,10) * wur(10,j,1) &
1759 + dxt(i,11) * wur(11,j,1)
1760
1761 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
1762 + dxt(i,2) * wvr(2,j,1) &
1763 + dxt(i,3) * wvr(3,j,1) &
1764 + dxt(i,4) * wvr(4,j,1) &
1765 + dxt(i,5) * wvr(5,j,1) &
1766 + dxt(i,6) * wvr(6,j,1) &
1767 + dxt(i,7) * wvr(7,j,1) &
1768 + dxt(i,8) * wvr(8,j,1) &
1769 + dxt(i,9) * wvr(9,j,1) &
1770 + dxt(i,10) * wvr(10,j,1) &
1771 + dxt(i,11) * wvr(11,j,1)
1772
1773 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
1774 + dxt(i,2) * wwr(2,j,1) &
1775 + dxt(i,3) * wwr(3,j,1) &
1776 + dxt(i,4) * wwr(4,j,1) &
1777 + dxt(i,5) * wwr(5,j,1) &
1778 + dxt(i,6) * wwr(6,j,1) &
1779 + dxt(i,7) * wwr(7,j,1) &
1780 + dxt(i,8) * wwr(8,j,1) &
1781 + dxt(i,9) * wwr(9,j,1) &
1782 + dxt(i,10) * wwr(10,j,1) &
1783 + dxt(i,11) * wwr(11,j,1)
1784 end do
1785 end do
1786
1787 do k = 1, lx
1788 do j = 1, lx
1789 do i = 1, lx
1790 au(i,j,k,e) = au(i,j,k,e) &
1791 + dyt(j,1) * wus(i,1,k) &
1792 + dyt(j,2) * wus(i,2,k) &
1793 + dyt(j,3) * wus(i,3,k) &
1794 + dyt(j,4) * wus(i,4,k) &
1795 + dyt(j,5) * wus(i,5,k) &
1796 + dyt(j,6) * wus(i,6,k) &
1797 + dyt(j,7) * wus(i,7,k) &
1798 + dyt(j,8) * wus(i,8,k) &
1799 + dyt(j,9) * wus(i,9,k) &
1800 + dyt(j,10) * wus(i,10,k) &
1801 + dyt(j,11) * wus(i,11,k)
1802
1803 av(i,j,k,e) = av(i,j,k,e) &
1804 + dyt(j,1) * wvs(i,1,k) &
1805 + dyt(j,2) * wvs(i,2,k) &
1806 + dyt(j,3) * wvs(i,3,k) &
1807 + dyt(j,4) * wvs(i,4,k) &
1808 + dyt(j,5) * wvs(i,5,k) &
1809 + dyt(j,6) * wvs(i,6,k) &
1810 + dyt(j,7) * wvs(i,7,k) &
1811 + dyt(j,8) * wvs(i,8,k) &
1812 + dyt(j,9) * wvs(i,9,k) &
1813 + dyt(j,10) * wvs(i,10,k) &
1814 + dyt(j,11) * wvs(i,11,k)
1815
1816 aw(i,j,k,e) = aw(i,j,k,e) &
1817 + dyt(j,1) * wws(i,1,k) &
1818 + dyt(j,2) * wws(i,2,k) &
1819 + dyt(j,3) * wws(i,3,k) &
1820 + dyt(j,4) * wws(i,4,k) &
1821 + dyt(j,5) * wws(i,5,k) &
1822 + dyt(j,6) * wws(i,6,k) &
1823 + dyt(j,7) * wws(i,7,k) &
1824 + dyt(j,8) * wws(i,8,k) &
1825 + dyt(j,9) * wws(i,9,k) &
1826 + dyt(j,10) * wws(i,10,k) &
1827 + dyt(j,11) * wws(i,11,k)
1828 end do
1829 end do
1830 end do
1831
1832 do k = 1, lx
1833 do i = 1, lx*lx
1834 au(i,1,k,e) = au(i,1,k,e) &
1835 + dzt(k,1) * wut(i,1,1) &
1836 + dzt(k,2) * wut(i,1,2) &
1837 + dzt(k,3) * wut(i,1,3) &
1838 + dzt(k,4) * wut(i,1,4) &
1839 + dzt(k,5) * wut(i,1,5) &
1840 + dzt(k,6) * wut(i,1,6) &
1841 + dzt(k,7) * wut(i,1,7) &
1842 + dzt(k,8) * wut(i,1,8) &
1843 + dzt(k,9) * wut(i,1,9) &
1844 + dzt(k,10) * wut(i,1,10) &
1845 + dzt(k,11) * wut(i,1,11)
1846
1847 av(i,1,k,e) = av(i,1,k,e) &
1848 + dzt(k,1) * wvt(i,1,1) &
1849 + dzt(k,2) * wvt(i,1,2) &
1850 + dzt(k,3) * wvt(i,1,3) &
1851 + dzt(k,4) * wvt(i,1,4) &
1852 + dzt(k,5) * wvt(i,1,5) &
1853 + dzt(k,6) * wvt(i,1,6) &
1854 + dzt(k,7) * wvt(i,1,7) &
1855 + dzt(k,8) * wvt(i,1,8) &
1856 + dzt(k,9) * wvt(i,1,9) &
1857 + dzt(k,10) * wvt(i,1,10) &
1858 + dzt(k,11) * wvt(i,1,11)
1859
1860 aw(i,1,k,e) = aw(i,1,k,e) &
1861 + dzt(k,1) * wwt(i,1,1) &
1862 + dzt(k,2) * wwt(i,1,2) &
1863 + dzt(k,3) * wwt(i,1,3) &
1864 + dzt(k,4) * wwt(i,1,4) &
1865 + dzt(k,5) * wwt(i,1,5) &
1866 + dzt(k,6) * wwt(i,1,6) &
1867 + dzt(k,7) * wwt(i,1,7) &
1868 + dzt(k,8) * wwt(i,1,8) &
1869 + dzt(k,9) * wwt(i,1,9) &
1870 + dzt(k,10) * wwt(i,1,10) &
1871 + dzt(k,11) * wwt(i,1,11)
1872 end do
1873 end do
1874
1875 end do
1876 !$omp end do
1877 end subroutine ax_helm_stress_lx11
1878
1879 subroutine ax_helm_stress_lx10(au, av, aw, u, v, w, Dx, Dy, Dz, &
1880 Dxt, Dyt, Dzt, &
1881 h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
1882 jacinv, weights3, n)
1883 integer, parameter :: lx = 10
1884 integer, intent(in) :: n
1885 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1886 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
1887 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
1888 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
1889 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
1890 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
1891 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1892 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
1893 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
1894 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
1895 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
1896 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
1897 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
1898 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
1899 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
1900 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
1901 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
1902 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
1903 real(kind=rp), intent(in) :: dx(lx, lx)
1904 real(kind=rp), intent(in) :: dy(lx, lx)
1905 real(kind=rp), intent(in) :: dz(lx, lx)
1906 real(kind=rp), intent(in) :: dxt(lx, lx)
1907 real(kind=rp), intent(in) :: dyt(lx, lx)
1908 real(kind=rp), intent(in) :: dzt(lx, lx)
1909 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
1910
1911 real(kind=rp) :: wur(lx, lx, lx)
1912 real(kind=rp) :: wus(lx, lx, lx)
1913 real(kind=rp) :: wut(lx, lx, lx)
1914 real(kind=rp) :: wvr(lx, lx, lx)
1915 real(kind=rp) :: wvs(lx, lx, lx)
1916 real(kind=rp) :: wvt(lx, lx, lx)
1917 real(kind=rp) :: wwr(lx, lx, lx)
1918 real(kind=rp) :: wws(lx, lx, lx)
1919 real(kind=rp) :: wwt(lx, lx, lx)
1920
1921 integer :: e, i, j, k, l
1922 real(kind=rp) :: dj
1923 real(kind=rp) :: s11, s22, s33, s12, s13, s23
1924 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
1925
1926 !$omp do
1927 do e = 1, n
1928
1929 do j = 1, lx * lx
1930 do i = 1, lx
1931 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1932 + dx(i,2) * u(2,j,1,e) &
1933 + dx(i,3) * u(3,j,1,e) &
1934 + dx(i,4) * u(4,j,1,e) &
1935 + dx(i,5) * u(5,j,1,e) &
1936 + dx(i,6) * u(6,j,1,e) &
1937 + dx(i,7) * u(7,j,1,e) &
1938 + dx(i,8) * u(8,j,1,e) &
1939 + dx(i,9) * u(9,j,1,e) &
1940 + dx(i,10) * u(10,j,1,e)
1941
1942 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
1943 + dx(i,2) * v(2,j,1,e) &
1944 + dx(i,3) * v(3,j,1,e) &
1945 + dx(i,4) * v(4,j,1,e) &
1946 + dx(i,5) * v(5,j,1,e) &
1947 + dx(i,6) * v(6,j,1,e) &
1948 + dx(i,7) * v(7,j,1,e) &
1949 + dx(i,8) * v(8,j,1,e) &
1950 + dx(i,9) * v(9,j,1,e) &
1951 + dx(i,10) * v(10,j,1,e)
1952
1953 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
1954 + dx(i,2) * w(2,j,1,e) &
1955 + dx(i,3) * w(3,j,1,e) &
1956 + dx(i,4) * w(4,j,1,e) &
1957 + dx(i,5) * w(5,j,1,e) &
1958 + dx(i,6) * w(6,j,1,e) &
1959 + dx(i,7) * w(7,j,1,e) &
1960 + dx(i,8) * w(8,j,1,e) &
1961 + dx(i,9) * w(9,j,1,e) &
1962 + dx(i,10) * w(10,j,1,e)
1963 end do
1964 end do
1965
1966 do k = 1, lx
1967 do j = 1, lx
1968 do i = 1, lx
1969 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1970 + dy(j,2) * u(i,2,k,e) &
1971 + dy(j,3) * u(i,3,k,e) &
1972 + dy(j,4) * u(i,4,k,e) &
1973 + dy(j,5) * u(i,5,k,e) &
1974 + dy(j,6) * u(i,6,k,e) &
1975 + dy(j,7) * u(i,7,k,e) &
1976 + dy(j,8) * u(i,8,k,e) &
1977 + dy(j,9) * u(i,9,k,e) &
1978 + dy(j,10) * u(i,10,k,e)
1979
1980 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
1981 + dy(j,2) * v(i,2,k,e) &
1982 + dy(j,3) * v(i,3,k,e) &
1983 + dy(j,4) * v(i,4,k,e) &
1984 + dy(j,5) * v(i,5,k,e) &
1985 + dy(j,6) * v(i,6,k,e) &
1986 + dy(j,7) * v(i,7,k,e) &
1987 + dy(j,8) * v(i,8,k,e) &
1988 + dy(j,9) * v(i,9,k,e) &
1989 + dy(j,10) * v(i,10,k,e)
1990
1991 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
1992 + dy(j,2) * w(i,2,k,e) &
1993 + dy(j,3) * w(i,3,k,e) &
1994 + dy(j,4) * w(i,4,k,e) &
1995 + dy(j,5) * w(i,5,k,e) &
1996 + dy(j,6) * w(i,6,k,e) &
1997 + dy(j,7) * w(i,7,k,e) &
1998 + dy(j,8) * w(i,8,k,e) &
1999 + dy(j,9) * w(i,9,k,e) &
2000 + dy(j,10) * w(i,10,k,e)
2001 end do
2002 end do
2003 end do
2004
2005 do k = 1, lx
2006 do i = 1, lx*lx
2007 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2008 + dz(k,2) * u(i,1,2,e) &
2009 + dz(k,3) * u(i,1,3,e) &
2010 + dz(k,4) * u(i,1,4,e) &
2011 + dz(k,5) * u(i,1,5,e) &
2012 + dz(k,6) * u(i,1,6,e) &
2013 + dz(k,7) * u(i,1,7,e) &
2014 + dz(k,8) * u(i,1,8,e) &
2015 + dz(k,9) * u(i,1,9,e) &
2016 + dz(k,10) * u(i,1,10,e)
2017
2018 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2019 + dz(k,2) * v(i,1,2,e) &
2020 + dz(k,3) * v(i,1,3,e) &
2021 + dz(k,4) * v(i,1,4,e) &
2022 + dz(k,5) * v(i,1,5,e) &
2023 + dz(k,6) * v(i,1,6,e) &
2024 + dz(k,7) * v(i,1,7,e) &
2025 + dz(k,8) * v(i,1,8,e) &
2026 + dz(k,9) * v(i,1,9,e) &
2027 + dz(k,10) * v(i,1,10,e)
2028
2029 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2030 + dz(k,2) * w(i,1,2,e) &
2031 + dz(k,3) * w(i,1,3,e) &
2032 + dz(k,4) * w(i,1,4,e) &
2033 + dz(k,5) * w(i,1,5,e) &
2034 + dz(k,6) * w(i,1,6,e) &
2035 + dz(k,7) * w(i,1,7,e) &
2036 + dz(k,8) * w(i,1,8,e) &
2037 + dz(k,9) * w(i,1,9,e) &
2038 + dz(k,10) * w(i,1,10,e)
2039 end do
2040 end do
2041
2042 do i = 1, lx*lx*lx
2043
2044 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
2045 + wut(i,1,1) * dtdx(i,1,1,e)
2046 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
2047 + wut(i,1,1) * dtdy(i,1,1,e)
2048 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
2049 + wut(i,1,1) * dtdz(i,1,1,e)
2050
2051 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
2052 + wvt(i,1,1) * dtdx(i,1,1,e)
2053 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
2054 + wvt(i,1,1) * dtdy(i,1,1,e)
2055 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
2056 + wvt(i,1,1) * dtdz(i,1,1,e)
2057
2058 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
2059 + wwt(i,1,1) * dtdx(i,1,1,e)
2060 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
2061 + wwt(i,1,1) * dtdy(i,1,1,e)
2062 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
2063 + wwt(i,1,1) * dtdz(i,1,1,e)
2064
2065 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
2066 s11 = dj*(u1 + u1)
2067 s22 = dj*(v2 + v2)
2068 s33 = dj*(w3 + w3)
2069 s12 = dj*(u2 + v1)
2070 s13 = dj*(u3 + w1)
2071 s23 = dj*(v3 + w2)
2072
2073 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
2074 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
2075 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
2076
2077 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
2078 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
2079 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
2080
2081 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
2082 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
2083 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
2084 end do
2085
2086 do j = 1, lx*lx
2087 do i = 1, lx
2088 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
2089 + dxt(i,2) * wur(2,j,1) &
2090 + dxt(i,3) * wur(3,j,1) &
2091 + dxt(i,4) * wur(4,j,1) &
2092 + dxt(i,5) * wur(5,j,1) &
2093 + dxt(i,6) * wur(6,j,1) &
2094 + dxt(i,7) * wur(7,j,1) &
2095 + dxt(i,8) * wur(8,j,1) &
2096 + dxt(i,9) * wur(9,j,1) &
2097 + dxt(i,10) * wur(10,j,1)
2098
2099 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
2100 + dxt(i,2) * wvr(2,j,1) &
2101 + dxt(i,3) * wvr(3,j,1) &
2102 + dxt(i,4) * wvr(4,j,1) &
2103 + dxt(i,5) * wvr(5,j,1) &
2104 + dxt(i,6) * wvr(6,j,1) &
2105 + dxt(i,7) * wvr(7,j,1) &
2106 + dxt(i,8) * wvr(8,j,1) &
2107 + dxt(i,9) * wvr(9,j,1) &
2108 + dxt(i,10) * wvr(10,j,1)
2109
2110 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
2111 + dxt(i,2) * wwr(2,j,1) &
2112 + dxt(i,3) * wwr(3,j,1) &
2113 + dxt(i,4) * wwr(4,j,1) &
2114 + dxt(i,5) * wwr(5,j,1) &
2115 + dxt(i,6) * wwr(6,j,1) &
2116 + dxt(i,7) * wwr(7,j,1) &
2117 + dxt(i,8) * wwr(8,j,1) &
2118 + dxt(i,9) * wwr(9,j,1) &
2119 + dxt(i,10) * wwr(10,j,1)
2120 end do
2121 end do
2122
2123 do k = 1, lx
2124 do j = 1, lx
2125 do i = 1, lx
2126 au(i,j,k,e) = au(i,j,k,e) &
2127 + dyt(j,1) * wus(i,1,k) &
2128 + dyt(j,2) * wus(i,2,k) &
2129 + dyt(j,3) * wus(i,3,k) &
2130 + dyt(j,4) * wus(i,4,k) &
2131 + dyt(j,5) * wus(i,5,k) &
2132 + dyt(j,6) * wus(i,6,k) &
2133 + dyt(j,7) * wus(i,7,k) &
2134 + dyt(j,8) * wus(i,8,k) &
2135 + dyt(j,9) * wus(i,9,k) &
2136 + dyt(j,10) * wus(i,10,k)
2137
2138 av(i,j,k,e) = av(i,j,k,e) &
2139 + dyt(j,1) * wvs(i,1,k) &
2140 + dyt(j,2) * wvs(i,2,k) &
2141 + dyt(j,3) * wvs(i,3,k) &
2142 + dyt(j,4) * wvs(i,4,k) &
2143 + dyt(j,5) * wvs(i,5,k) &
2144 + dyt(j,6) * wvs(i,6,k) &
2145 + dyt(j,7) * wvs(i,7,k) &
2146 + dyt(j,8) * wvs(i,8,k) &
2147 + dyt(j,9) * wvs(i,9,k) &
2148 + dyt(j,10) * wvs(i,10,k)
2149
2150 aw(i,j,k,e) = aw(i,j,k,e) &
2151 + dyt(j,1) * wws(i,1,k) &
2152 + dyt(j,2) * wws(i,2,k) &
2153 + dyt(j,3) * wws(i,3,k) &
2154 + dyt(j,4) * wws(i,4,k) &
2155 + dyt(j,5) * wws(i,5,k) &
2156 + dyt(j,6) * wws(i,6,k) &
2157 + dyt(j,7) * wws(i,7,k) &
2158 + dyt(j,8) * wws(i,8,k) &
2159 + dyt(j,9) * wws(i,9,k) &
2160 + dyt(j,10) * wws(i,10,k)
2161 end do
2162 end do
2163 end do
2164
2165 do k = 1, lx
2166 do i = 1, lx*lx
2167 au(i,1,k,e) = au(i,1,k,e) &
2168 + dzt(k,1) * wut(i,1,1) &
2169 + dzt(k,2) * wut(i,1,2) &
2170 + dzt(k,3) * wut(i,1,3) &
2171 + dzt(k,4) * wut(i,1,4) &
2172 + dzt(k,5) * wut(i,1,5) &
2173 + dzt(k,6) * wut(i,1,6) &
2174 + dzt(k,7) * wut(i,1,7) &
2175 + dzt(k,8) * wut(i,1,8) &
2176 + dzt(k,9) * wut(i,1,9) &
2177 + dzt(k,10) * wut(i,1,10)
2178
2179 av(i,1,k,e) = av(i,1,k,e) &
2180 + dzt(k,1) * wvt(i,1,1) &
2181 + dzt(k,2) * wvt(i,1,2) &
2182 + dzt(k,3) * wvt(i,1,3) &
2183 + dzt(k,4) * wvt(i,1,4) &
2184 + dzt(k,5) * wvt(i,1,5) &
2185 + dzt(k,6) * wvt(i,1,6) &
2186 + dzt(k,7) * wvt(i,1,7) &
2187 + dzt(k,8) * wvt(i,1,8) &
2188 + dzt(k,9) * wvt(i,1,9) &
2189 + dzt(k,10) * wvt(i,1,10)
2190
2191 aw(i,1,k,e) = aw(i,1,k,e) &
2192 + dzt(k,1) * wwt(i,1,1) &
2193 + dzt(k,2) * wwt(i,1,2) &
2194 + dzt(k,3) * wwt(i,1,3) &
2195 + dzt(k,4) * wwt(i,1,4) &
2196 + dzt(k,5) * wwt(i,1,5) &
2197 + dzt(k,6) * wwt(i,1,6) &
2198 + dzt(k,7) * wwt(i,1,7) &
2199 + dzt(k,8) * wwt(i,1,8) &
2200 + dzt(k,9) * wwt(i,1,9) &
2201 + dzt(k,10) * wwt(i,1,10)
2202 end do
2203 end do
2204
2205 end do
2206 !$omp end do
2207 end subroutine ax_helm_stress_lx10
2208
2209 subroutine ax_helm_stress_lx9(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
2210 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
2211 jacinv, weights3, n)
2212 integer, parameter :: lx = 9
2213 integer, intent(in) :: n
2214 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2215 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
2216 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
2217 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2218 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2219 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2220 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
2221 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
2222 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
2223 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
2224 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
2225 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
2226 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
2227 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
2228 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
2229 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
2230 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
2231 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
2232 real(kind=rp), intent(in) :: dx(lx, lx)
2233 real(kind=rp), intent(in) :: dy(lx, lx)
2234 real(kind=rp), intent(in) :: dz(lx, lx)
2235 real(kind=rp), intent(in) :: dxt(lx, lx)
2236 real(kind=rp), intent(in) :: dyt(lx, lx)
2237 real(kind=rp), intent(in) :: dzt(lx, lx)
2238 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
2239
2240 real(kind=rp) :: wur(lx, lx, lx)
2241 real(kind=rp) :: wus(lx, lx, lx)
2242 real(kind=rp) :: wut(lx, lx, lx)
2243 real(kind=rp) :: wvr(lx, lx, lx)
2244 real(kind=rp) :: wvs(lx, lx, lx)
2245 real(kind=rp) :: wvt(lx, lx, lx)
2246 real(kind=rp) :: wwr(lx, lx, lx)
2247 real(kind=rp) :: wws(lx, lx, lx)
2248 real(kind=rp) :: wwt(lx, lx, lx)
2249
2250 integer :: e, i, j, k, l
2251 real(kind=rp) :: dj
2252 real(kind=rp) :: s11, s22, s33, s12, s13, s23
2253 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
2254
2255 !$omp do
2256 do e = 1, n
2257
2258 do j = 1, lx * lx
2259 do i = 1, lx
2260 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
2261 + dx(i,2) * u(2,j,1,e) &
2262 + dx(i,3) * u(3,j,1,e) &
2263 + dx(i,4) * u(4,j,1,e) &
2264 + dx(i,5) * u(5,j,1,e) &
2265 + dx(i,6) * u(6,j,1,e) &
2266 + dx(i,7) * u(7,j,1,e) &
2267 + dx(i,8) * u(8,j,1,e) &
2268 + dx(i,9) * u(9,j,1,e)
2269
2270 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
2271 + dx(i,2) * v(2,j,1,e) &
2272 + dx(i,3) * v(3,j,1,e) &
2273 + dx(i,4) * v(4,j,1,e) &
2274 + dx(i,5) * v(5,j,1,e) &
2275 + dx(i,6) * v(6,j,1,e) &
2276 + dx(i,7) * v(7,j,1,e) &
2277 + dx(i,8) * v(8,j,1,e) &
2278 + dx(i,9) * v(9,j,1,e)
2279
2280 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
2281 + dx(i,2) * w(2,j,1,e) &
2282 + dx(i,3) * w(3,j,1,e) &
2283 + dx(i,4) * w(4,j,1,e) &
2284 + dx(i,5) * w(5,j,1,e) &
2285 + dx(i,6) * w(6,j,1,e) &
2286 + dx(i,7) * w(7,j,1,e) &
2287 + dx(i,8) * w(8,j,1,e) &
2288 + dx(i,9) * w(9,j,1,e)
2289 end do
2290 end do
2291
2292 do k = 1, lx
2293 do j = 1, lx
2294 do i = 1, lx
2295 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
2296 + dy(j,2) * u(i,2,k,e) &
2297 + dy(j,3) * u(i,3,k,e) &
2298 + dy(j,4) * u(i,4,k,e) &
2299 + dy(j,5) * u(i,5,k,e) &
2300 + dy(j,6) * u(i,6,k,e) &
2301 + dy(j,7) * u(i,7,k,e) &
2302 + dy(j,8) * u(i,8,k,e) &
2303 + dy(j,9) * u(i,9,k,e)
2304
2305 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
2306 + dy(j,2) * v(i,2,k,e) &
2307 + dy(j,3) * v(i,3,k,e) &
2308 + dy(j,4) * v(i,4,k,e) &
2309 + dy(j,5) * v(i,5,k,e) &
2310 + dy(j,6) * v(i,6,k,e) &
2311 + dy(j,7) * v(i,7,k,e) &
2312 + dy(j,8) * v(i,8,k,e) &
2313 + dy(j,9) * v(i,9,k,e)
2314
2315 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
2316 + dy(j,2) * w(i,2,k,e) &
2317 + dy(j,3) * w(i,3,k,e) &
2318 + dy(j,4) * w(i,4,k,e) &
2319 + dy(j,5) * w(i,5,k,e) &
2320 + dy(j,6) * w(i,6,k,e) &
2321 + dy(j,7) * w(i,7,k,e) &
2322 + dy(j,8) * w(i,8,k,e) &
2323 + dy(j,9) * w(i,9,k,e)
2324 end do
2325 end do
2326 end do
2327
2328 do k = 1, lx
2329 do i = 1, lx*lx
2330 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2331 + dz(k,2) * u(i,1,2,e) &
2332 + dz(k,3) * u(i,1,3,e) &
2333 + dz(k,4) * u(i,1,4,e) &
2334 + dz(k,5) * u(i,1,5,e) &
2335 + dz(k,6) * u(i,1,6,e) &
2336 + dz(k,7) * u(i,1,7,e) &
2337 + dz(k,8) * u(i,1,8,e) &
2338 + dz(k,9) * u(i,1,9,e)
2339
2340 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2341 + dz(k,2) * v(i,1,2,e) &
2342 + dz(k,3) * v(i,1,3,e) &
2343 + dz(k,4) * v(i,1,4,e) &
2344 + dz(k,5) * v(i,1,5,e) &
2345 + dz(k,6) * v(i,1,6,e) &
2346 + dz(k,7) * v(i,1,7,e) &
2347 + dz(k,8) * v(i,1,8,e) &
2348 + dz(k,9) * v(i,1,9,e)
2349
2350 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2351 + dz(k,2) * w(i,1,2,e) &
2352 + dz(k,3) * w(i,1,3,e) &
2353 + dz(k,4) * w(i,1,4,e) &
2354 + dz(k,5) * w(i,1,5,e) &
2355 + dz(k,6) * w(i,1,6,e) &
2356 + dz(k,7) * w(i,1,7,e) &
2357 + dz(k,8) * w(i,1,8,e) &
2358 + dz(k,9) * w(i,1,9,e)
2359 end do
2360 end do
2361
2362 do i = 1, lx*lx*lx
2363
2364 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
2365 + wut(i,1,1) * dtdx(i,1,1,e)
2366 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
2367 + wut(i,1,1) * dtdy(i,1,1,e)
2368 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
2369 + wut(i,1,1) * dtdz(i,1,1,e)
2370
2371 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
2372 + wvt(i,1,1) * dtdx(i,1,1,e)
2373 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
2374 + wvt(i,1,1) * dtdy(i,1,1,e)
2375 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
2376 + wvt(i,1,1) * dtdz(i,1,1,e)
2377
2378 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
2379 + wwt(i,1,1) * dtdx(i,1,1,e)
2380 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
2381 + wwt(i,1,1) * dtdy(i,1,1,e)
2382 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
2383 + wwt(i,1,1) * dtdz(i,1,1,e)
2384
2385 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
2386 s11 = dj*(u1 + u1)
2387 s22 = dj*(v2 + v2)
2388 s33 = dj*(w3 + w3)
2389 s12 = dj*(u2 + v1)
2390 s13 = dj*(u3 + w1)
2391 s23 = dj*(v3 + w2)
2392
2393 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
2394 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
2395 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
2396
2397 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
2398 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
2399 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
2400
2401 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
2402 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
2403 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
2404 end do
2405
2406 do j = 1, lx*lx
2407 do i = 1, lx
2408 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
2409 + dxt(i,2) * wur(2,j,1) &
2410 + dxt(i,3) * wur(3,j,1) &
2411 + dxt(i,4) * wur(4,j,1) &
2412 + dxt(i,5) * wur(5,j,1) &
2413 + dxt(i,6) * wur(6,j,1) &
2414 + dxt(i,7) * wur(7,j,1) &
2415 + dxt(i,8) * wur(8,j,1) &
2416 + dxt(i,9) * wur(9,j,1)
2417
2418 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
2419 + dxt(i,2) * wvr(2,j,1) &
2420 + dxt(i,3) * wvr(3,j,1) &
2421 + dxt(i,4) * wvr(4,j,1) &
2422 + dxt(i,5) * wvr(5,j,1) &
2423 + dxt(i,6) * wvr(6,j,1) &
2424 + dxt(i,7) * wvr(7,j,1) &
2425 + dxt(i,8) * wvr(8,j,1) &
2426 + dxt(i,9) * wvr(9,j,1)
2427
2428 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
2429 + dxt(i,2) * wwr(2,j,1) &
2430 + dxt(i,3) * wwr(3,j,1) &
2431 + dxt(i,4) * wwr(4,j,1) &
2432 + dxt(i,5) * wwr(5,j,1) &
2433 + dxt(i,6) * wwr(6,j,1) &
2434 + dxt(i,7) * wwr(7,j,1) &
2435 + dxt(i,8) * wwr(8,j,1) &
2436 + dxt(i,9) * wwr(9,j,1)
2437 end do
2438 end do
2439
2440 do k = 1, lx
2441 do j = 1, lx
2442 do i = 1, lx
2443 au(i,j,k,e) = au(i,j,k,e) &
2444 + dyt(j,1) * wus(i,1,k) &
2445 + dyt(j,2) * wus(i,2,k) &
2446 + dyt(j,3) * wus(i,3,k) &
2447 + dyt(j,4) * wus(i,4,k) &
2448 + dyt(j,5) * wus(i,5,k) &
2449 + dyt(j,6) * wus(i,6,k) &
2450 + dyt(j,7) * wus(i,7,k) &
2451 + dyt(j,8) * wus(i,8,k) &
2452 + dyt(j,9) * wus(i,9,k)
2453
2454 av(i,j,k,e) = av(i,j,k,e) &
2455 + dyt(j,1) * wvs(i,1,k) &
2456 + dyt(j,2) * wvs(i,2,k) &
2457 + dyt(j,3) * wvs(i,3,k) &
2458 + dyt(j,4) * wvs(i,4,k) &
2459 + dyt(j,5) * wvs(i,5,k) &
2460 + dyt(j,6) * wvs(i,6,k) &
2461 + dyt(j,7) * wvs(i,7,k) &
2462 + dyt(j,8) * wvs(i,8,k) &
2463 + dyt(j,9) * wvs(i,9,k)
2464
2465 aw(i,j,k,e) = aw(i,j,k,e) &
2466 + dyt(j,1) * wws(i,1,k) &
2467 + dyt(j,2) * wws(i,2,k) &
2468 + dyt(j,3) * wws(i,3,k) &
2469 + dyt(j,4) * wws(i,4,k) &
2470 + dyt(j,5) * wws(i,5,k) &
2471 + dyt(j,6) * wws(i,6,k) &
2472 + dyt(j,7) * wws(i,7,k) &
2473 + dyt(j,8) * wws(i,8,k) &
2474 + dyt(j,9) * wws(i,9,k)
2475 end do
2476 end do
2477 end do
2478
2479 do k = 1, lx
2480 do i = 1, lx*lx
2481 au(i,1,k,e) = au(i,1,k,e) &
2482 + dzt(k,1) * wut(i,1,1) &
2483 + dzt(k,2) * wut(i,1,2) &
2484 + dzt(k,3) * wut(i,1,3) &
2485 + dzt(k,4) * wut(i,1,4) &
2486 + dzt(k,5) * wut(i,1,5) &
2487 + dzt(k,6) * wut(i,1,6) &
2488 + dzt(k,7) * wut(i,1,7) &
2489 + dzt(k,8) * wut(i,1,8) &
2490 + dzt(k,9) * wut(i,1,9)
2491
2492 av(i,1,k,e) = av(i,1,k,e) &
2493 + dzt(k,1) * wvt(i,1,1) &
2494 + dzt(k,2) * wvt(i,1,2) &
2495 + dzt(k,3) * wvt(i,1,3) &
2496 + dzt(k,4) * wvt(i,1,4) &
2497 + dzt(k,5) * wvt(i,1,5) &
2498 + dzt(k,6) * wvt(i,1,6) &
2499 + dzt(k,7) * wvt(i,1,7) &
2500 + dzt(k,8) * wvt(i,1,8) &
2501 + dzt(k,9) * wvt(i,1,9)
2502
2503 aw(i,1,k,e) = aw(i,1,k,e) &
2504 + dzt(k,1) * wwt(i,1,1) &
2505 + dzt(k,2) * wwt(i,1,2) &
2506 + dzt(k,3) * wwt(i,1,3) &
2507 + dzt(k,4) * wwt(i,1,4) &
2508 + dzt(k,5) * wwt(i,1,5) &
2509 + dzt(k,6) * wwt(i,1,6) &
2510 + dzt(k,7) * wwt(i,1,7) &
2511 + dzt(k,8) * wwt(i,1,8) &
2512 + dzt(k,9) * wwt(i,1,9)
2513 end do
2514 end do
2515
2516 end do
2517 !$omp end do
2518 end subroutine ax_helm_stress_lx9
2519
2520 subroutine ax_helm_stress_lx8(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
2521 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
2522 jacinv, weights3, n)
2523 integer, parameter :: lx = 8
2524 integer, intent(in) :: n
2525 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2526 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
2527 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
2528 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2529 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2530 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2531 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
2532 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
2533 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
2534 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
2535 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
2536 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
2537 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
2538 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
2539 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
2540 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
2541 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
2542 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
2543 real(kind=rp), intent(in) :: dx(lx, lx)
2544 real(kind=rp), intent(in) :: dy(lx, lx)
2545 real(kind=rp), intent(in) :: dz(lx, lx)
2546 real(kind=rp), intent(in) :: dxt(lx, lx)
2547 real(kind=rp), intent(in) :: dyt(lx, lx)
2548 real(kind=rp), intent(in) :: dzt(lx, lx)
2549 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
2550
2551 real(kind=rp) :: wur(lx, lx, lx)
2552 real(kind=rp) :: wus(lx, lx, lx)
2553 real(kind=rp) :: wut(lx, lx, lx)
2554 real(kind=rp) :: wvr(lx, lx, lx)
2555 real(kind=rp) :: wvs(lx, lx, lx)
2556 real(kind=rp) :: wvt(lx, lx, lx)
2557 real(kind=rp) :: wwr(lx, lx, lx)
2558 real(kind=rp) :: wws(lx, lx, lx)
2559 real(kind=rp) :: wwt(lx, lx, lx)
2560
2561 integer :: e, i, j, k, l
2562 real(kind=rp) :: dj
2563 real(kind=rp) :: s11, s22, s33, s12, s13, s23
2564 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
2565
2566 !$omp do
2567 do e = 1, n
2568
2569 do j = 1, lx * lx
2570 do i = 1, lx
2571 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
2572 + dx(i,2) * u(2,j,1,e) &
2573 + dx(i,3) * u(3,j,1,e) &
2574 + dx(i,4) * u(4,j,1,e) &
2575 + dx(i,5) * u(5,j,1,e) &
2576 + dx(i,6) * u(6,j,1,e) &
2577 + dx(i,7) * u(7,j,1,e) &
2578 + dx(i,8) * u(8,j,1,e)
2579
2580 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
2581 + dx(i,2) * v(2,j,1,e) &
2582 + dx(i,3) * v(3,j,1,e) &
2583 + dx(i,4) * v(4,j,1,e) &
2584 + dx(i,5) * v(5,j,1,e) &
2585 + dx(i,6) * v(6,j,1,e) &
2586 + dx(i,7) * v(7,j,1,e) &
2587 + dx(i,8) * v(8,j,1,e)
2588
2589 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
2590 + dx(i,2) * w(2,j,1,e) &
2591 + dx(i,3) * w(3,j,1,e) &
2592 + dx(i,4) * w(4,j,1,e) &
2593 + dx(i,5) * w(5,j,1,e) &
2594 + dx(i,6) * w(6,j,1,e) &
2595 + dx(i,7) * w(7,j,1,e) &
2596 + dx(i,8) * w(8,j,1,e)
2597 end do
2598 end do
2599
2600 do k = 1, lx
2601 do j = 1, lx
2602 do i = 1, lx
2603 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
2604 + dy(j,2) * u(i,2,k,e) &
2605 + dy(j,3) * u(i,3,k,e) &
2606 + dy(j,4) * u(i,4,k,e) &
2607 + dy(j,5) * u(i,5,k,e) &
2608 + dy(j,6) * u(i,6,k,e) &
2609 + dy(j,7) * u(i,7,k,e) &
2610 + dy(j,8) * u(i,8,k,e)
2611
2612 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
2613 + dy(j,2) * v(i,2,k,e) &
2614 + dy(j,3) * v(i,3,k,e) &
2615 + dy(j,4) * v(i,4,k,e) &
2616 + dy(j,5) * v(i,5,k,e) &
2617 + dy(j,6) * v(i,6,k,e) &
2618 + dy(j,7) * v(i,7,k,e) &
2619 + dy(j,8) * v(i,8,k,e)
2620
2621 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
2622 + dy(j,2) * w(i,2,k,e) &
2623 + dy(j,3) * w(i,3,k,e) &
2624 + dy(j,4) * w(i,4,k,e) &
2625 + dy(j,5) * w(i,5,k,e) &
2626 + dy(j,6) * w(i,6,k,e) &
2627 + dy(j,7) * w(i,7,k,e) &
2628 + dy(j,8) * w(i,8,k,e)
2629 end do
2630 end do
2631 end do
2632
2633 do k = 1, lx
2634 do i = 1, lx*lx
2635 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2636 + dz(k,2) * u(i,1,2,e) &
2637 + dz(k,3) * u(i,1,3,e) &
2638 + dz(k,4) * u(i,1,4,e) &
2639 + dz(k,5) * u(i,1,5,e) &
2640 + dz(k,6) * u(i,1,6,e) &
2641 + dz(k,7) * u(i,1,7,e) &
2642 + dz(k,8) * u(i,1,8,e)
2643
2644 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2645 + dz(k,2) * v(i,1,2,e) &
2646 + dz(k,3) * v(i,1,3,e) &
2647 + dz(k,4) * v(i,1,4,e) &
2648 + dz(k,5) * v(i,1,5,e) &
2649 + dz(k,6) * v(i,1,6,e) &
2650 + dz(k,7) * v(i,1,7,e) &
2651 + dz(k,8) * v(i,1,8,e)
2652
2653 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2654 + dz(k,2) * w(i,1,2,e) &
2655 + dz(k,3) * w(i,1,3,e) &
2656 + dz(k,4) * w(i,1,4,e) &
2657 + dz(k,5) * w(i,1,5,e) &
2658 + dz(k,6) * w(i,1,6,e) &
2659 + dz(k,7) * w(i,1,7,e) &
2660 + dz(k,8) * w(i,1,8,e)
2661 end do
2662 end do
2663
2664 do i = 1, lx*lx*lx
2665
2666 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
2667 + wut(i,1,1) * dtdx(i,1,1,e)
2668 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
2669 + wut(i,1,1) * dtdy(i,1,1,e)
2670 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
2671 + wut(i,1,1) * dtdz(i,1,1,e)
2672
2673 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
2674 + wvt(i,1,1) * dtdx(i,1,1,e)
2675 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
2676 + wvt(i,1,1) * dtdy(i,1,1,e)
2677 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
2678 + wvt(i,1,1) * dtdz(i,1,1,e)
2679
2680 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
2681 + wwt(i,1,1) * dtdx(i,1,1,e)
2682 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
2683 + wwt(i,1,1) * dtdy(i,1,1,e)
2684 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
2685 + wwt(i,1,1) * dtdz(i,1,1,e)
2686
2687 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
2688 s11 = dj*(u1 + u1)
2689 s22 = dj*(v2 + v2)
2690 s33 = dj*(w3 + w3)
2691 s12 = dj*(u2 + v1)
2692 s13 = dj*(u3 + w1)
2693 s23 = dj*(v3 + w2)
2694
2695 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
2696 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
2697 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
2698
2699 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
2700 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
2701 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
2702
2703 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
2704 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
2705 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
2706 end do
2707
2708 do j = 1, lx*lx
2709 do i = 1, lx
2710 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
2711 + dxt(i,2) * wur(2,j,1) &
2712 + dxt(i,3) * wur(3,j,1) &
2713 + dxt(i,4) * wur(4,j,1) &
2714 + dxt(i,5) * wur(5,j,1) &
2715 + dxt(i,6) * wur(6,j,1) &
2716 + dxt(i,7) * wur(7,j,1) &
2717 + dxt(i,8) * wur(8,j,1)
2718
2719 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
2720 + dxt(i,2) * wvr(2,j,1) &
2721 + dxt(i,3) * wvr(3,j,1) &
2722 + dxt(i,4) * wvr(4,j,1) &
2723 + dxt(i,5) * wvr(5,j,1) &
2724 + dxt(i,6) * wvr(6,j,1) &
2725 + dxt(i,7) * wvr(7,j,1) &
2726 + dxt(i,8) * wvr(8,j,1)
2727
2728 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
2729 + dxt(i,2) * wwr(2,j,1) &
2730 + dxt(i,3) * wwr(3,j,1) &
2731 + dxt(i,4) * wwr(4,j,1) &
2732 + dxt(i,5) * wwr(5,j,1) &
2733 + dxt(i,6) * wwr(6,j,1) &
2734 + dxt(i,7) * wwr(7,j,1) &
2735 + dxt(i,8) * wwr(8,j,1)
2736
2737 end do
2738 end do
2739
2740 do k = 1, lx
2741 do j = 1, lx
2742 do i = 1, lx
2743 au(i,j,k,e) = au(i,j,k,e) &
2744 + dyt(j,1) * wus(i,1,k) &
2745 + dyt(j,2) * wus(i,2,k) &
2746 + dyt(j,3) * wus(i,3,k) &
2747 + dyt(j,4) * wus(i,4,k) &
2748 + dyt(j,5) * wus(i,5,k) &
2749 + dyt(j,6) * wus(i,6,k) &
2750 + dyt(j,7) * wus(i,7,k) &
2751 + dyt(j,8) * wus(i,8,k)
2752
2753 av(i,j,k,e) = av(i,j,k,e) &
2754 + dyt(j,1) * wvs(i,1,k) &
2755 + dyt(j,2) * wvs(i,2,k) &
2756 + dyt(j,3) * wvs(i,3,k) &
2757 + dyt(j,4) * wvs(i,4,k) &
2758 + dyt(j,5) * wvs(i,5,k) &
2759 + dyt(j,6) * wvs(i,6,k) &
2760 + dyt(j,7) * wvs(i,7,k) &
2761 + dyt(j,8) * wvs(i,8,k)
2762
2763 aw(i,j,k,e) = aw(i,j,k,e) &
2764 + dyt(j,1) * wws(i,1,k) &
2765 + dyt(j,2) * wws(i,2,k) &
2766 + dyt(j,3) * wws(i,3,k) &
2767 + dyt(j,4) * wws(i,4,k) &
2768 + dyt(j,5) * wws(i,5,k) &
2769 + dyt(j,6) * wws(i,6,k) &
2770 + dyt(j,7) * wws(i,7,k) &
2771 + dyt(j,8) * wws(i,8,k)
2772 end do
2773 end do
2774 end do
2775
2776 do k = 1, lx
2777 do i = 1, lx*lx
2778 au(i,1,k,e) = au(i,1,k,e) &
2779 + dzt(k,1) * wut(i,1,1) &
2780 + dzt(k,2) * wut(i,1,2) &
2781 + dzt(k,3) * wut(i,1,3) &
2782 + dzt(k,4) * wut(i,1,4) &
2783 + dzt(k,5) * wut(i,1,5) &
2784 + dzt(k,6) * wut(i,1,6) &
2785 + dzt(k,7) * wut(i,1,7) &
2786 + dzt(k,8) * wut(i,1,8)
2787
2788 av(i,1,k,e) = av(i,1,k,e) &
2789 + dzt(k,1) * wvt(i,1,1) &
2790 + dzt(k,2) * wvt(i,1,2) &
2791 + dzt(k,3) * wvt(i,1,3) &
2792 + dzt(k,4) * wvt(i,1,4) &
2793 + dzt(k,5) * wvt(i,1,5) &
2794 + dzt(k,6) * wvt(i,1,6) &
2795 + dzt(k,7) * wvt(i,1,7) &
2796 + dzt(k,8) * wvt(i,1,8)
2797
2798 aw(i,1,k,e) = aw(i,1,k,e) &
2799 + dzt(k,1) * wwt(i,1,1) &
2800 + dzt(k,2) * wwt(i,1,2) &
2801 + dzt(k,3) * wwt(i,1,3) &
2802 + dzt(k,4) * wwt(i,1,4) &
2803 + dzt(k,5) * wwt(i,1,5) &
2804 + dzt(k,6) * wwt(i,1,6) &
2805 + dzt(k,7) * wwt(i,1,7) &
2806 + dzt(k,8) * wwt(i,1,8)
2807 end do
2808 end do
2809
2810 end do
2811 !$omp end do
2812 end subroutine ax_helm_stress_lx8
2813
2814 subroutine ax_helm_stress_lx7(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
2815 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
2816 jacinv, weights3, n)
2817 integer, parameter :: lx = 7
2818 integer, intent(in) :: n
2819 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2820 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
2821 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
2822 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2823 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2824 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2825 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
2826 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
2827 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
2828 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
2829 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
2830 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
2831 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
2832 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
2833 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
2834 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
2835 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
2836 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
2837 real(kind=rp), intent(in) :: dx(lx, lx)
2838 real(kind=rp), intent(in) :: dy(lx, lx)
2839 real(kind=rp), intent(in) :: dz(lx, lx)
2840 real(kind=rp), intent(in) :: dxt(lx, lx)
2841 real(kind=rp), intent(in) :: dyt(lx, lx)
2842 real(kind=rp), intent(in) :: dzt(lx, lx)
2843 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
2844
2845 real(kind=rp) :: wur(lx, lx, lx)
2846 real(kind=rp) :: wus(lx, lx, lx)
2847 real(kind=rp) :: wut(lx, lx, lx)
2848 real(kind=rp) :: wvr(lx, lx, lx)
2849 real(kind=rp) :: wvs(lx, lx, lx)
2850 real(kind=rp) :: wvt(lx, lx, lx)
2851 real(kind=rp) :: wwr(lx, lx, lx)
2852 real(kind=rp) :: wws(lx, lx, lx)
2853 real(kind=rp) :: wwt(lx, lx, lx)
2854
2855 integer :: e, i, j, k, l
2856 real(kind=rp) :: dj
2857 real(kind=rp) :: s11, s22, s33, s12, s13, s23
2858 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
2859
2860 !$omp do
2861 do e = 1, n
2862
2863 do j = 1, lx * lx
2864 do i = 1, lx
2865 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
2866 + dx(i,2) * u(2,j,1,e) &
2867 + dx(i,3) * u(3,j,1,e) &
2868 + dx(i,4) * u(4,j,1,e) &
2869 + dx(i,5) * u(5,j,1,e) &
2870 + dx(i,6) * u(6,j,1,e) &
2871 + dx(i,7) * u(7,j,1,e)
2872
2873 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
2874 + dx(i,2) * v(2,j,1,e) &
2875 + dx(i,3) * v(3,j,1,e) &
2876 + dx(i,4) * v(4,j,1,e) &
2877 + dx(i,5) * v(5,j,1,e) &
2878 + dx(i,6) * v(6,j,1,e) &
2879 + dx(i,7) * v(7,j,1,e)
2880
2881 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
2882 + dx(i,2) * w(2,j,1,e) &
2883 + dx(i,3) * w(3,j,1,e) &
2884 + dx(i,4) * w(4,j,1,e) &
2885 + dx(i,5) * w(5,j,1,e) &
2886 + dx(i,6) * w(6,j,1,e) &
2887 + dx(i,7) * w(7,j,1,e)
2888 end do
2889 end do
2890
2891 do k = 1, lx
2892 do j = 1, lx
2893 do i = 1, lx
2894 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
2895 + dy(j,2) * u(i,2,k,e) &
2896 + dy(j,3) * u(i,3,k,e) &
2897 + dy(j,4) * u(i,4,k,e) &
2898 + dy(j,5) * u(i,5,k,e) &
2899 + dy(j,6) * u(i,6,k,e) &
2900 + dy(j,7) * u(i,7,k,e)
2901
2902 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
2903 + dy(j,2) * v(i,2,k,e) &
2904 + dy(j,3) * v(i,3,k,e) &
2905 + dy(j,4) * v(i,4,k,e) &
2906 + dy(j,5) * v(i,5,k,e) &
2907 + dy(j,6) * v(i,6,k,e) &
2908 + dy(j,7) * v(i,7,k,e)
2909
2910 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
2911 + dy(j,2) * w(i,2,k,e) &
2912 + dy(j,3) * w(i,3,k,e) &
2913 + dy(j,4) * w(i,4,k,e) &
2914 + dy(j,5) * w(i,5,k,e) &
2915 + dy(j,6) * w(i,6,k,e) &
2916 + dy(j,7) * w(i,7,k,e)
2917
2918 end do
2919 end do
2920 end do
2921
2922 do k = 1, lx
2923 do i = 1, lx*lx
2924 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2925 + dz(k,2) * u(i,1,2,e) &
2926 + dz(k,3) * u(i,1,3,e) &
2927 + dz(k,4) * u(i,1,4,e) &
2928 + dz(k,5) * u(i,1,5,e) &
2929 + dz(k,6) * u(i,1,6,e) &
2930 + dz(k,7) * u(i,1,7,e)
2931
2932 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2933 + dz(k,2) * v(i,1,2,e) &
2934 + dz(k,3) * v(i,1,3,e) &
2935 + dz(k,4) * v(i,1,4,e) &
2936 + dz(k,5) * v(i,1,5,e) &
2937 + dz(k,6) * v(i,1,6,e) &
2938 + dz(k,7) * v(i,1,7,e)
2939
2940 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2941 + dz(k,2) * w(i,1,2,e) &
2942 + dz(k,3) * w(i,1,3,e) &
2943 + dz(k,4) * w(i,1,4,e) &
2944 + dz(k,5) * w(i,1,5,e) &
2945 + dz(k,6) * w(i,1,6,e) &
2946 + dz(k,7) * w(i,1,7,e)
2947 end do
2948 end do
2949
2950 do i = 1, lx*lx*lx
2951
2952 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
2953 + wut(i,1,1) * dtdx(i,1,1,e)
2954 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
2955 + wut(i,1,1) * dtdy(i,1,1,e)
2956 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
2957 + wut(i,1,1) * dtdz(i,1,1,e)
2958
2959 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
2960 + wvt(i,1,1) * dtdx(i,1,1,e)
2961 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
2962 + wvt(i,1,1) * dtdy(i,1,1,e)
2963 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
2964 + wvt(i,1,1) * dtdz(i,1,1,e)
2965
2966 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
2967 + wwt(i,1,1) * dtdx(i,1,1,e)
2968 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
2969 + wwt(i,1,1) * dtdy(i,1,1,e)
2970 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
2971 + wwt(i,1,1) * dtdz(i,1,1,e)
2972
2973 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
2974 s11 = dj*(u1 + u1)
2975 s22 = dj*(v2 + v2)
2976 s33 = dj*(w3 + w3)
2977 s12 = dj*(u2 + v1)
2978 s13 = dj*(u3 + w1)
2979 s23 = dj*(v3 + w2)
2980
2981 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
2982 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
2983 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
2984
2985 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
2986 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
2987 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
2988
2989 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
2990 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
2991 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
2992 end do
2993
2994 do j = 1, lx*lx
2995 do i = 1, lx
2996 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
2997 + dxt(i,2) * wur(2,j,1) &
2998 + dxt(i,3) * wur(3,j,1) &
2999 + dxt(i,4) * wur(4,j,1) &
3000 + dxt(i,5) * wur(5,j,1) &
3001 + dxt(i,6) * wur(6,j,1) &
3002 + dxt(i,7) * wur(7,j,1)
3003
3004 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
3005 + dxt(i,2) * wvr(2,j,1) &
3006 + dxt(i,3) * wvr(3,j,1) &
3007 + dxt(i,4) * wvr(4,j,1) &
3008 + dxt(i,5) * wvr(5,j,1) &
3009 + dxt(i,6) * wvr(6,j,1) &
3010 + dxt(i,7) * wvr(7,j,1)
3011
3012 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
3013 + dxt(i,2) * wwr(2,j,1) &
3014 + dxt(i,3) * wwr(3,j,1) &
3015 + dxt(i,4) * wwr(4,j,1) &
3016 + dxt(i,5) * wwr(5,j,1) &
3017 + dxt(i,6) * wwr(6,j,1) &
3018 + dxt(i,7) * wwr(7,j,1)
3019
3020 end do
3021 end do
3022
3023 do k = 1, lx
3024 do j = 1, lx
3025 do i = 1, lx
3026 au(i,j,k,e) = au(i,j,k,e) &
3027 + dyt(j,1) * wus(i,1,k) &
3028 + dyt(j,2) * wus(i,2,k) &
3029 + dyt(j,3) * wus(i,3,k) &
3030 + dyt(j,4) * wus(i,4,k) &
3031 + dyt(j,5) * wus(i,5,k) &
3032 + dyt(j,6) * wus(i,6,k) &
3033 + dyt(j,7) * wus(i,7,k)
3034
3035 av(i,j,k,e) = av(i,j,k,e) &
3036 + dyt(j,1) * wvs(i,1,k) &
3037 + dyt(j,2) * wvs(i,2,k) &
3038 + dyt(j,3) * wvs(i,3,k) &
3039 + dyt(j,4) * wvs(i,4,k) &
3040 + dyt(j,5) * wvs(i,5,k) &
3041 + dyt(j,6) * wvs(i,6,k) &
3042 + dyt(j,7) * wvs(i,7,k)
3043
3044 aw(i,j,k,e) = aw(i,j,k,e) &
3045 + dyt(j,1) * wws(i,1,k) &
3046 + dyt(j,2) * wws(i,2,k) &
3047 + dyt(j,3) * wws(i,3,k) &
3048 + dyt(j,4) * wws(i,4,k) &
3049 + dyt(j,5) * wws(i,5,k) &
3050 + dyt(j,6) * wws(i,6,k) &
3051 + dyt(j,7) * wws(i,7,k)
3052 end do
3053 end do
3054 end do
3055
3056 do k = 1, lx
3057 do i = 1, lx*lx
3058 au(i,1,k,e) = au(i,1,k,e) &
3059 + dzt(k,1) * wut(i,1,1) &
3060 + dzt(k,2) * wut(i,1,2) &
3061 + dzt(k,3) * wut(i,1,3) &
3062 + dzt(k,4) * wut(i,1,4) &
3063 + dzt(k,5) * wut(i,1,5) &
3064 + dzt(k,6) * wut(i,1,6) &
3065 + dzt(k,7) * wut(i,1,7)
3066
3067 av(i,1,k,e) = av(i,1,k,e) &
3068 + dzt(k,1) * wvt(i,1,1) &
3069 + dzt(k,2) * wvt(i,1,2) &
3070 + dzt(k,3) * wvt(i,1,3) &
3071 + dzt(k,4) * wvt(i,1,4) &
3072 + dzt(k,5) * wvt(i,1,5) &
3073 + dzt(k,6) * wvt(i,1,6) &
3074 + dzt(k,7) * wvt(i,1,7)
3075
3076 aw(i,1,k,e) = aw(i,1,k,e) &
3077 + dzt(k,1) * wwt(i,1,1) &
3078 + dzt(k,2) * wwt(i,1,2) &
3079 + dzt(k,3) * wwt(i,1,3) &
3080 + dzt(k,4) * wwt(i,1,4) &
3081 + dzt(k,5) * wwt(i,1,5) &
3082 + dzt(k,6) * wwt(i,1,6) &
3083 + dzt(k,7) * wwt(i,1,7)
3084 end do
3085 end do
3086
3087 end do
3088 !$omp end do
3089 end subroutine ax_helm_stress_lx7
3090
3091 subroutine ax_helm_stress_lx6(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
3092 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
3093 jacinv, weights3, n)
3094 integer, parameter :: lx = 6
3095 integer, intent(in) :: n
3096 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3097 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3098 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3099 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3100 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3101 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3102 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3103 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
3104 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
3105 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
3106 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
3107 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
3108 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
3109 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
3110 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
3111 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
3112 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
3113 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
3114 real(kind=rp), intent(in) :: dx(lx, lx)
3115 real(kind=rp), intent(in) :: dy(lx, lx)
3116 real(kind=rp), intent(in) :: dz(lx, lx)
3117 real(kind=rp), intent(in) :: dxt(lx, lx)
3118 real(kind=rp), intent(in) :: dyt(lx, lx)
3119 real(kind=rp), intent(in) :: dzt(lx, lx)
3120 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
3121
3122 real(kind=rp) :: wur(lx, lx, lx)
3123 real(kind=rp) :: wus(lx, lx, lx)
3124 real(kind=rp) :: wut(lx, lx, lx)
3125 real(kind=rp) :: wvr(lx, lx, lx)
3126 real(kind=rp) :: wvs(lx, lx, lx)
3127 real(kind=rp) :: wvt(lx, lx, lx)
3128 real(kind=rp) :: wwr(lx, lx, lx)
3129 real(kind=rp) :: wws(lx, lx, lx)
3130 real(kind=rp) :: wwt(lx, lx, lx)
3131
3132 integer :: e, i, j, k, l
3133 real(kind=rp) :: dj
3134 real(kind=rp) :: s11, s22, s33, s12, s13, s23
3135 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
3136
3137 !$omp do
3138 do e = 1, n
3139
3140 do j = 1, lx * lx
3141 do i = 1, lx
3142 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3143 + dx(i,2) * u(2,j,1,e) &
3144 + dx(i,3) * u(3,j,1,e) &
3145 + dx(i,4) * u(4,j,1,e) &
3146 + dx(i,5) * u(5,j,1,e) &
3147 + dx(i,6) * u(6,j,1,e)
3148
3149 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3150 + dx(i,2) * v(2,j,1,e) &
3151 + dx(i,3) * v(3,j,1,e) &
3152 + dx(i,4) * v(4,j,1,e) &
3153 + dx(i,5) * v(5,j,1,e) &
3154 + dx(i,6) * v(6,j,1,e)
3155
3156 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3157 + dx(i,2) * w(2,j,1,e) &
3158 + dx(i,3) * w(3,j,1,e) &
3159 + dx(i,4) * w(4,j,1,e) &
3160 + dx(i,5) * w(5,j,1,e) &
3161 + dx(i,6) * w(6,j,1,e)
3162 end do
3163 end do
3164
3165 do k = 1, lx
3166 do j = 1, lx
3167 do i = 1, lx
3168 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3169 + dy(j,2) * u(i,2,k,e) &
3170 + dy(j,3) * u(i,3,k,e) &
3171 + dy(j,4) * u(i,4,k,e) &
3172 + dy(j,5) * u(i,5,k,e) &
3173 + dy(j,6) * u(i,6,k,e)
3174
3175 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3176 + dy(j,2) * v(i,2,k,e) &
3177 + dy(j,3) * v(i,3,k,e) &
3178 + dy(j,4) * v(i,4,k,e) &
3179 + dy(j,5) * v(i,5,k,e) &
3180 + dy(j,6) * v(i,6,k,e)
3181
3182 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3183 + dy(j,2) * w(i,2,k,e) &
3184 + dy(j,3) * w(i,3,k,e) &
3185 + dy(j,4) * w(i,4,k,e) &
3186 + dy(j,5) * w(i,5,k,e) &
3187 + dy(j,6) * w(i,6,k,e)
3188 end do
3189 end do
3190 end do
3191
3192 do k = 1, lx
3193 do i = 1, lx*lx
3194 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3195 + dz(k,2) * u(i,1,2,e) &
3196 + dz(k,3) * u(i,1,3,e) &
3197 + dz(k,4) * u(i,1,4,e) &
3198 + dz(k,5) * u(i,1,5,e) &
3199 + dz(k,6) * u(i,1,6,e)
3200
3201 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3202 + dz(k,2) * v(i,1,2,e) &
3203 + dz(k,3) * v(i,1,3,e) &
3204 + dz(k,4) * v(i,1,4,e) &
3205 + dz(k,5) * v(i,1,5,e) &
3206 + dz(k,6) * v(i,1,6,e)
3207
3208 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3209 + dz(k,2) * w(i,1,2,e) &
3210 + dz(k,3) * w(i,1,3,e) &
3211 + dz(k,4) * w(i,1,4,e) &
3212 + dz(k,5) * w(i,1,5,e) &
3213 + dz(k,6) * w(i,1,6,e)
3214 end do
3215 end do
3216
3217 do i = 1, lx*lx*lx
3218
3219 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
3220 + wut(i,1,1) * dtdx(i,1,1,e)
3221 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
3222 + wut(i,1,1) * dtdy(i,1,1,e)
3223 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
3224 + wut(i,1,1) * dtdz(i,1,1,e)
3225
3226 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
3227 + wvt(i,1,1) * dtdx(i,1,1,e)
3228 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
3229 + wvt(i,1,1) * dtdy(i,1,1,e)
3230 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
3231 + wvt(i,1,1) * dtdz(i,1,1,e)
3232
3233 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
3234 + wwt(i,1,1) * dtdx(i,1,1,e)
3235 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
3236 + wwt(i,1,1) * dtdy(i,1,1,e)
3237 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
3238 + wwt(i,1,1) * dtdz(i,1,1,e)
3239
3240 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
3241 s11 = dj*(u1 + u1)
3242 s22 = dj*(v2 + v2)
3243 s33 = dj*(w3 + w3)
3244 s12 = dj*(u2 + v1)
3245 s13 = dj*(u3 + w1)
3246 s23 = dj*(v3 + w2)
3247
3248 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
3249 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
3250 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
3251
3252 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
3253 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
3254 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
3255
3256 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
3257 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
3258 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
3259 end do
3260
3261 do j = 1, lx*lx
3262 do i = 1, lx
3263 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
3264 + dxt(i,2) * wur(2,j,1) &
3265 + dxt(i,3) * wur(3,j,1) &
3266 + dxt(i,4) * wur(4,j,1) &
3267 + dxt(i,5) * wur(5,j,1) &
3268 + dxt(i,6) * wur(6,j,1)
3269
3270 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
3271 + dxt(i,2) * wvr(2,j,1) &
3272 + dxt(i,3) * wvr(3,j,1) &
3273 + dxt(i,4) * wvr(4,j,1) &
3274 + dxt(i,5) * wvr(5,j,1) &
3275 + dxt(i,6) * wvr(6,j,1)
3276
3277 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
3278 + dxt(i,2) * wwr(2,j,1) &
3279 + dxt(i,3) * wwr(3,j,1) &
3280 + dxt(i,4) * wwr(4,j,1) &
3281 + dxt(i,5) * wwr(5,j,1) &
3282 + dxt(i,6) * wwr(6,j,1)
3283 end do
3284 end do
3285
3286 do k = 1, lx
3287 do j = 1, lx
3288 do i = 1, lx
3289 au(i,j,k,e) = au(i,j,k,e) &
3290 + dyt(j,1) * wus(i,1,k) &
3291 + dyt(j,2) * wus(i,2,k) &
3292 + dyt(j,3) * wus(i,3,k) &
3293 + dyt(j,4) * wus(i,4,k) &
3294 + dyt(j,5) * wus(i,5,k) &
3295 + dyt(j,6) * wus(i,6,k)
3296
3297 av(i,j,k,e) = av(i,j,k,e) &
3298 + dyt(j,1) * wvs(i,1,k) &
3299 + dyt(j,2) * wvs(i,2,k) &
3300 + dyt(j,3) * wvs(i,3,k) &
3301 + dyt(j,4) * wvs(i,4,k) &
3302 + dyt(j,5) * wvs(i,5,k) &
3303 + dyt(j,6) * wvs(i,6,k)
3304
3305 aw(i,j,k,e) = aw(i,j,k,e) &
3306 + dyt(j,1) * wws(i,1,k) &
3307 + dyt(j,2) * wws(i,2,k) &
3308 + dyt(j,3) * wws(i,3,k) &
3309 + dyt(j,4) * wws(i,4,k) &
3310 + dyt(j,5) * wws(i,5,k) &
3311 + dyt(j,6) * wws(i,6,k)
3312 end do
3313 end do
3314 end do
3315
3316 do k = 1, lx
3317 do i = 1, lx*lx
3318 au(i,1,k,e) = au(i,1,k,e) &
3319 + dzt(k,1) * wut(i,1,1) &
3320 + dzt(k,2) * wut(i,1,2) &
3321 + dzt(k,3) * wut(i,1,3) &
3322 + dzt(k,4) * wut(i,1,4) &
3323 + dzt(k,5) * wut(i,1,5) &
3324 + dzt(k,6) * wut(i,1,6)
3325
3326 av(i,1,k,e) = av(i,1,k,e) &
3327 + dzt(k,1) * wvt(i,1,1) &
3328 + dzt(k,2) * wvt(i,1,2) &
3329 + dzt(k,3) * wvt(i,1,3) &
3330 + dzt(k,4) * wvt(i,1,4) &
3331 + dzt(k,5) * wvt(i,1,5) &
3332 + dzt(k,6) * wvt(i,1,6)
3333
3334 aw(i,1,k,e) = aw(i,1,k,e) &
3335 + dzt(k,1) * wwt(i,1,1) &
3336 + dzt(k,2) * wwt(i,1,2) &
3337 + dzt(k,3) * wwt(i,1,3) &
3338 + dzt(k,4) * wwt(i,1,4) &
3339 + dzt(k,5) * wwt(i,1,5) &
3340 + dzt(k,6) * wwt(i,1,6)
3341 end do
3342 end do
3343
3344 end do
3345 !$omp end do
3346 end subroutine ax_helm_stress_lx6
3347
3348 subroutine ax_helm_stress_lx5(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
3349 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
3350 jacinv, weights3, n)
3351 integer, parameter :: lx = 5
3352 integer, intent(in) :: n
3353 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3354 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3355 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3356 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3357 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3358 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3359 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3360 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
3361 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
3362 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
3363 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
3364 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
3365 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
3366 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
3367 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
3368 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
3369 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
3370 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
3371 real(kind=rp), intent(in) :: dx(lx, lx)
3372 real(kind=rp), intent(in) :: dy(lx, lx)
3373 real(kind=rp), intent(in) :: dz(lx, lx)
3374 real(kind=rp), intent(in) :: dxt(lx, lx)
3375 real(kind=rp), intent(in) :: dyt(lx, lx)
3376 real(kind=rp), intent(in) :: dzt(lx, lx)
3377 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
3378
3379 real(kind=rp) :: wur(lx, lx, lx)
3380 real(kind=rp) :: wus(lx, lx, lx)
3381 real(kind=rp) :: wut(lx, lx, lx)
3382 real(kind=rp) :: wvr(lx, lx, lx)
3383 real(kind=rp) :: wvs(lx, lx, lx)
3384 real(kind=rp) :: wvt(lx, lx, lx)
3385 real(kind=rp) :: wwr(lx, lx, lx)
3386 real(kind=rp) :: wws(lx, lx, lx)
3387 real(kind=rp) :: wwt(lx, lx, lx)
3388
3389 integer :: e, i, j, k, l
3390 real(kind=rp) :: dj
3391 real(kind=rp) :: s11, s22, s33, s12, s13, s23
3392 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
3393
3394 !$omp do
3395 do e = 1, n
3396
3397 do j = 1, lx * lx
3398 do i = 1, lx
3399 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3400 + dx(i,2) * u(2,j,1,e) &
3401 + dx(i,3) * u(3,j,1,e) &
3402 + dx(i,4) * u(4,j,1,e) &
3403 + dx(i,5) * u(5,j,1,e)
3404
3405 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3406 + dx(i,2) * v(2,j,1,e) &
3407 + dx(i,3) * v(3,j,1,e) &
3408 + dx(i,4) * v(4,j,1,e) &
3409 + dx(i,5) * v(5,j,1,e)
3410
3411 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3412 + dx(i,2) * w(2,j,1,e) &
3413 + dx(i,3) * w(3,j,1,e) &
3414 + dx(i,4) * w(4,j,1,e) &
3415 + dx(i,5) * w(5,j,1,e)
3416 end do
3417 end do
3418
3419 do k = 1, lx
3420 do j = 1, lx
3421 do i = 1, lx
3422 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3423 + dy(j,2) * u(i,2,k,e) &
3424 + dy(j,3) * u(i,3,k,e) &
3425 + dy(j,4) * u(i,4,k,e) &
3426 + dy(j,5) * u(i,5,k,e)
3427
3428 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3429 + dy(j,2) * v(i,2,k,e) &
3430 + dy(j,3) * v(i,3,k,e) &
3431 + dy(j,4) * v(i,4,k,e) &
3432 + dy(j,5) * v(i,5,k,e)
3433
3434 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3435 + dy(j,2) * w(i,2,k,e) &
3436 + dy(j,3) * w(i,3,k,e) &
3437 + dy(j,4) * w(i,4,k,e) &
3438 + dy(j,5) * w(i,5,k,e)
3439 end do
3440 end do
3441 end do
3442
3443 do k = 1, lx
3444 do i = 1, lx*lx
3445 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3446 + dz(k,2) * u(i,1,2,e) &
3447 + dz(k,3) * u(i,1,3,e) &
3448 + dz(k,4) * u(i,1,4,e) &
3449 + dz(k,5) * u(i,1,5,e)
3450
3451 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3452 + dz(k,2) * v(i,1,2,e) &
3453 + dz(k,3) * v(i,1,3,e) &
3454 + dz(k,4) * v(i,1,4,e) &
3455 + dz(k,5) * v(i,1,5,e)
3456
3457 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3458 + dz(k,2) * w(i,1,2,e) &
3459 + dz(k,3) * w(i,1,3,e) &
3460 + dz(k,4) * w(i,1,4,e) &
3461 + dz(k,5) * w(i,1,5,e)
3462 end do
3463 end do
3464
3465 do i = 1, lx*lx*lx
3466
3467 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
3468 + wut(i,1,1) * dtdx(i,1,1,e)
3469 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
3470 + wut(i,1,1) * dtdy(i,1,1,e)
3471 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
3472 + wut(i,1,1) * dtdz(i,1,1,e)
3473
3474 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
3475 + wvt(i,1,1) * dtdx(i,1,1,e)
3476 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
3477 + wvt(i,1,1) * dtdy(i,1,1,e)
3478 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
3479 + wvt(i,1,1) * dtdz(i,1,1,e)
3480
3481 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
3482 + wwt(i,1,1) * dtdx(i,1,1,e)
3483 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
3484 + wwt(i,1,1) * dtdy(i,1,1,e)
3485 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
3486 + wwt(i,1,1) * dtdz(i,1,1,e)
3487
3488 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
3489 s11 = dj*(u1 + u1)
3490 s22 = dj*(v2 + v2)
3491 s33 = dj*(w3 + w3)
3492 s12 = dj*(u2 + v1)
3493 s13 = dj*(u3 + w1)
3494 s23 = dj*(v3 + w2)
3495
3496 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
3497 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
3498 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
3499
3500 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
3501 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
3502 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
3503
3504 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
3505 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
3506 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
3507 end do
3508
3509 do j = 1, lx*lx
3510 do i = 1, lx
3511 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
3512 + dxt(i,2) * wur(2,j,1) &
3513 + dxt(i,3) * wur(3,j,1) &
3514 + dxt(i,4) * wur(4,j,1) &
3515 + dxt(i,5) * wur(5,j,1)
3516
3517 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
3518 + dxt(i,2) * wvr(2,j,1) &
3519 + dxt(i,3) * wvr(3,j,1) &
3520 + dxt(i,4) * wvr(4,j,1) &
3521 + dxt(i,5) * wvr(5,j,1)
3522
3523 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
3524 + dxt(i,2) * wwr(2,j,1) &
3525 + dxt(i,3) * wwr(3,j,1) &
3526 + dxt(i,4) * wwr(4,j,1) &
3527 + dxt(i,5) * wwr(5,j,1)
3528
3529 end do
3530 end do
3531
3532 do k = 1, lx
3533 do j = 1, lx
3534 do i = 1, lx
3535 au(i,j,k,e) = au(i,j,k,e) &
3536 + dyt(j,1) * wus(i,1,k) &
3537 + dyt(j,2) * wus(i,2,k) &
3538 + dyt(j,3) * wus(i,3,k) &
3539 + dyt(j,4) * wus(i,4,k) &
3540 + dyt(j,5) * wus(i,5,k)
3541
3542 av(i,j,k,e) = av(i,j,k,e) &
3543 + dyt(j,1) * wvs(i,1,k) &
3544 + dyt(j,2) * wvs(i,2,k) &
3545 + dyt(j,3) * wvs(i,3,k) &
3546 + dyt(j,4) * wvs(i,4,k) &
3547 + dyt(j,5) * wvs(i,5,k)
3548
3549 aw(i,j,k,e) = aw(i,j,k,e) &
3550 + dyt(j,1) * wws(i,1,k) &
3551 + dyt(j,2) * wws(i,2,k) &
3552 + dyt(j,3) * wws(i,3,k) &
3553 + dyt(j,4) * wws(i,4,k) &
3554 + dyt(j,5) * wws(i,5,k)
3555 end do
3556 end do
3557 end do
3558
3559 do k = 1, lx
3560 do i = 1, lx*lx
3561 au(i,1,k,e) = au(i,1,k,e) &
3562 + dzt(k,1) * wut(i,1,1) &
3563 + dzt(k,2) * wut(i,1,2) &
3564 + dzt(k,3) * wut(i,1,3) &
3565 + dzt(k,4) * wut(i,1,4) &
3566 + dzt(k,5) * wut(i,1,5)
3567
3568 av(i,1,k,e) = av(i,1,k,e) &
3569 + dzt(k,1) * wvt(i,1,1) &
3570 + dzt(k,2) * wvt(i,1,2) &
3571 + dzt(k,3) * wvt(i,1,3) &
3572 + dzt(k,4) * wvt(i,1,4) &
3573 + dzt(k,5) * wvt(i,1,5)
3574
3575 aw(i,1,k,e) = aw(i,1,k,e) &
3576 + dzt(k,1) * wwt(i,1,1) &
3577 + dzt(k,2) * wwt(i,1,2) &
3578 + dzt(k,3) * wwt(i,1,3) &
3579 + dzt(k,4) * wwt(i,1,4) &
3580 + dzt(k,5) * wwt(i,1,5)
3581 end do
3582 end do
3583
3584 end do
3585 !$omp end do
3586 end subroutine ax_helm_stress_lx5
3587
3588 subroutine ax_helm_stress_lx4(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
3589 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
3590 jacinv, weights3, n)
3591 integer, parameter :: lx = 4
3592 integer, intent(in) :: n
3593 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3594 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3595 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3596 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3597 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3598 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3599 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3600 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
3601 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
3602 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
3603 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
3604 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
3605 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
3606 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
3607 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
3608 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
3609 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
3610 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
3611 real(kind=rp), intent(in) :: dx(lx, lx)
3612 real(kind=rp), intent(in) :: dy(lx, lx)
3613 real(kind=rp), intent(in) :: dz(lx, lx)
3614 real(kind=rp), intent(in) :: dxt(lx, lx)
3615 real(kind=rp), intent(in) :: dyt(lx, lx)
3616 real(kind=rp), intent(in) :: dzt(lx, lx)
3617 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
3618
3619 real(kind=rp) :: wur(lx, lx, lx)
3620 real(kind=rp) :: wus(lx, lx, lx)
3621 real(kind=rp) :: wut(lx, lx, lx)
3622 real(kind=rp) :: wvr(lx, lx, lx)
3623 real(kind=rp) :: wvs(lx, lx, lx)
3624 real(kind=rp) :: wvt(lx, lx, lx)
3625 real(kind=rp) :: wwr(lx, lx, lx)
3626 real(kind=rp) :: wws(lx, lx, lx)
3627 real(kind=rp) :: wwt(lx, lx, lx)
3628
3629 integer :: e, i, j, k, l
3630 real(kind=rp) :: dj
3631 real(kind=rp) :: s11, s22, s33, s12, s13, s23
3632 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
3633
3634 !$omp do
3635 do e = 1, n
3636
3637 do j = 1, lx * lx
3638 do i = 1, lx
3639 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3640 + dx(i,2) * u(2,j,1,e) &
3641 + dx(i,3) * u(3,j,1,e) &
3642 + dx(i,4) * u(4,j,1,e)
3643
3644 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3645 + dx(i,2) * v(2,j,1,e) &
3646 + dx(i,3) * v(3,j,1,e) &
3647 + dx(i,4) * v(4,j,1,e)
3648
3649 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3650 + dx(i,2) * w(2,j,1,e) &
3651 + dx(i,3) * w(3,j,1,e) &
3652 + dx(i,4) * w(4,j,1,e)
3653 end do
3654 end do
3655
3656 do k = 1, lx
3657 do j = 1, lx
3658 do i = 1, lx
3659 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3660 + dy(j,2) * u(i,2,k,e) &
3661 + dy(j,3) * u(i,3,k,e) &
3662 + dy(j,4) * u(i,4,k,e)
3663
3664 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3665 + dy(j,2) * v(i,2,k,e) &
3666 + dy(j,3) * v(i,3,k,e) &
3667 + dy(j,4) * v(i,4,k,e)
3668
3669 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3670 + dy(j,2) * w(i,2,k,e) &
3671 + dy(j,3) * w(i,3,k,e) &
3672 + dy(j,4) * w(i,4,k,e)
3673 end do
3674 end do
3675 end do
3676
3677 do k = 1, lx
3678 do i = 1, lx*lx
3679 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3680 + dz(k,2) * u(i,1,2,e) &
3681 + dz(k,3) * u(i,1,3,e) &
3682 + dz(k,4) * u(i,1,4,e)
3683
3684 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3685 + dz(k,2) * v(i,1,2,e) &
3686 + dz(k,3) * v(i,1,3,e) &
3687 + dz(k,4) * v(i,1,4,e)
3688
3689 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3690 + dz(k,2) * w(i,1,2,e) &
3691 + dz(k,3) * w(i,1,3,e) &
3692 + dz(k,4) * w(i,1,4,e)
3693 end do
3694 end do
3695
3696 do i = 1, lx*lx*lx
3697
3698 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
3699 + wut(i,1,1) * dtdx(i,1,1,e)
3700 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
3701 + wut(i,1,1) * dtdy(i,1,1,e)
3702 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
3703 + wut(i,1,1) * dtdz(i,1,1,e)
3704
3705 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
3706 + wvt(i,1,1) * dtdx(i,1,1,e)
3707 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
3708 + wvt(i,1,1) * dtdy(i,1,1,e)
3709 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
3710 + wvt(i,1,1) * dtdz(i,1,1,e)
3711
3712 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
3713 + wwt(i,1,1) * dtdx(i,1,1,e)
3714 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
3715 + wwt(i,1,1) * dtdy(i,1,1,e)
3716 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
3717 + wwt(i,1,1) * dtdz(i,1,1,e)
3718
3719 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
3720 s11 = dj*(u1 + u1)
3721 s22 = dj*(v2 + v2)
3722 s33 = dj*(w3 + w3)
3723 s12 = dj*(u2 + v1)
3724 s13 = dj*(u3 + w1)
3725 s23 = dj*(v3 + w2)
3726
3727 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
3728 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
3729 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
3730
3731 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
3732 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
3733 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
3734
3735 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
3736 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
3737 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
3738 end do
3739
3740 do j = 1, lx*lx
3741 do i = 1, lx
3742 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
3743 + dxt(i,2) * wur(2,j,1) &
3744 + dxt(i,3) * wur(3,j,1) &
3745 + dxt(i,4) * wur(4,j,1)
3746
3747 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
3748 + dxt(i,2) * wvr(2,j,1) &
3749 + dxt(i,3) * wvr(3,j,1) &
3750 + dxt(i,4) * wvr(4,j,1)
3751
3752 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
3753 + dxt(i,2) * wwr(2,j,1) &
3754 + dxt(i,3) * wwr(3,j,1) &
3755 + dxt(i,4) * wwr(4,j,1)
3756 end do
3757 end do
3758
3759 do k = 1, lx
3760 do j = 1, lx
3761 do i = 1, lx
3762 au(i,j,k,e) = au(i,j,k,e) &
3763 + dyt(j,1) * wus(i,1,k) &
3764 + dyt(j,2) * wus(i,2,k) &
3765 + dyt(j,3) * wus(i,3,k) &
3766 + dyt(j,4) * wus(i,4,k)
3767
3768 av(i,j,k,e) = av(i,j,k,e) &
3769 + dyt(j,1) * wvs(i,1,k) &
3770 + dyt(j,2) * wvs(i,2,k) &
3771 + dyt(j,3) * wvs(i,3,k) &
3772 + dyt(j,4) * wvs(i,4,k)
3773
3774 aw(i,j,k,e) = aw(i,j,k,e) &
3775 + dyt(j,1) * wws(i,1,k) &
3776 + dyt(j,2) * wws(i,2,k) &
3777 + dyt(j,3) * wws(i,3,k) &
3778 + dyt(j,4) * wws(i,4,k)
3779 end do
3780 end do
3781 end do
3782
3783 do k = 1, lx
3784 do i = 1, lx*lx
3785 au(i,1,k,e) = au(i,1,k,e) &
3786 + dzt(k,1) * wut(i,1,1) &
3787 + dzt(k,2) * wut(i,1,2) &
3788 + dzt(k,3) * wut(i,1,3) &
3789 + dzt(k,4) * wut(i,1,4)
3790
3791 av(i,1,k,e) = av(i,1,k,e) &
3792 + dzt(k,1) * wvt(i,1,1) &
3793 + dzt(k,2) * wvt(i,1,2) &
3794 + dzt(k,3) * wvt(i,1,3) &
3795 + dzt(k,4) * wvt(i,1,4)
3796
3797 aw(i,1,k,e) = aw(i,1,k,e) &
3798 + dzt(k,1) * wwt(i,1,1) &
3799 + dzt(k,2) * wwt(i,1,2) &
3800 + dzt(k,3) * wwt(i,1,3) &
3801 + dzt(k,4) * wwt(i,1,4)
3802 end do
3803 end do
3804
3805 end do
3806 !$omp end do
3807 end subroutine ax_helm_stress_lx4
3808
3809 subroutine ax_helm_stress_lx3(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
3810 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
3811 jacinv, weights3, n)
3812 integer, parameter :: lx = 3
3813 integer, intent(in) :: n
3814 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3815 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3816 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3817 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3818 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3819 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3820 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3821 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
3822 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
3823 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
3824 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
3825 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
3826 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
3827 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
3828 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
3829 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
3830 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
3831 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
3832 real(kind=rp), intent(in) :: dx(lx, lx)
3833 real(kind=rp), intent(in) :: dy(lx, lx)
3834 real(kind=rp), intent(in) :: dz(lx, lx)
3835 real(kind=rp), intent(in) :: dxt(lx, lx)
3836 real(kind=rp), intent(in) :: dyt(lx, lx)
3837 real(kind=rp), intent(in) :: dzt(lx, lx)
3838 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
3839
3840 real(kind=rp) :: wur(lx, lx, lx)
3841 real(kind=rp) :: wus(lx, lx, lx)
3842 real(kind=rp) :: wut(lx, lx, lx)
3843 real(kind=rp) :: wvr(lx, lx, lx)
3844 real(kind=rp) :: wvs(lx, lx, lx)
3845 real(kind=rp) :: wvt(lx, lx, lx)
3846 real(kind=rp) :: wwr(lx, lx, lx)
3847 real(kind=rp) :: wws(lx, lx, lx)
3848 real(kind=rp) :: wwt(lx, lx, lx)
3849
3850 integer :: e, i, j, k, l
3851 real(kind=rp) :: dj
3852 real(kind=rp) :: s11, s22, s33, s12, s13, s23
3853 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
3854
3855 !$omp do
3856 do e = 1, n
3857
3858 do j = 1, lx * lx
3859 do i = 1, lx
3860 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3861 + dx(i,2) * u(2,j,1,e) &
3862 + dx(i,3) * u(3,j,1,e)
3863
3864 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3865 + dx(i,2) * v(2,j,1,e) &
3866 + dx(i,3) * v(3,j,1,e)
3867
3868 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3869 + dx(i,2) * w(2,j,1,e) &
3870 + dx(i,3) * w(3,j,1,e)
3871 end do
3872 end do
3873
3874 do k = 1, lx
3875 do j = 1, lx
3876 do i = 1, lx
3877 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3878 + dy(j,2) * u(i,2,k,e) &
3879 + dy(j,3) * u(i,3,k,e)
3880
3881 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3882 + dy(j,2) * v(i,2,k,e) &
3883 + dy(j,3) * v(i,3,k,e)
3884
3885 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3886 + dy(j,2) * w(i,2,k,e) &
3887 + dy(j,3) * w(i,3,k,e)
3888 end do
3889 end do
3890 end do
3891
3892 do k = 1, lx
3893 do i = 1, lx*lx
3894 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3895 + dz(k,2) * u(i,1,2,e) &
3896 + dz(k,3) * u(i,1,3,e)
3897
3898 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3899 + dz(k,2) * v(i,1,2,e) &
3900 + dz(k,3) * v(i,1,3,e)
3901
3902 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3903 + dz(k,2) * w(i,1,2,e) &
3904 + dz(k,3) * w(i,1,3,e)
3905 end do
3906 end do
3907
3908 do i = 1, lx*lx*lx
3909
3910 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
3911 + wut(i,1,1) * dtdx(i,1,1,e)
3912 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
3913 + wut(i,1,1) * dtdy(i,1,1,e)
3914 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
3915 + wut(i,1,1) * dtdz(i,1,1,e)
3916
3917 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
3918 + wvt(i,1,1) * dtdx(i,1,1,e)
3919 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
3920 + wvt(i,1,1) * dtdy(i,1,1,e)
3921 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
3922 + wvt(i,1,1) * dtdz(i,1,1,e)
3923
3924 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
3925 + wwt(i,1,1) * dtdx(i,1,1,e)
3926 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
3927 + wwt(i,1,1) * dtdy(i,1,1,e)
3928 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
3929 + wwt(i,1,1) * dtdz(i,1,1,e)
3930
3931 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
3932 s11 = dj*(u1 + u1)
3933 s22 = dj*(v2 + v2)
3934 s33 = dj*(w3 + w3)
3935 s12 = dj*(u2 + v1)
3936 s13 = dj*(u3 + w1)
3937 s23 = dj*(v3 + w2)
3938
3939 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
3940 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
3941 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
3942
3943 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
3944 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
3945 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
3946
3947 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
3948 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
3949 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
3950 end do
3951
3952 do j = 1, lx*lx
3953 do i = 1, lx
3954 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
3955 + dxt(i,2) * wur(2,j,1) &
3956 + dxt(i,3) * wur(3,j,1)
3957
3958 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
3959 + dxt(i,2) * wvr(2,j,1) &
3960 + dxt(i,3) * wvr(3,j,1)
3961
3962 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
3963 + dxt(i,2) * wwr(2,j,1) &
3964 + dxt(i,3) * wwr(3,j,1)
3965 end do
3966 end do
3967
3968 do k = 1, lx
3969 do j = 1, lx
3970 do i = 1, lx
3971 au(i,j,k,e) = au(i,j,k,e) &
3972 + dyt(j,1) * wus(i,1,k) &
3973 + dyt(j,2) * wus(i,2,k) &
3974 + dyt(j,3) * wus(i,3,k)
3975
3976 av(i,j,k,e) = av(i,j,k,e) &
3977 + dyt(j,1) * wvs(i,1,k) &
3978 + dyt(j,2) * wvs(i,2,k) &
3979 + dyt(j,3) * wvs(i,3,k)
3980
3981 aw(i,j,k,e) = aw(i,j,k,e) &
3982 + dyt(j,1) * wws(i,1,k) &
3983 + dyt(j,2) * wws(i,2,k) &
3984 + dyt(j,3) * wws(i,3,k)
3985 end do
3986 end do
3987 end do
3988
3989 do k = 1, lx
3990 do i = 1, lx*lx
3991 au(i,1,k,e) = au(i,1,k,e) &
3992 + dzt(k,1) * wut(i,1,1) &
3993 + dzt(k,2) * wut(i,1,2) &
3994 + dzt(k,3) * wut(i,1,3)
3995
3996 av(i,1,k,e) = av(i,1,k,e) &
3997 + dzt(k,1) * wvt(i,1,1) &
3998 + dzt(k,2) * wvt(i,1,2) &
3999 + dzt(k,3) * wvt(i,1,3)
4000
4001 aw(i,1,k,e) = aw(i,1,k,e) &
4002 + dzt(k,1) * wwt(i,1,1) &
4003 + dzt(k,2) * wwt(i,1,2) &
4004 + dzt(k,3) * wwt(i,1,3)
4005 end do
4006 end do
4007
4008 end do
4009 !$omp end do
4010 end subroutine ax_helm_stress_lx3
4011
4012 subroutine ax_helm_stress_lx2(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, &
4013 Dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
4014 jacinv, weights3, n)
4015 integer, parameter :: lx = 2
4016 integer, intent(in) :: n
4017 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
4018 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
4019 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
4020 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
4021 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
4022 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
4023 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
4024 real(kind=rp), intent(in) :: h2(lx, lx, lx, n)
4025 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
4026 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
4027 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
4028 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
4029 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
4030 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
4031 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
4032 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
4033 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
4034 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
4035 real(kind=rp), intent(in) :: dx(lx, lx)
4036 real(kind=rp), intent(in) :: dy(lx, lx)
4037 real(kind=rp), intent(in) :: dz(lx, lx)
4038 real(kind=rp), intent(in) :: dxt(lx, lx)
4039 real(kind=rp), intent(in) :: dyt(lx, lx)
4040 real(kind=rp), intent(in) :: dzt(lx, lx)
4041 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
4042
4043 real(kind=rp) :: wur(lx, lx, lx)
4044 real(kind=rp) :: wus(lx, lx, lx)
4045 real(kind=rp) :: wut(lx, lx, lx)
4046 real(kind=rp) :: wvr(lx, lx, lx)
4047 real(kind=rp) :: wvs(lx, lx, lx)
4048 real(kind=rp) :: wvt(lx, lx, lx)
4049 real(kind=rp) :: wwr(lx, lx, lx)
4050 real(kind=rp) :: wws(lx, lx, lx)
4051 real(kind=rp) :: wwt(lx, lx, lx)
4052
4053 integer :: e, i, j, k
4054 real(kind=rp) :: dj
4055 real(kind=rp) :: s11, s22, s33, s12, s13, s23
4056 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
4057
4058 !$omp do
4059 do e = 1, n
4060
4061 do j = 1, lx * lx
4062 do i = 1, lx
4063 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
4064 + dx(i,2) * u(2,j,1,e)
4065
4066 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
4067 + dx(i,2) * v(2,j,1,e)
4068
4069 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
4070 + dx(i,2) * w(2,j,1,e)
4071 end do
4072 end do
4073
4074 do k = 1, lx
4075 do j = 1, lx
4076 do i = 1, lx
4077 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
4078 + dy(j,2) * u(i,2,k,e)
4079
4080 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
4081 + dy(j,2) * v(i,2,k,e)
4082
4083 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
4084 + dy(j,2) * w(i,2,k,e)
4085 end do
4086 end do
4087 end do
4088
4089 do k = 1, lx
4090 do i = 1, lx*lx
4091 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
4092 + dz(k,2) * u(i,1,2,e)
4093
4094 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
4095 + dz(k,2) * v(i,1,2,e)
4096
4097 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
4098 + dz(k,2) * w(i,1,2,e)
4099 end do
4100 end do
4101
4102 do i = 1, lx*lx*lx
4103
4104 u1 = wur(i,1,1) * drdx(i,1,1,e) + wus(i,1,1) * dsdx(i,1,1,e) &
4105 + wut(i,1,1) * dtdx(i,1,1,e)
4106 u2 = wur(i,1,1) * drdy(i,1,1,e) + wus(i,1,1) * dsdy(i,1,1,e) &
4107 + wut(i,1,1) * dtdy(i,1,1,e)
4108 u3 = wur(i,1,1) * drdz(i,1,1,e) + wus(i,1,1) * dsdz(i,1,1,e) &
4109 + wut(i,1,1) * dtdz(i,1,1,e)
4110
4111 v1 = wvr(i,1,1) * drdx(i,1,1,e) + wvs(i,1,1) * dsdx(i,1,1,e) &
4112 + wvt(i,1,1) * dtdx(i,1,1,e)
4113 v2 = wvr(i,1,1) * drdy(i,1,1,e) + wvs(i,1,1) * dsdy(i,1,1,e) &
4114 + wvt(i,1,1) * dtdy(i,1,1,e)
4115 v3 = wvr(i,1,1) * drdz(i,1,1,e) + wvs(i,1,1) * dsdz(i,1,1,e) &
4116 + wvt(i,1,1) * dtdz(i,1,1,e)
4117
4118 w1 = wwr(i,1,1) * drdx(i,1,1,e) + wws(i,1,1) * dsdx(i,1,1,e) &
4119 + wwt(i,1,1) * dtdx(i,1,1,e)
4120 w2 = wwr(i,1,1) * drdy(i,1,1,e) + wws(i,1,1) * dsdy(i,1,1,e) &
4121 + wwt(i,1,1) * dtdy(i,1,1,e)
4122 w3 = wwr(i,1,1) * drdz(i,1,1,e) + wws(i,1,1) * dsdz(i,1,1,e) &
4123 + wwt(i,1,1) * dtdz(i,1,1,e)
4124
4125 dj = h1(i,1,1,e) * weights3(i,1,1) * jacinv(i,1,1,e)
4126 s11 = dj*(u1 + u1)
4127 s22 = dj*(v2 + v2)
4128 s33 = dj*(w3 + w3)
4129 s12 = dj*(u2 + v1)
4130 s13 = dj*(u3 + w1)
4131 s23 = dj*(v3 + w2)
4132
4133 wur(i,1,1) = drdx(i,1,1,e)*s11 + drdy(i,1,1,e)*s12 + drdz(i,1,1,e)*s13
4134 wus(i,1,1) = dsdx(i,1,1,e)*s11 + dsdy(i,1,1,e)*s12 + dsdz(i,1,1,e)*s13
4135 wut(i,1,1) = dtdx(i,1,1,e)*s11 + dtdy(i,1,1,e)*s12 + dtdz(i,1,1,e)*s13
4136
4137 wvr(i,1,1) = drdx(i,1,1,e)*s12 + drdy(i,1,1,e)*s22 + drdz(i,1,1,e)*s23
4138 wvs(i,1,1) = dsdx(i,1,1,e)*s12 + dsdy(i,1,1,e)*s22 + dsdz(i,1,1,e)*s23
4139 wvt(i,1,1) = dtdx(i,1,1,e)*s12 + dtdy(i,1,1,e)*s22 + dtdz(i,1,1,e)*s23
4140
4141 wwr(i,1,1) = drdx(i,1,1,e)*s13 + drdy(i,1,1,e)*s23 + drdz(i,1,1,e)*s33
4142 wws(i,1,1) = dsdx(i,1,1,e)*s13 + dsdy(i,1,1,e)*s23 + dsdz(i,1,1,e)*s33
4143 wwt(i,1,1) = dtdx(i,1,1,e)*s13 + dtdy(i,1,1,e)*s23 + dtdz(i,1,1,e)*s33
4144 end do
4145
4146 do j = 1, lx*lx
4147 do i = 1, lx
4148 au(i,j,1,e) = dxt(i,1) * wur(1,j,1) &
4149 + dxt(i,2) * wur(2,j,1)
4150
4151 av(i,j,1,e) = dxt(i,1) * wvr(1,j,1) &
4152 + dxt(i,2) * wvr(2,j,1)
4153
4154 aw(i,j,1,e) = dxt(i,1) * wwr(1,j,1) &
4155 + dxt(i,2) * wwr(2,j,1)
4156 end do
4157 end do
4158
4159 do k = 1, lx
4160 do j = 1, lx
4161 do i = 1, lx
4162 au(i,j,k,e) = au(i,j,k,e) &
4163 + dyt(j,1) * wus(i,1,k) &
4164 + dyt(j,2) * wus(i,2,k)
4165
4166 av(i,j,k,e) = av(i,j,k,e) &
4167 + dyt(j,1) * wvs(i,1,k) &
4168 + dyt(j,2) * wvs(i,2,k)
4169
4170 aw(i,j,k,e) = aw(i,j,k,e) &
4171 + dyt(j,1) * wws(i,1,k) &
4172 + dyt(j,2) * wws(i,2,k)
4173 end do
4174 end do
4175 end do
4176
4177 do k = 1, lx
4178 do i = 1, lx*lx
4179 au(i,1,k,e) = au(i,1,k,e) &
4180 + dzt(k,1) * wut(i,1,1) &
4181 + dzt(k,2) * wut(i,1,2)
4182
4183 av(i,1,k,e) = av(i,1,k,e) &
4184 + dzt(k,1) * wvt(i,1,1) &
4185 + dzt(k,2) * wvt(i,1,2)
4186
4187 aw(i,1,k,e) = aw(i,1,k,e) &
4188 + dzt(k,1) * wwt(i,1,1) &
4189 + dzt(k,2) * wwt(i,1,2)
4190 end do
4191 end do
4192
4193 end do
4194 !$omp end do
4195 end subroutine ax_helm_stress_lx2
4196
4197end module ax_helm_full_cpu
ax_helm_full_cpu
Definition ax_helm_full_cpu.f90:33
ax_helm_full_cpu::ax_helm_full_compute_vector
subroutine ax_helm_full_compute_vector(this, au, av, aw, u, v, w, coef, msh, xh)
Compute inside a Krylov method, taking 3 components of a vector field in a coupled manner.
Definition ax_helm_full_cpu.f90:64
ax_helm_full_cpu::ax_helm_stress_lx10
subroutine ax_helm_stress_lx10(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:1883
ax_helm_full_cpu::ax_helm_stress_lx7
subroutine ax_helm_stress_lx7(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:2817
ax_helm_full_cpu::ax_helm_stress_lx
subroutine ax_helm_stress_lx(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n, lx)
Definition ax_helm_full_cpu.f90:184
ax_helm_full_cpu::ax_helm_stress_lx6
subroutine ax_helm_stress_lx6(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:3094
ax_helm_full_cpu::ax_helm_stress_lx11
subroutine ax_helm_stress_lx11(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:1535
ax_helm_full_cpu::ax_helm_stress_lx4
subroutine ax_helm_stress_lx4(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:3591
ax_helm_full_cpu::ax_helm_stress_lx2
subroutine ax_helm_stress_lx2(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:4015
ax_helm_full_cpu::ax_helm_stress_lx12
subroutine ax_helm_stress_lx12(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:1168
ax_helm_full_cpu::ax_helm_stress_lx5
subroutine ax_helm_stress_lx5(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:3351
ax_helm_full_cpu::ax_helm_stress_lx13
subroutine ax_helm_stress_lx13(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:784
ax_helm_full_cpu::ax_helm_stress_lx8
subroutine ax_helm_stress_lx8(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:2523
ax_helm_full_cpu::ax_helm_stress_lx14
subroutine ax_helm_stress_lx14(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:381
ax_helm_full_cpu::ax_helm_stress_lx9
subroutine ax_helm_stress_lx9(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:2212
ax_helm_full_cpu::ax_helm_stress_lx3
subroutine ax_helm_stress_lx3(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, h2, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, n)
Definition ax_helm_full_cpu.f90:3812
ax_helm_full
Definition ax_helm_full.f90:33
coefs
Coefficients.
Definition coef.f90:34
mesh
Defines a mesh.
Definition mesh.f90:34
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
space
Defines a function space.
Definition space.f90:34
ax_helm_full::ax_helm_full_t
Matrix-vector product for a Helmholtz problem.
Definition ax_helm_full.f90:45
ax_helm_full_cpu::ax_helm_full_cpu_t
CPU matrix-vector product for a Helmholtz problem with full stress tensor.
Definition ax_helm_full_cpu.f90:43
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:61
mesh::mesh_t
Definition mesh.f90:71
space::space_t
The function space for the SEM solution fields.
Definition space.f90:63