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