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