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