Neko 1.99.9
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
ax_helm_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 2021-2026, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
33module ax_helm_cpu
34 use ax_helm, only : ax_helm_t
35 use num_types, only : rp
36 use coefs, only : coef_t
37 use space, only : space_t
38 use mesh, only : mesh_t
39 implicit none
40 private
41
43 type, public, extends(ax_helm_t) :: ax_helm_cpu_t
44 contains
46 procedure, pass(this) :: compute => ax_helm_compute
48 procedure, pass(this) :: compute_vector => ax_helm_compute_vector
49 end type ax_helm_cpu_t
50
51 interface
52
63 module subroutine ax_helm_compute_vector(this, au, av, aw, &
64 u, v, w, coef, msh, xh)
65 class(ax_helm_cpu_t), intent(in) :: this
66 type(mesh_t), intent(in) :: msh
67 type(space_t), intent(in) :: xh
68 type(coef_t), intent(in) :: coef
69 real(kind=rp), intent(inout) :: au(xh%lx, xh%ly, xh%lz, msh%nelv)
70 real(kind=rp), intent(inout) :: av(xh%lx, xh%ly, xh%lz, msh%nelv)
71 real(kind=rp), intent(inout) :: aw(xh%lx, xh%ly, xh%lz, msh%nelv)
72 real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, msh%nelv)
73 real(kind=rp), intent(in) :: v(xh%lx, xh%ly, xh%lz, msh%nelv)
74 real(kind=rp), intent(in) :: w(xh%lx, xh%ly, xh%lz, msh%nelv)
75 end subroutine ax_helm_compute_vector
76 end interface
77
78contains
79
88 subroutine ax_helm_compute(this, w, u, coef, msh, Xh)
89 class(ax_helm_cpu_t), intent(in) :: this
90 type(mesh_t), intent(in) :: msh
91 type(space_t), intent(in) :: Xh
92 type(coef_t), intent(in) :: coef
93 real(kind=rp), intent(inout) :: w(xh%lx, xh%ly, xh%lz, msh%nelv)
94 real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, msh%nelv)
95 integer :: i
96
97 !$omp parallel
98 select case(xh%lx)
99 case (14)
100 call ax_helm_lx14(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
101 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
102 coef%G23, msh%nelv)
103 case (13)
104 call ax_helm_lx13(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
105 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
106 coef%G23, msh%nelv)
107 case (12)
108 call ax_helm_lx12(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
109 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
110 coef%G23, msh%nelv)
111 case (11)
112 call ax_helm_lx11(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
113 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
114 coef%G23, msh%nelv)
115 case (10)
116 call ax_helm_lx10(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
117 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
118 coef%G23, msh%nelv)
119 case (9)
120 call ax_helm_lx9(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
121 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
122 coef%G23, msh%nelv)
123 case (8)
124 call ax_helm_lx8(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
125 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
126 coef%G23, msh%nelv)
127 case (7)
128 call ax_helm_lx7(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
129 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
130 coef%G23, msh%nelv)
131 case (6)
132 call ax_helm_lx6(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
133 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
134 coef%G23, msh%nelv)
135 case (5)
136 call ax_helm_lx5(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
137 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
138 coef%G23, msh%nelv)
139 case (4)
140 call ax_helm_lx4(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
141 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
142 coef%G23, msh%nelv)
143 case (3)
144 call ax_helm_lx3(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
145 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
146 coef%G23, msh%nelv)
147 case (2)
148 call ax_helm_lx2(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
149 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
150 coef%G23, msh%nelv)
151 case default
152 call ax_helm_lx(w, u, xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
153 coef%h1, coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, &
154 coef%G23, msh%nelv, xh%lx)
155 end select
156
157 if (coef%ifh2) then
158 !$omp do private(i)
159 do i = 1, coef%dof%size()
160 w(i,1,1,1) = w(i,1,1,1) + &
161 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * u(i,1,1,1)
162 end do
163 !$omp end do
164 end if
165 !$omp end parallel
166
167 end subroutine ax_helm_compute
168
181 subroutine ax_helm_lx(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
182 h1, G11, G22, G33, G12, G13, G23, n, lx)
183 integer, intent(in) :: n, lx
184 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
185 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
186 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
187 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
188 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
189 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
190 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
191 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
192 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
193 real(kind=rp), intent(in) :: dx(lx, lx)
194 real(kind=rp), intent(in) :: dy(lx, lx)
195 real(kind=rp), intent(in) :: dz(lx, lx)
196 real(kind=rp), intent(in) :: dxt(lx, lx)
197 real(kind=rp), intent(in) :: dyt(lx, lx)
198 real(kind=rp), intent(in) :: dzt(lx, lx)
199 real(kind=rp) :: ur(lx, lx, lx)
200 real(kind=rp) :: us(lx, lx, lx)
201 real(kind=rp) :: ut(lx, lx, lx)
202 real(kind=rp) :: wur(lx, lx, lx)
203 real(kind=rp) :: wus(lx, lx, lx)
204 real(kind=rp) :: wut(lx, lx, lx)
205 real(kind=rp) :: tmp
206 integer :: e, i, j, k, l
207
208 !$omp do
209 do e = 1, n
210 do j = 1, lx * lx
211 do i = 1, lx
212 tmp = 0.0_rp
213 do k = 1, lx
214 tmp = tmp + dx(i,k) * u(k,j,1,e)
215 end do
216 wur(i,j,1) = tmp
217 end do
218 end do
219
220 do k = 1, lx
221 do j = 1, lx
222 do i = 1, lx
223 tmp = 0.0_rp
224 do l = 1, lx
225 tmp = tmp + dy(j,l) * u(i,l,k,e)
226 end do
227 wus(i,j,k) = tmp
228 end do
229 end do
230 end do
231
232 do k = 1, lx
233 do i = 1, lx*lx
234 tmp = 0.0_rp
235 do l = 1, lx
236 tmp = tmp + dz(k,l) * u(i,1,l,e)
237 end do
238 wut(i,1,k) = tmp
239 end do
240 end do
241
242 do i = 1, lx*lx*lx
243 ur(i,1,1) = h1(i,1,1,e) &
244 * ( g11(i,1,1,e) * wur(i,1,1) &
245 + g12(i,1,1,e) * wus(i,1,1) &
246 + g13(i,1,1,e) * wut(i,1,1) )
247 us(i,1,1) = h1(i,1,1,e) &
248 * ( g12(i,1,1,e) * wur(i,1,1) &
249 + g22(i,1,1,e) * wus(i,1,1) &
250 + g23(i,1,1,e) * wut(i,1,1) )
251 ut(i,1,1) = h1(i,1,1,e) &
252 * ( g13(i,1,1,e) * wur(i,1,1) &
253 + g23(i,1,1,e) * wus(i,1,1) &
254 + g33(i,1,1,e) * wut(i,1,1) )
255 end do
256
257 do j = 1, lx*lx
258 do i = 1, lx
259 tmp = 0.0_rp
260 do k = 1, lx
261 tmp = tmp + dxt(i,k) * ur(k,j,1)
262 end do
263 w(i,j,1,e) = tmp
264 end do
265 end do
266
267 do k = 1, lx
268 do j = 1, lx
269 do i = 1, lx
270 tmp = 0.0_rp
271 do l = 1, lx
272 tmp = tmp + dyt(j,l) * us(i,l,k)
273 end do
274 w(i,j,k,e) = w(i,j,k,e) + tmp
275 end do
276 end do
277 end do
278
279 do k = 1, lx
280 do i = 1, lx*lx
281 tmp = 0.0_rp
282 do l = 1, lx
283 tmp = tmp + dzt(k,l) * ut(i,1,l)
284 end do
285 w(i,1,k,e) = w(i,1,k,e) + tmp
286 end do
287 end do
288
289 end do
290 !$omp end do
291 end subroutine ax_helm_lx
292
293 subroutine ax_helm_lx14(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
294 h1, G11, G22, G33, G12, G13, G23, n)
295 integer, parameter :: lx = 14
296 integer, intent(in) :: n
297 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
298 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
299 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
300 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
301 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
302 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
303 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
304 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
305 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
306 real(kind=rp), intent(in) :: dx(lx, lx)
307 real(kind=rp), intent(in) :: dy(lx, lx)
308 real(kind=rp), intent(in) :: dz(lx, lx)
309 real(kind=rp), intent(in) :: dxt(lx, lx)
310 real(kind=rp), intent(in) :: dyt(lx, lx)
311 real(kind=rp), intent(in) :: dzt(lx, lx)
312 real(kind=rp) :: ur(lx, lx, lx)
313 real(kind=rp) :: us(lx, lx, lx)
314 real(kind=rp) :: ut(lx, lx, lx)
315 real(kind=rp) :: wur(lx, lx, lx)
316 real(kind=rp) :: wus(lx, lx, lx)
317 real(kind=rp) :: wut(lx, lx, lx)
318 integer :: e, i, j, k
319
320 !$omp do
321 do e = 1, n
322 do j = 1, lx * lx
323 do i = 1, lx
324 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
325 + dx(i,2) * u(2,j,1,e) &
326 + dx(i,3) * u(3,j,1,e) &
327 + dx(i,4) * u(4,j,1,e) &
328 + dx(i,5) * u(5,j,1,e) &
329 + dx(i,6) * u(6,j,1,e) &
330 + dx(i,7) * u(7,j,1,e) &
331 + dx(i,8) * u(8,j,1,e) &
332 + dx(i,9) * u(9,j,1,e) &
333 + dx(i,10) * u(10,j,1,e) &
334 + dx(i,11) * u(11,j,1,e) &
335 + dx(i,12) * u(12,j,1,e) &
336 + dx(i,13) * u(13,j,1,e) &
337 + dx(i,14) * u(14,j,1,e)
338 end do
339 end do
340
341 do k = 1, lx
342 do j = 1, lx
343 do i = 1, lx
344 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
345 + dy(j,2) * u(i,2,k,e) &
346 + dy(j,3) * u(i,3,k,e) &
347 + dy(j,4) * u(i,4,k,e) &
348 + dy(j,5) * u(i,5,k,e) &
349 + dy(j,6) * u(i,6,k,e) &
350 + dy(j,7) * u(i,7,k,e) &
351 + dy(j,8) * u(i,8,k,e) &
352 + dy(j,9) * u(i,9,k,e) &
353 + dy(j,10) * u(i,10,k,e) &
354 + dy(j,11) * u(i,11,k,e) &
355 + dy(j,12) * u(i,12,k,e) &
356 + dy(j,13) * u(i,13,k,e) &
357 + dy(j,14) * u(i,14,k,e)
358 end do
359 end do
360 end do
361
362 do k = 1, lx
363 do i = 1, lx*lx
364 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
365 + dz(k,2) * u(i,1,2,e) &
366 + dz(k,3) * u(i,1,3,e) &
367 + dz(k,4) * u(i,1,4,e) &
368 + dz(k,5) * u(i,1,5,e) &
369 + dz(k,6) * u(i,1,6,e) &
370 + dz(k,7) * u(i,1,7,e) &
371 + dz(k,8) * u(i,1,8,e) &
372 + dz(k,9) * u(i,1,9,e) &
373 + dz(k,10) * u(i,1,10,e) &
374 + dz(k,11) * u(i,1,11,e) &
375 + dz(k,12) * u(i,1,12,e) &
376 + dz(k,13) * u(i,1,13,e) &
377 + dz(k,14) * u(i,1,14,e)
378 end do
379 end do
380
381 do i = 1, lx*lx*lx
382 ur(i,1,1) = h1(i,1,1,e) &
383 * ( g11(i,1,1,e) * wur(i,1,1) &
384 + g12(i,1,1,e) * wus(i,1,1) &
385 + g13(i,1,1,e) * wut(i,1,1) )
386 us(i,1,1) = h1(i,1,1,e) &
387 * ( g12(i,1,1,e) * wur(i,1,1) &
388 + g22(i,1,1,e) * wus(i,1,1) &
389 + g23(i,1,1,e) * wut(i,1,1) )
390 ut(i,1,1) = h1(i,1,1,e) &
391 * ( g13(i,1,1,e) * wur(i,1,1) &
392 + g23(i,1,1,e) * wus(i,1,1) &
393 + g33(i,1,1,e) * wut(i,1,1) )
394 end do
395
396 do j = 1, lx*lx
397 do i = 1, lx
398 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
399 + dxt(i,2) * ur(2,j,1) &
400 + dxt(i,3) * ur(3,j,1) &
401 + dxt(i,4) * ur(4,j,1) &
402 + dxt(i,5) * ur(5,j,1) &
403 + dxt(i,6) * ur(6,j,1) &
404 + dxt(i,7) * ur(7,j,1) &
405 + dxt(i,8) * ur(8,j,1) &
406 + dxt(i,9) * ur(9,j,1) &
407 + dxt(i,10) * ur(10,j,1) &
408 + dxt(i,11) * ur(11,j,1) &
409 + dxt(i,12) * ur(12,j,1) &
410 + dxt(i,13) * ur(13,j,1) &
411 + dxt(i,14) * ur(14,j,1)
412 end do
413 end do
414
415 do k = 1, lx
416 do j = 1, lx
417 do i = 1, lx
418 w(i,j,k,e) = w(i,j,k,e) &
419 + dyt(j,1) * us(i,1,k) &
420 + dyt(j,2) * us(i,2,k) &
421 + dyt(j,3) * us(i,3,k) &
422 + dyt(j,4) * us(i,4,k) &
423 + dyt(j,5) * us(i,5,k) &
424 + dyt(j,6) * us(i,6,k) &
425 + dyt(j,7) * us(i,7,k) &
426 + dyt(j,8) * us(i,8,k) &
427 + dyt(j,9) * us(i,9,k) &
428 + dyt(j,10) * us(i,10,k) &
429 + dyt(j,11) * us(i,11,k) &
430 + dyt(j,12) * us(i,12,k) &
431 + dyt(j,13) * us(i,13,k) &
432 + dyt(j,14) * us(i,14,k)
433 end do
434 end do
435 end do
436
437 do k = 1, lx
438 do i = 1, lx*lx
439 w(i,1,k,e) = w(i,1,k,e) &
440 + dzt(k,1) * ut(i,1,1) &
441 + dzt(k,2) * ut(i,1,2) &
442 + dzt(k,3) * ut(i,1,3) &
443 + dzt(k,4) * ut(i,1,4) &
444 + dzt(k,5) * ut(i,1,5) &
445 + dzt(k,6) * ut(i,1,6) &
446 + dzt(k,7) * ut(i,1,7) &
447 + dzt(k,8) * ut(i,1,8) &
448 + dzt(k,9) * ut(i,1,9) &
449 + dzt(k,10) * ut(i,1,10) &
450 + dzt(k,11) * ut(i,1,11) &
451 + dzt(k,12) * ut(i,1,12) &
452 + dzt(k,13) * ut(i,1,13) &
453 + dzt(k,14) * ut(i,1,14)
454 end do
455 end do
456
457 end do
458 !$omp end do
459 end subroutine ax_helm_lx14
460
461 subroutine ax_helm_lx13(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
462 h1, G11, G22, G33, G12, G13, G23, n)
463 integer, parameter :: lx = 13
464 integer, intent(in) :: n
465 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
466 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
467 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
468 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
469 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
470 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
471 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
472 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
473 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
474 real(kind=rp), intent(in) :: dx(lx, lx)
475 real(kind=rp), intent(in) :: dy(lx, lx)
476 real(kind=rp), intent(in) :: dz(lx, lx)
477 real(kind=rp), intent(in) :: dxt(lx, lx)
478 real(kind=rp), intent(in) :: dyt(lx, lx)
479 real(kind=rp), intent(in) :: dzt(lx, lx)
480 real(kind=rp) :: ur(lx, lx, lx)
481 real(kind=rp) :: us(lx, lx, lx)
482 real(kind=rp) :: ut(lx, lx, lx)
483 real(kind=rp) :: wur(lx, lx, lx)
484 real(kind=rp) :: wus(lx, lx, lx)
485 real(kind=rp) :: wut(lx, lx, lx)
486 integer :: e, i, j, k
487
488 !$omp do
489 do e = 1, n
490 do j = 1, lx * lx
491 do i = 1, lx
492 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
493 + dx(i,2) * u(2,j,1,e) &
494 + dx(i,3) * u(3,j,1,e) &
495 + dx(i,4) * u(4,j,1,e) &
496 + dx(i,5) * u(5,j,1,e) &
497 + dx(i,6) * u(6,j,1,e) &
498 + dx(i,7) * u(7,j,1,e) &
499 + dx(i,8) * u(8,j,1,e) &
500 + dx(i,9) * u(9,j,1,e) &
501 + dx(i,10) * u(10,j,1,e) &
502 + dx(i,11) * u(11,j,1,e) &
503 + dx(i,12) * u(12,j,1,e) &
504 + dx(i,13) * u(13,j,1,e)
505
506 end do
507 end do
508
509 do k = 1, lx
510 do j = 1, lx
511 do i = 1, lx
512 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
513 + dy(j,2) * u(i,2,k,e) &
514 + dy(j,3) * u(i,3,k,e) &
515 + dy(j,4) * u(i,4,k,e) &
516 + dy(j,5) * u(i,5,k,e) &
517 + dy(j,6) * u(i,6,k,e) &
518 + dy(j,7) * u(i,7,k,e) &
519 + dy(j,8) * u(i,8,k,e) &
520 + dy(j,9) * u(i,9,k,e) &
521 + dy(j,10) * u(i,10,k,e) &
522 + dy(j,11) * u(i,11,k,e) &
523 + dy(j,12) * u(i,12,k,e) &
524 + dy(j,13) * u(i,13,k,e)
525 end do
526 end do
527 end do
528
529 do k = 1, lx
530 do i = 1, lx*lx
531 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
532 + dz(k,2) * u(i,1,2,e) &
533 + dz(k,3) * u(i,1,3,e) &
534 + dz(k,4) * u(i,1,4,e) &
535 + dz(k,5) * u(i,1,5,e) &
536 + dz(k,6) * u(i,1,6,e) &
537 + dz(k,7) * u(i,1,7,e) &
538 + dz(k,8) * u(i,1,8,e) &
539 + dz(k,9) * u(i,1,9,e) &
540 + dz(k,10) * u(i,1,10,e) &
541 + dz(k,11) * u(i,1,11,e) &
542 + dz(k,12) * u(i,1,12,e) &
543 + dz(k,13) * u(i,1,13,e)
544 end do
545 end do
546
547 do i = 1, lx*lx*lx
548 ur(i,1,1) = h1(i,1,1,e) &
549 * ( g11(i,1,1,e) * wur(i,1,1) &
550 + g12(i,1,1,e) * wus(i,1,1) &
551 + g13(i,1,1,e) * wut(i,1,1) )
552 us(i,1,1) = h1(i,1,1,e) &
553 * ( g12(i,1,1,e) * wur(i,1,1) &
554 + g22(i,1,1,e) * wus(i,1,1) &
555 + g23(i,1,1,e) * wut(i,1,1) )
556 ut(i,1,1) = h1(i,1,1,e) &
557 * ( g13(i,1,1,e) * wur(i,1,1) &
558 + g23(i,1,1,e) * wus(i,1,1) &
559 + g33(i,1,1,e) * wut(i,1,1) )
560 end do
561
562 do j = 1, lx*lx
563 do i = 1, lx
564 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
565 + dxt(i,2) * ur(2,j,1) &
566 + dxt(i,3) * ur(3,j,1) &
567 + dxt(i,4) * ur(4,j,1) &
568 + dxt(i,5) * ur(5,j,1) &
569 + dxt(i,6) * ur(6,j,1) &
570 + dxt(i,7) * ur(7,j,1) &
571 + dxt(i,8) * ur(8,j,1) &
572 + dxt(i,9) * ur(9,j,1) &
573 + dxt(i,10) * ur(10,j,1) &
574 + dxt(i,11) * ur(11,j,1) &
575 + dxt(i,12) * ur(12,j,1) &
576 + dxt(i,13) * ur(13,j,1)
577 end do
578 end do
579
580 do k = 1, lx
581 do j = 1, lx
582 do i = 1, lx
583 w(i,j,k,e) = w(i,j,k,e) &
584 + dyt(j,1) * us(i,1,k) &
585 + dyt(j,2) * us(i,2,k) &
586 + dyt(j,3) * us(i,3,k) &
587 + dyt(j,4) * us(i,4,k) &
588 + dyt(j,5) * us(i,5,k) &
589 + dyt(j,6) * us(i,6,k) &
590 + dyt(j,7) * us(i,7,k) &
591 + dyt(j,8) * us(i,8,k) &
592 + dyt(j,9) * us(i,9,k) &
593 + dyt(j,10) * us(i,10,k) &
594 + dyt(j,11) * us(i,11,k) &
595 + dyt(j,12) * us(i,12,k) &
596 + dyt(j,13) * us(i,13,k)
597 end do
598 end do
599 end do
600
601 do k = 1, lx
602 do i = 1, lx*lx
603 w(i,1,k,e) = w(i,1,k,e) &
604 + dzt(k,1) * ut(i,1,1) &
605 + dzt(k,2) * ut(i,1,2) &
606 + dzt(k,3) * ut(i,1,3) &
607 + dzt(k,4) * ut(i,1,4) &
608 + dzt(k,5) * ut(i,1,5) &
609 + dzt(k,6) * ut(i,1,6) &
610 + dzt(k,7) * ut(i,1,7) &
611 + dzt(k,8) * ut(i,1,8) &
612 + dzt(k,9) * ut(i,1,9) &
613 + dzt(k,10) * ut(i,1,10) &
614 + dzt(k,11) * ut(i,1,11) &
615 + dzt(k,12) * ut(i,1,12) &
616 + dzt(k,13) * ut(i,1,13)
617 end do
618 end do
619
620 end do
621 !$omp end do
622 end subroutine ax_helm_lx13
623
624 subroutine ax_helm_lx12(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
625 h1, G11, G22, G33, G12, G13, G23, n)
626 integer, parameter :: lx = 12
627 integer, intent(in) :: n
628 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
629 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
630 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
631 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
632 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
633 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
634 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
635 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
636 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
637 real(kind=rp), intent(in) :: dx(lx, lx)
638 real(kind=rp), intent(in) :: dy(lx, lx)
639 real(kind=rp), intent(in) :: dz(lx, lx)
640 real(kind=rp), intent(in) :: dxt(lx, lx)
641 real(kind=rp), intent(in) :: dyt(lx, lx)
642 real(kind=rp), intent(in) :: dzt(lx, lx)
643 real(kind=rp) :: ur(lx, lx, lx)
644 real(kind=rp) :: us(lx, lx, lx)
645 real(kind=rp) :: ut(lx, lx, lx)
646 real(kind=rp) :: wur(lx, lx, lx)
647 real(kind=rp) :: wus(lx, lx, lx)
648 real(kind=rp) :: wut(lx, lx, lx)
649 integer :: e, i, j, k
650
651 !$omp do
652 do e = 1, n
653 do j = 1, lx * lx
654 do i = 1, lx
655 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
656 + dx(i,2) * u(2,j,1,e) &
657 + dx(i,3) * u(3,j,1,e) &
658 + dx(i,4) * u(4,j,1,e) &
659 + dx(i,5) * u(5,j,1,e) &
660 + dx(i,6) * u(6,j,1,e) &
661 + dx(i,7) * u(7,j,1,e) &
662 + dx(i,8) * u(8,j,1,e) &
663 + dx(i,9) * u(9,j,1,e) &
664 + dx(i,10) * u(10,j,1,e) &
665 + dx(i,11) * u(11,j,1,e) &
666 + dx(i,12) * u(12,j,1,e)
667 end do
668 end do
669
670 do k = 1, lx
671 do j = 1, lx
672 do i = 1, lx
673 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
674 + dy(j,2) * u(i,2,k,e) &
675 + dy(j,3) * u(i,3,k,e) &
676 + dy(j,4) * u(i,4,k,e) &
677 + dy(j,5) * u(i,5,k,e) &
678 + dy(j,6) * u(i,6,k,e) &
679 + dy(j,7) * u(i,7,k,e) &
680 + dy(j,8) * u(i,8,k,e) &
681 + dy(j,9) * u(i,9,k,e) &
682 + dy(j,10) * u(i,10,k,e) &
683 + dy(j,11) * u(i,11,k,e) &
684 + dy(j,12) * u(i,12,k,e)
685 end do
686 end do
687 end do
688
689 do k = 1, lx
690 do i = 1, lx*lx
691 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
692 + dz(k,2) * u(i,1,2,e) &
693 + dz(k,3) * u(i,1,3,e) &
694 + dz(k,4) * u(i,1,4,e) &
695 + dz(k,5) * u(i,1,5,e) &
696 + dz(k,6) * u(i,1,6,e) &
697 + dz(k,7) * u(i,1,7,e) &
698 + dz(k,8) * u(i,1,8,e) &
699 + dz(k,9) * u(i,1,9,e) &
700 + dz(k,10) * u(i,1,10,e) &
701 + dz(k,11) * u(i,1,11,e) &
702 + dz(k,12) * u(i,1,12,e)
703 end do
704 end do
705
706 do i = 1, lx*lx*lx
707 ur(i,1,1) = h1(i,1,1,e) &
708 * ( g11(i,1,1,e) * wur(i,1,1) &
709 + g12(i,1,1,e) * wus(i,1,1) &
710 + g13(i,1,1,e) * wut(i,1,1) )
711 us(i,1,1) = h1(i,1,1,e) &
712 * ( g12(i,1,1,e) * wur(i,1,1) &
713 + g22(i,1,1,e) * wus(i,1,1) &
714 + g23(i,1,1,e) * wut(i,1,1) )
715 ut(i,1,1) = h1(i,1,1,e) &
716 * ( g13(i,1,1,e) * wur(i,1,1) &
717 + g23(i,1,1,e) * wus(i,1,1) &
718 + g33(i,1,1,e) * wut(i,1,1) )
719 end do
720
721 do j = 1, lx*lx
722 do i = 1, lx
723 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
724 + dxt(i,2) * ur(2,j,1) &
725 + dxt(i,3) * ur(3,j,1) &
726 + dxt(i,4) * ur(4,j,1) &
727 + dxt(i,5) * ur(5,j,1) &
728 + dxt(i,6) * ur(6,j,1) &
729 + dxt(i,7) * ur(7,j,1) &
730 + dxt(i,8) * ur(8,j,1) &
731 + dxt(i,9) * ur(9,j,1) &
732 + dxt(i,10) * ur(10,j,1) &
733 + dxt(i,11) * ur(11,j,1) &
734 + dxt(i,12) * ur(12,j,1)
735 end do
736 end do
737
738 do k = 1, lx
739 do j = 1, lx
740 do i = 1, lx
741 w(i,j,k,e) = w(i,j,k,e) &
742 + dyt(j,1) * us(i,1,k) &
743 + dyt(j,2) * us(i,2,k) &
744 + dyt(j,3) * us(i,3,k) &
745 + dyt(j,4) * us(i,4,k) &
746 + dyt(j,5) * us(i,5,k) &
747 + dyt(j,6) * us(i,6,k) &
748 + dyt(j,7) * us(i,7,k) &
749 + dyt(j,8) * us(i,8,k) &
750 + dyt(j,9) * us(i,9,k) &
751 + dyt(j,10) * us(i,10,k) &
752 + dyt(j,11) * us(i,11,k) &
753 + dyt(j,12) * us(i,12,k)
754 end do
755 end do
756 end do
757
758 do k = 1, lx
759 do i = 1, lx*lx
760 w(i,1,k,e) = w(i,1,k,e) &
761 + dzt(k,1) * ut(i,1,1) &
762 + dzt(k,2) * ut(i,1,2) &
763 + dzt(k,3) * ut(i,1,3) &
764 + dzt(k,4) * ut(i,1,4) &
765 + dzt(k,5) * ut(i,1,5) &
766 + dzt(k,6) * ut(i,1,6) &
767 + dzt(k,7) * ut(i,1,7) &
768 + dzt(k,8) * ut(i,1,8) &
769 + dzt(k,9) * ut(i,1,9) &
770 + dzt(k,10) * ut(i,1,10) &
771 + dzt(k,11) * ut(i,1,11) &
772 + dzt(k,12) * ut(i,1,12)
773 end do
774 end do
775
776 end do
777 !$omp end do
778 end subroutine ax_helm_lx12
779
780 subroutine ax_helm_lx11(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
781 h1, G11, G22, G33, G12, G13, G23, n)
782 integer, parameter :: lx = 11
783 integer, intent(in) :: n
784 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
785 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
786 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
787 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
788 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
789 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
790 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
791 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
792 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
793 real(kind=rp), intent(in) :: dx(lx, lx)
794 real(kind=rp), intent(in) :: dy(lx, lx)
795 real(kind=rp), intent(in) :: dz(lx, lx)
796 real(kind=rp), intent(in) :: dxt(lx, lx)
797 real(kind=rp), intent(in) :: dyt(lx, lx)
798 real(kind=rp), intent(in) :: dzt(lx, lx)
799 real(kind=rp) :: ur(lx, lx, lx)
800 real(kind=rp) :: us(lx, lx, lx)
801 real(kind=rp) :: ut(lx, lx, lx)
802 real(kind=rp) :: wur(lx, lx, lx)
803 real(kind=rp) :: wus(lx, lx, lx)
804 real(kind=rp) :: wut(lx, lx, lx)
805 integer :: e, i, j, k
806
807 !$omp do
808 do e = 1, n
809 do j = 1, lx * lx
810 do i = 1, lx
811 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
812 + dx(i,2) * u(2,j,1,e) &
813 + dx(i,3) * u(3,j,1,e) &
814 + dx(i,4) * u(4,j,1,e) &
815 + dx(i,5) * u(5,j,1,e) &
816 + dx(i,6) * u(6,j,1,e) &
817 + dx(i,7) * u(7,j,1,e) &
818 + dx(i,8) * u(8,j,1,e) &
819 + dx(i,9) * u(9,j,1,e) &
820 + dx(i,10) * u(10,j,1,e) &
821 + dx(i,11) * u(11,j,1,e)
822 end do
823 end do
824
825 do k = 1, lx
826 do j = 1, lx
827 do i = 1, lx
828 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
829 + dy(j,2) * u(i,2,k,e) &
830 + dy(j,3) * u(i,3,k,e) &
831 + dy(j,4) * u(i,4,k,e) &
832 + dy(j,5) * u(i,5,k,e) &
833 + dy(j,6) * u(i,6,k,e) &
834 + dy(j,7) * u(i,7,k,e) &
835 + dy(j,8) * u(i,8,k,e) &
836 + dy(j,9) * u(i,9,k,e) &
837 + dy(j,10) * u(i,10,k,e) &
838 + dy(j,11) * u(i,11,k,e)
839 end do
840 end do
841 end do
842
843 do k = 1, lx
844 do i = 1, lx*lx
845 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
846 + dz(k,2) * u(i,1,2,e) &
847 + dz(k,3) * u(i,1,3,e) &
848 + dz(k,4) * u(i,1,4,e) &
849 + dz(k,5) * u(i,1,5,e) &
850 + dz(k,6) * u(i,1,6,e) &
851 + dz(k,7) * u(i,1,7,e) &
852 + dz(k,8) * u(i,1,8,e) &
853 + dz(k,9) * u(i,1,9,e) &
854 + dz(k,10) * u(i,1,10,e) &
855 + dz(k,11) * u(i,1,11,e)
856 end do
857 end do
858
859 do i = 1, lx*lx*lx
860 ur(i,1,1) = h1(i,1,1,e) &
861 * ( g11(i,1,1,e) * wur(i,1,1) &
862 + g12(i,1,1,e) * wus(i,1,1) &
863 + g13(i,1,1,e) * wut(i,1,1) )
864 us(i,1,1) = h1(i,1,1,e) &
865 * ( g12(i,1,1,e) * wur(i,1,1) &
866 + g22(i,1,1,e) * wus(i,1,1) &
867 + g23(i,1,1,e) * wut(i,1,1) )
868 ut(i,1,1) = h1(i,1,1,e) &
869 * ( g13(i,1,1,e) * wur(i,1,1) &
870 + g23(i,1,1,e) * wus(i,1,1) &
871 + g33(i,1,1,e) * wut(i,1,1) )
872 end do
873
874 do j = 1, lx*lx
875 do i = 1, lx
876 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
877 + dxt(i,2) * ur(2,j,1) &
878 + dxt(i,3) * ur(3,j,1) &
879 + dxt(i,4) * ur(4,j,1) &
880 + dxt(i,5) * ur(5,j,1) &
881 + dxt(i,6) * ur(6,j,1) &
882 + dxt(i,7) * ur(7,j,1) &
883 + dxt(i,8) * ur(8,j,1) &
884 + dxt(i,9) * ur(9,j,1) &
885 + dxt(i,10) * ur(10,j,1) &
886 + dxt(i,11) * ur(11,j,1)
887 end do
888 end do
889
890 do k = 1, lx
891 do j = 1, lx
892 do i = 1, lx
893 w(i,j,k,e) = w(i,j,k,e) &
894 + dyt(j,1) * us(i,1,k) &
895 + dyt(j,2) * us(i,2,k) &
896 + dyt(j,3) * us(i,3,k) &
897 + dyt(j,4) * us(i,4,k) &
898 + dyt(j,5) * us(i,5,k) &
899 + dyt(j,6) * us(i,6,k) &
900 + dyt(j,7) * us(i,7,k) &
901 + dyt(j,8) * us(i,8,k) &
902 + dyt(j,9) * us(i,9,k) &
903 + dyt(j,10) * us(i,10,k) &
904 + dyt(j,11) * us(i,11,k)
905 end do
906 end do
907 end do
908
909 do k = 1, lx
910 do i = 1, lx*lx
911 w(i,1,k,e) = w(i,1,k,e) &
912 + dzt(k,1) * ut(i,1,1) &
913 + dzt(k,2) * ut(i,1,2) &
914 + dzt(k,3) * ut(i,1,3) &
915 + dzt(k,4) * ut(i,1,4) &
916 + dzt(k,5) * ut(i,1,5) &
917 + dzt(k,6) * ut(i,1,6) &
918 + dzt(k,7) * ut(i,1,7) &
919 + dzt(k,8) * ut(i,1,8) &
920 + dzt(k,9) * ut(i,1,9) &
921 + dzt(k,10) * ut(i,1,10) &
922 + dzt(k,11) * ut(i,1,11)
923 end do
924 end do
925
926 end do
927 !$omp end do
928 end subroutine ax_helm_lx11
929
930 subroutine ax_helm_lx10(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
931 h1, G11, G22, G33, G12, G13, G23, n)
932 integer, parameter :: lx = 10
933 integer, intent(in) :: n
934 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
935 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
936 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
937 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
938 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
939 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
940 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
941 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
942 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
943 real(kind=rp), intent(in) :: dx(lx, lx)
944 real(kind=rp), intent(in) :: dy(lx, lx)
945 real(kind=rp), intent(in) :: dz(lx, lx)
946 real(kind=rp), intent(in) :: dxt(lx, lx)
947 real(kind=rp), intent(in) :: dyt(lx, lx)
948 real(kind=rp), intent(in) :: dzt(lx, lx)
949 real(kind=rp) :: ur(lx, lx, lx)
950 real(kind=rp) :: us(lx, lx, lx)
951 real(kind=rp) :: ut(lx, lx, lx)
952 real(kind=rp) :: wur(lx, lx, lx)
953 real(kind=rp) :: wus(lx, lx, lx)
954 real(kind=rp) :: wut(lx, lx, lx)
955 integer :: e, i, j, k
956
957 !$omp do
958 do e = 1, n
959 do j = 1, lx * lx
960 do i = 1, lx
961 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
962 + dx(i,2) * u(2,j,1,e) &
963 + dx(i,3) * u(3,j,1,e) &
964 + dx(i,4) * u(4,j,1,e) &
965 + dx(i,5) * u(5,j,1,e) &
966 + dx(i,6) * u(6,j,1,e) &
967 + dx(i,7) * u(7,j,1,e) &
968 + dx(i,8) * u(8,j,1,e) &
969 + dx(i,9) * u(9,j,1,e) &
970 + dx(i,10) * u(10,j,1,e)
971 end do
972 end do
973
974 do k = 1, lx
975 do j = 1, lx
976 do i = 1, lx
977 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
978 + dy(j,2) * u(i,2,k,e) &
979 + dy(j,3) * u(i,3,k,e) &
980 + dy(j,4) * u(i,4,k,e) &
981 + dy(j,5) * u(i,5,k,e) &
982 + dy(j,6) * u(i,6,k,e) &
983 + dy(j,7) * u(i,7,k,e) &
984 + dy(j,8) * u(i,8,k,e) &
985 + dy(j,9) * u(i,9,k,e) &
986 + dy(j,10) * u(i,10,k,e)
987 end do
988 end do
989 end do
990
991 do k = 1, lx
992 do i = 1, lx*lx
993 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
994 + dz(k,2) * u(i,1,2,e) &
995 + dz(k,3) * u(i,1,3,e) &
996 + dz(k,4) * u(i,1,4,e) &
997 + dz(k,5) * u(i,1,5,e) &
998 + dz(k,6) * u(i,1,6,e) &
999 + dz(k,7) * u(i,1,7,e) &
1000 + dz(k,8) * u(i,1,8,e) &
1001 + dz(k,9) * u(i,1,9,e) &
1002 + dz(k,10) * u(i,1,10,e)
1003 end do
1004 end do
1005
1006 do i = 1, lx*lx*lx
1007 ur(i,1,1) = h1(i,1,1,e) &
1008 * ( g11(i,1,1,e) * wur(i,1,1) &
1009 + g12(i,1,1,e) * wus(i,1,1) &
1010 + g13(i,1,1,e) * wut(i,1,1) )
1011 us(i,1,1) = h1(i,1,1,e) &
1012 * ( g12(i,1,1,e) * wur(i,1,1) &
1013 + g22(i,1,1,e) * wus(i,1,1) &
1014 + g23(i,1,1,e) * wut(i,1,1) )
1015 ut(i,1,1) = h1(i,1,1,e) &
1016 * ( g13(i,1,1,e) * wur(i,1,1) &
1017 + g23(i,1,1,e) * wus(i,1,1) &
1018 + g33(i,1,1,e) * wut(i,1,1) )
1019 end do
1020
1021 do j = 1, lx*lx
1022 do i = 1, lx
1023 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1024 + dxt(i,2) * ur(2,j,1) &
1025 + dxt(i,3) * ur(3,j,1) &
1026 + dxt(i,4) * ur(4,j,1) &
1027 + dxt(i,5) * ur(5,j,1) &
1028 + dxt(i,6) * ur(6,j,1) &
1029 + dxt(i,7) * ur(7,j,1) &
1030 + dxt(i,8) * ur(8,j,1) &
1031 + dxt(i,9) * ur(9,j,1) &
1032 + dxt(i,10) * ur(10,j,1)
1033 end do
1034 end do
1035
1036 do k = 1, lx
1037 do j = 1, lx
1038 do i = 1, lx
1039 w(i,j,k,e) = w(i,j,k,e) &
1040 + dyt(j,1) * us(i,1,k) &
1041 + dyt(j,2) * us(i,2,k) &
1042 + dyt(j,3) * us(i,3,k) &
1043 + dyt(j,4) * us(i,4,k) &
1044 + dyt(j,5) * us(i,5,k) &
1045 + dyt(j,6) * us(i,6,k) &
1046 + dyt(j,7) * us(i,7,k) &
1047 + dyt(j,8) * us(i,8,k) &
1048 + dyt(j,9) * us(i,9,k) &
1049 + dyt(j,10) * us(i,10,k)
1050 end do
1051 end do
1052 end do
1053
1054 do k = 1, lx
1055 do i = 1, lx*lx
1056 w(i,1,k,e) = w(i,1,k,e) &
1057 + dzt(k,1) * ut(i,1,1) &
1058 + dzt(k,2) * ut(i,1,2) &
1059 + dzt(k,3) * ut(i,1,3) &
1060 + dzt(k,4) * ut(i,1,4) &
1061 + dzt(k,5) * ut(i,1,5) &
1062 + dzt(k,6) * ut(i,1,6) &
1063 + dzt(k,7) * ut(i,1,7) &
1064 + dzt(k,8) * ut(i,1,8) &
1065 + dzt(k,9) * ut(i,1,9) &
1066 + dzt(k,10) * ut(i,1,10)
1067 end do
1068 end do
1069
1070 end do
1071 !$omp end do
1072 end subroutine ax_helm_lx10
1073
1074 subroutine ax_helm_lx9(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1075 h1, G11, G22, G33, G12, G13, G23, n)
1076 integer, parameter :: lx = 9
1077 integer, intent(in) :: n
1078 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1079 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1080 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1081 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1082 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1083 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1084 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1085 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1086 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1087 real(kind=rp), intent(in) :: dx(lx, lx)
1088 real(kind=rp), intent(in) :: dy(lx, lx)
1089 real(kind=rp), intent(in) :: dz(lx, lx)
1090 real(kind=rp), intent(in) :: dxt(lx, lx)
1091 real(kind=rp), intent(in) :: dyt(lx, lx)
1092 real(kind=rp), intent(in) :: dzt(lx, lx)
1093 real(kind=rp) :: ur(lx, lx, lx)
1094 real(kind=rp) :: us(lx, lx, lx)
1095 real(kind=rp) :: ut(lx, lx, lx)
1096 real(kind=rp) :: wur(lx, lx, lx)
1097 real(kind=rp) :: wus(lx, lx, lx)
1098 real(kind=rp) :: wut(lx, lx, lx)
1099 integer :: e, i, j, k
1100
1101 !$omp do
1102 do e = 1, n
1103 do j = 1, lx * lx
1104 do i = 1, lx
1105 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1106 + dx(i,2) * u(2,j,1,e) &
1107 + dx(i,3) * u(3,j,1,e) &
1108 + dx(i,4) * u(4,j,1,e) &
1109 + dx(i,5) * u(5,j,1,e) &
1110 + dx(i,6) * u(6,j,1,e) &
1111 + dx(i,7) * u(7,j,1,e) &
1112 + dx(i,8) * u(8,j,1,e) &
1113 + dx(i,9) * u(9,j,1,e)
1114 end do
1115 end do
1116
1117 do k = 1, lx
1118 do j = 1, lx
1119 do i = 1, lx
1120 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1121 + dy(j,2) * u(i,2,k,e) &
1122 + dy(j,3) * u(i,3,k,e) &
1123 + dy(j,4) * u(i,4,k,e) &
1124 + dy(j,5) * u(i,5,k,e) &
1125 + dy(j,6) * u(i,6,k,e) &
1126 + dy(j,7) * u(i,7,k,e) &
1127 + dy(j,8) * u(i,8,k,e) &
1128 + dy(j,9) * u(i,9,k,e)
1129 end do
1130 end do
1131 end do
1132
1133 do k = 1, lx
1134 do i = 1, lx*lx
1135 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1136 + dz(k,2) * u(i,1,2,e) &
1137 + dz(k,3) * u(i,1,3,e) &
1138 + dz(k,4) * u(i,1,4,e) &
1139 + dz(k,5) * u(i,1,5,e) &
1140 + dz(k,6) * u(i,1,6,e) &
1141 + dz(k,7) * u(i,1,7,e) &
1142 + dz(k,8) * u(i,1,8,e) &
1143 + dz(k,9) * u(i,1,9,e)
1144 end do
1145 end do
1146
1147 do i = 1, lx*lx*lx
1148 ur(i,1,1) = h1(i,1,1,e) &
1149 * ( g11(i,1,1,e) * wur(i,1,1) &
1150 + g12(i,1,1,e) * wus(i,1,1) &
1151 + g13(i,1,1,e) * wut(i,1,1) )
1152 us(i,1,1) = h1(i,1,1,e) &
1153 * ( g12(i,1,1,e) * wur(i,1,1) &
1154 + g22(i,1,1,e) * wus(i,1,1) &
1155 + g23(i,1,1,e) * wut(i,1,1) )
1156 ut(i,1,1) = h1(i,1,1,e) &
1157 * ( g13(i,1,1,e) * wur(i,1,1) &
1158 + g23(i,1,1,e) * wus(i,1,1) &
1159 + g33(i,1,1,e) * wut(i,1,1) )
1160 end do
1161
1162 do j = 1, lx*lx
1163 do i = 1, lx
1164 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1165 + dxt(i,2) * ur(2,j,1) &
1166 + dxt(i,3) * ur(3,j,1) &
1167 + dxt(i,4) * ur(4,j,1) &
1168 + dxt(i,5) * ur(5,j,1) &
1169 + dxt(i,6) * ur(6,j,1) &
1170 + dxt(i,7) * ur(7,j,1) &
1171 + dxt(i,8) * ur(8,j,1) &
1172 + dxt(i,9) * ur(9,j,1)
1173 end do
1174 end do
1175
1176 do k = 1, lx
1177 do j = 1, lx
1178 do i = 1, lx
1179 w(i,j,k,e) = w(i,j,k,e) &
1180 + dyt(j,1) * us(i,1,k) &
1181 + dyt(j,2) * us(i,2,k) &
1182 + dyt(j,3) * us(i,3,k) &
1183 + dyt(j,4) * us(i,4,k) &
1184 + dyt(j,5) * us(i,5,k) &
1185 + dyt(j,6) * us(i,6,k) &
1186 + dyt(j,7) * us(i,7,k) &
1187 + dyt(j,8) * us(i,8,k) &
1188 + dyt(j,9) * us(i,9,k)
1189 end do
1190 end do
1191 end do
1192
1193 do k = 1, lx
1194 do i = 1, lx*lx
1195 w(i,1,k,e) = w(i,1,k,e) &
1196 + dzt(k,1) * ut(i,1,1) &
1197 + dzt(k,2) * ut(i,1,2) &
1198 + dzt(k,3) * ut(i,1,3) &
1199 + dzt(k,4) * ut(i,1,4) &
1200 + dzt(k,5) * ut(i,1,5) &
1201 + dzt(k,6) * ut(i,1,6) &
1202 + dzt(k,7) * ut(i,1,7) &
1203 + dzt(k,8) * ut(i,1,8) &
1204 + dzt(k,9) * ut(i,1,9)
1205 end do
1206 end do
1207
1208 end do
1209 !$omp end do
1210 end subroutine ax_helm_lx9
1211
1212 subroutine ax_helm_lx8(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1213 h1, G11, G22, G33, G12, G13, G23, n)
1214 integer, parameter :: lx = 8
1215 integer, intent(in) :: n
1216 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1217 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1218 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1219 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1220 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1221 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1222 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1223 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1224 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1225 real(kind=rp), intent(in) :: dx(lx, lx)
1226 real(kind=rp), intent(in) :: dy(lx, lx)
1227 real(kind=rp), intent(in) :: dz(lx, lx)
1228 real(kind=rp), intent(in) :: dxt(lx, lx)
1229 real(kind=rp), intent(in) :: dyt(lx, lx)
1230 real(kind=rp), intent(in) :: dzt(lx, lx)
1231 real(kind=rp) :: ur(lx, lx, lx)
1232 real(kind=rp) :: us(lx, lx, lx)
1233 real(kind=rp) :: ut(lx, lx, lx)
1234 real(kind=rp) :: wur(lx, lx, lx)
1235 real(kind=rp) :: wus(lx, lx, lx)
1236 real(kind=rp) :: wut(lx, lx, lx)
1237 integer :: e, i, j, k
1238
1239 !$omp do
1240 do e = 1, n
1241 do j = 1, lx * lx
1242 do i = 1, lx
1243 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1244 + dx(i,2) * u(2,j,1,e) &
1245 + dx(i,3) * u(3,j,1,e) &
1246 + dx(i,4) * u(4,j,1,e) &
1247 + dx(i,5) * u(5,j,1,e) &
1248 + dx(i,6) * u(6,j,1,e) &
1249 + dx(i,7) * u(7,j,1,e) &
1250 + dx(i,8) * u(8,j,1,e)
1251 end do
1252 end do
1253
1254 do k = 1, lx
1255 do j = 1, lx
1256 do i = 1, lx
1257 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1258 + dy(j,2) * u(i,2,k,e) &
1259 + dy(j,3) * u(i,3,k,e) &
1260 + dy(j,4) * u(i,4,k,e) &
1261 + dy(j,5) * u(i,5,k,e) &
1262 + dy(j,6) * u(i,6,k,e) &
1263 + dy(j,7) * u(i,7,k,e) &
1264 + dy(j,8) * u(i,8,k,e)
1265 end do
1266 end do
1267 end do
1268
1269 do k = 1, lx
1270 do i = 1, lx*lx
1271 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1272 + dz(k,2) * u(i,1,2,e) &
1273 + dz(k,3) * u(i,1,3,e) &
1274 + dz(k,4) * u(i,1,4,e) &
1275 + dz(k,5) * u(i,1,5,e) &
1276 + dz(k,6) * u(i,1,6,e) &
1277 + dz(k,7) * u(i,1,7,e) &
1278 + dz(k,8) * u(i,1,8,e)
1279 end do
1280 end do
1281
1282 do i = 1, lx*lx*lx
1283 ur(i,1,1) = h1(i,1,1,e) &
1284 * ( g11(i,1,1,e) * wur(i,1,1) &
1285 + g12(i,1,1,e) * wus(i,1,1) &
1286 + g13(i,1,1,e) * wut(i,1,1) )
1287 us(i,1,1) = h1(i,1,1,e) &
1288 * ( g12(i,1,1,e) * wur(i,1,1) &
1289 + g22(i,1,1,e) * wus(i,1,1) &
1290 + g23(i,1,1,e) * wut(i,1,1) )
1291 ut(i,1,1) = h1(i,1,1,e) &
1292 * ( g13(i,1,1,e) * wur(i,1,1) &
1293 + g23(i,1,1,e) * wus(i,1,1) &
1294 + g33(i,1,1,e) * wut(i,1,1) )
1295 end do
1296
1297 do j = 1, lx*lx
1298 do i = 1, lx
1299 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1300 + dxt(i,2) * ur(2,j,1) &
1301 + dxt(i,3) * ur(3,j,1) &
1302 + dxt(i,4) * ur(4,j,1) &
1303 + dxt(i,5) * ur(5,j,1) &
1304 + dxt(i,6) * ur(6,j,1) &
1305 + dxt(i,7) * ur(7,j,1) &
1306 + dxt(i,8) * ur(8,j,1)
1307 end do
1308 end do
1309
1310 do k = 1, lx
1311 do j = 1, lx
1312 do i = 1, lx
1313 w(i,j,k,e) = w(i,j,k,e) &
1314 + dyt(j,1) * us(i,1,k) &
1315 + dyt(j,2) * us(i,2,k) &
1316 + dyt(j,3) * us(i,3,k) &
1317 + dyt(j,4) * us(i,4,k) &
1318 + dyt(j,5) * us(i,5,k) &
1319 + dyt(j,6) * us(i,6,k) &
1320 + dyt(j,7) * us(i,7,k) &
1321 + dyt(j,8) * us(i,8,k)
1322 end do
1323 end do
1324 end do
1325
1326 do k = 1, lx
1327 do i = 1, lx*lx
1328 w(i,1,k,e) = w(i,1,k,e) &
1329 + dzt(k,1) * ut(i,1,1) &
1330 + dzt(k,2) * ut(i,1,2) &
1331 + dzt(k,3) * ut(i,1,3) &
1332 + dzt(k,4) * ut(i,1,4) &
1333 + dzt(k,5) * ut(i,1,5) &
1334 + dzt(k,6) * ut(i,1,6) &
1335 + dzt(k,7) * ut(i,1,7) &
1336 + dzt(k,8) * ut(i,1,8)
1337 end do
1338 end do
1339
1340 end do
1341 !$omp end do
1342 end subroutine ax_helm_lx8
1343
1344 subroutine ax_helm_lx7(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1345 h1, G11, G22, G33, G12, G13, G23, n)
1346 integer, parameter :: lx = 7
1347 integer, intent(in) :: n
1348 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1349 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1350 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1351 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1352 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1353 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1354 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1355 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1356 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1357 real(kind=rp), intent(in) :: dx(lx, lx)
1358 real(kind=rp), intent(in) :: dy(lx, lx)
1359 real(kind=rp), intent(in) :: dz(lx, lx)
1360 real(kind=rp), intent(in) :: dxt(lx, lx)
1361 real(kind=rp), intent(in) :: dyt(lx, lx)
1362 real(kind=rp), intent(in) :: dzt(lx, lx)
1363 real(kind=rp) :: ur(lx, lx, lx)
1364 real(kind=rp) :: us(lx, lx, lx)
1365 real(kind=rp) :: ut(lx, lx, lx)
1366 real(kind=rp) :: wur(lx, lx, lx)
1367 real(kind=rp) :: wus(lx, lx, lx)
1368 real(kind=rp) :: wut(lx, lx, lx)
1369 integer :: e, i, j, k
1370
1371 !$omp do
1372 do e = 1, n
1373 do j = 1, lx * lx
1374 do i = 1, lx
1375 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1376 + dx(i,2) * u(2,j,1,e) &
1377 + dx(i,3) * u(3,j,1,e) &
1378 + dx(i,4) * u(4,j,1,e) &
1379 + dx(i,5) * u(5,j,1,e) &
1380 + dx(i,6) * u(6,j,1,e) &
1381 + dx(i,7) * u(7,j,1,e)
1382 end do
1383 end do
1384
1385 do k = 1, lx
1386 do j = 1, lx
1387 do i = 1, lx
1388 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1389 + dy(j,2) * u(i,2,k,e) &
1390 + dy(j,3) * u(i,3,k,e) &
1391 + dy(j,4) * u(i,4,k,e) &
1392 + dy(j,5) * u(i,5,k,e) &
1393 + dy(j,6) * u(i,6,k,e) &
1394 + dy(j,7) * u(i,7,k,e)
1395 end do
1396 end do
1397 end do
1398
1399 do k = 1, lx
1400 do i = 1, lx*lx
1401 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1402 + dz(k,2) * u(i,1,2,e) &
1403 + dz(k,3) * u(i,1,3,e) &
1404 + dz(k,4) * u(i,1,4,e) &
1405 + dz(k,5) * u(i,1,5,e) &
1406 + dz(k,6) * u(i,1,6,e) &
1407 + dz(k,7) * u(i,1,7,e)
1408 end do
1409 end do
1410
1411 do i = 1, lx*lx*lx
1412 ur(i,1,1) = h1(i,1,1,e) &
1413 * ( g11(i,1,1,e) * wur(i,1,1) &
1414 + g12(i,1,1,e) * wus(i,1,1) &
1415 + g13(i,1,1,e) * wut(i,1,1) )
1416 us(i,1,1) = h1(i,1,1,e) &
1417 * ( g12(i,1,1,e) * wur(i,1,1) &
1418 + g22(i,1,1,e) * wus(i,1,1) &
1419 + g23(i,1,1,e) * wut(i,1,1) )
1420 ut(i,1,1) = h1(i,1,1,e) &
1421 * ( g13(i,1,1,e) * wur(i,1,1) &
1422 + g23(i,1,1,e) * wus(i,1,1) &
1423 + g33(i,1,1,e) * wut(i,1,1) )
1424 end do
1425
1426 do j = 1, lx*lx
1427 do i = 1, lx
1428 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1429 + dxt(i,2) * ur(2,j,1) &
1430 + dxt(i,3) * ur(3,j,1) &
1431 + dxt(i,4) * ur(4,j,1) &
1432 + dxt(i,5) * ur(5,j,1) &
1433 + dxt(i,6) * ur(6,j,1) &
1434 + dxt(i,7) * ur(7,j,1)
1435 end do
1436 end do
1437
1438 do k = 1, lx
1439 do j = 1, lx
1440 do i = 1, lx
1441 w(i,j,k,e) = w(i,j,k,e) &
1442 + dyt(j,1) * us(i,1,k) &
1443 + dyt(j,2) * us(i,2,k) &
1444 + dyt(j,3) * us(i,3,k) &
1445 + dyt(j,4) * us(i,4,k) &
1446 + dyt(j,5) * us(i,5,k) &
1447 + dyt(j,6) * us(i,6,k) &
1448 + dyt(j,7) * us(i,7,k)
1449 end do
1450 end do
1451 end do
1452
1453 do k = 1, lx
1454 do i = 1, lx*lx
1455 w(i,1,k,e) = w(i,1,k,e) &
1456 + dzt(k,1) * ut(i,1,1) &
1457 + dzt(k,2) * ut(i,1,2) &
1458 + dzt(k,3) * ut(i,1,3) &
1459 + dzt(k,4) * ut(i,1,4) &
1460 + dzt(k,5) * ut(i,1,5) &
1461 + dzt(k,6) * ut(i,1,6) &
1462 + dzt(k,7) * ut(i,1,7)
1463 end do
1464 end do
1465
1466 end do
1467 !$omp end do
1468 end subroutine ax_helm_lx7
1469
1470 subroutine ax_helm_lx6(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1471 h1, G11, G22, G33, G12, G13, G23, n)
1472 integer, parameter :: lx = 6
1473 integer, intent(in) :: n
1474 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1475 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1476 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1477 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1478 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1479 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1480 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1481 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1482 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1483 real(kind=rp), intent(in) :: dx(lx, lx)
1484 real(kind=rp), intent(in) :: dy(lx, lx)
1485 real(kind=rp), intent(in) :: dz(lx, lx)
1486 real(kind=rp), intent(in) :: dxt(lx, lx)
1487 real(kind=rp), intent(in) :: dyt(lx, lx)
1488 real(kind=rp), intent(in) :: dzt(lx, lx)
1489 real(kind=rp) :: ur(lx, lx, lx)
1490 real(kind=rp) :: us(lx, lx, lx)
1491 real(kind=rp) :: ut(lx, lx, lx)
1492 real(kind=rp) :: wur(lx, lx, lx)
1493 real(kind=rp) :: wus(lx, lx, lx)
1494 real(kind=rp) :: wut(lx, lx, lx)
1495 integer :: e, i, j, k
1496
1497 !$omp do
1498 do e = 1, n
1499 do j = 1, lx * lx
1500 do i = 1, lx
1501 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1502 + dx(i,2) * u(2,j,1,e) &
1503 + dx(i,3) * u(3,j,1,e) &
1504 + dx(i,4) * u(4,j,1,e) &
1505 + dx(i,5) * u(5,j,1,e) &
1506 + dx(i,6) * u(6,j,1,e)
1507 end do
1508 end do
1509
1510 do k = 1, lx
1511 do j = 1, lx
1512 do i = 1, lx
1513 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1514 + dy(j,2) * u(i,2,k,e) &
1515 + dy(j,3) * u(i,3,k,e) &
1516 + dy(j,4) * u(i,4,k,e) &
1517 + dy(j,5) * u(i,5,k,e) &
1518 + dy(j,6) * u(i,6,k,e)
1519 end do
1520 end do
1521 end do
1522
1523 do k = 1, lx
1524 do i = 1, lx*lx
1525 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1526 + dz(k,2) * u(i,1,2,e) &
1527 + dz(k,3) * u(i,1,3,e) &
1528 + dz(k,4) * u(i,1,4,e) &
1529 + dz(k,5) * u(i,1,5,e) &
1530 + dz(k,6) * u(i,1,6,e)
1531 end do
1532 end do
1533
1534 do i = 1, lx*lx*lx
1535 ur(i,1,1) = h1(i,1,1,e) &
1536 * ( g11(i,1,1,e) * wur(i,1,1) &
1537 + g12(i,1,1,e) * wus(i,1,1) &
1538 + g13(i,1,1,e) * wut(i,1,1) )
1539 us(i,1,1) = h1(i,1,1,e) &
1540 * ( g12(i,1,1,e) * wur(i,1,1) &
1541 + g22(i,1,1,e) * wus(i,1,1) &
1542 + g23(i,1,1,e) * wut(i,1,1) )
1543 ut(i,1,1) = h1(i,1,1,e) &
1544 * ( g13(i,1,1,e) * wur(i,1,1) &
1545 + g23(i,1,1,e) * wus(i,1,1) &
1546 + g33(i,1,1,e) * wut(i,1,1) )
1547 end do
1548
1549 do j = 1, lx*lx
1550 do i = 1, lx
1551 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1552 + dxt(i,2) * ur(2,j,1) &
1553 + dxt(i,3) * ur(3,j,1) &
1554 + dxt(i,4) * ur(4,j,1) &
1555 + dxt(i,5) * ur(5,j,1) &
1556 + dxt(i,6) * ur(6,j,1)
1557 end do
1558 end do
1559
1560 do k = 1, lx
1561 do j = 1, lx
1562 do i = 1, lx
1563 w(i,j,k,e) = w(i,j,k,e) &
1564 + dyt(j,1) * us(i,1,k) &
1565 + dyt(j,2) * us(i,2,k) &
1566 + dyt(j,3) * us(i,3,k) &
1567 + dyt(j,4) * us(i,4,k) &
1568 + dyt(j,5) * us(i,5,k) &
1569 + dyt(j,6) * us(i,6,k)
1570 end do
1571 end do
1572 end do
1573
1574 do k = 1, lx
1575 do i = 1, lx*lx
1576 w(i,1,k,e) = w(i,1,k,e) &
1577 + dzt(k,1) * ut(i,1,1) &
1578 + dzt(k,2) * ut(i,1,2) &
1579 + dzt(k,3) * ut(i,1,3) &
1580 + dzt(k,4) * ut(i,1,4) &
1581 + dzt(k,5) * ut(i,1,5) &
1582 + dzt(k,6) * ut(i,1,6)
1583 end do
1584 end do
1585
1586 end do
1587 !$omp end do
1588 end subroutine ax_helm_lx6
1589
1590 subroutine ax_helm_lx5(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1591 h1, G11, G22, G33, G12, G13, G23, n)
1592 integer, parameter :: lx = 5
1593 integer, intent(in) :: n
1594 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1595 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1596 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1597 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1598 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1599 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1600 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1601 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1602 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1603 real(kind=rp), intent(in) :: dx(lx, lx)
1604 real(kind=rp), intent(in) :: dy(lx, lx)
1605 real(kind=rp), intent(in) :: dz(lx, lx)
1606 real(kind=rp), intent(in) :: dxt(lx, lx)
1607 real(kind=rp), intent(in) :: dyt(lx, lx)
1608 real(kind=rp), intent(in) :: dzt(lx, lx)
1609 real(kind=rp) :: ur(lx, lx, lx)
1610 real(kind=rp) :: us(lx, lx, lx)
1611 real(kind=rp) :: ut(lx, lx, lx)
1612 real(kind=rp) :: wur(lx, lx, lx)
1613 real(kind=rp) :: wus(lx, lx, lx)
1614 real(kind=rp) :: wut(lx, lx, lx)
1615 integer :: e, i, j, k
1616
1617 !$omp do
1618 do e = 1, n
1619 do j = 1, lx * lx
1620 do i = 1, lx
1621 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1622 + dx(i,2) * u(2,j,1,e) &
1623 + dx(i,3) * u(3,j,1,e) &
1624 + dx(i,4) * u(4,j,1,e) &
1625 + dx(i,5) * u(5,j,1,e)
1626 end do
1627 end do
1628
1629 do k = 1, lx
1630 do j = 1, lx
1631 do i = 1, lx
1632 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1633 + dy(j,2) * u(i,2,k,e) &
1634 + dy(j,3) * u(i,3,k,e) &
1635 + dy(j,4) * u(i,4,k,e) &
1636 + dy(j,5) * u(i,5,k,e)
1637 end do
1638 end do
1639 end do
1640
1641 do k = 1, lx
1642 do i = 1, lx*lx
1643 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1644 + dz(k,2) * u(i,1,2,e) &
1645 + dz(k,3) * u(i,1,3,e) &
1646 + dz(k,4) * u(i,1,4,e) &
1647 + dz(k,5) * u(i,1,5,e)
1648 end do
1649 end do
1650
1651 do i = 1, lx*lx*lx
1652 ur(i,1,1) = h1(i,1,1,e) &
1653 * ( g11(i,1,1,e) * wur(i,1,1) &
1654 + g12(i,1,1,e) * wus(i,1,1) &
1655 + g13(i,1,1,e) * wut(i,1,1) )
1656 us(i,1,1) = h1(i,1,1,e) &
1657 * ( g12(i,1,1,e) * wur(i,1,1) &
1658 + g22(i,1,1,e) * wus(i,1,1) &
1659 + g23(i,1,1,e) * wut(i,1,1) )
1660 ut(i,1,1) = h1(i,1,1,e) &
1661 * ( g13(i,1,1,e) * wur(i,1,1) &
1662 + g23(i,1,1,e) * wus(i,1,1) &
1663 + g33(i,1,1,e) * wut(i,1,1) )
1664 end do
1665
1666 do j = 1, lx*lx
1667 do i = 1, lx
1668 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1669 + dxt(i,2) * ur(2,j,1) &
1670 + dxt(i,3) * ur(3,j,1) &
1671 + dxt(i,4) * ur(4,j,1) &
1672 + dxt(i,5) * ur(5,j,1)
1673 end do
1674 end do
1675
1676 do k = 1, lx
1677 do j = 1, lx
1678 do i = 1, lx
1679 w(i,j,k,e) = w(i,j,k,e) &
1680 + dyt(j,1) * us(i,1,k) &
1681 + dyt(j,2) * us(i,2,k) &
1682 + dyt(j,3) * us(i,3,k) &
1683 + dyt(j,4) * us(i,4,k) &
1684 + dyt(j,5) * us(i,5,k)
1685 end do
1686 end do
1687 end do
1688
1689 do k = 1, lx
1690 do i = 1, lx*lx
1691 w(i,1,k,e) = w(i,1,k,e) &
1692 + dzt(k,1) * ut(i,1,1) &
1693 + dzt(k,2) * ut(i,1,2) &
1694 + dzt(k,3) * ut(i,1,3) &
1695 + dzt(k,4) * ut(i,1,4) &
1696 + dzt(k,5) * ut(i,1,5)
1697 end do
1698 end do
1699
1700 end do
1701 !$omp end do
1702 end subroutine ax_helm_lx5
1703
1704 subroutine ax_helm_lx4(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1705 h1, G11, G22, G33, G12, G13, G23, n)
1706 integer, parameter :: lx = 4
1707 integer, intent(in) :: n
1708 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1709 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1710 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1711 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1712 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1713 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1714 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1715 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1716 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1717 real(kind=rp), intent(in) :: dx(lx, lx)
1718 real(kind=rp), intent(in) :: dy(lx, lx)
1719 real(kind=rp), intent(in) :: dz(lx, lx)
1720 real(kind=rp), intent(in) :: dxt(lx, lx)
1721 real(kind=rp), intent(in) :: dyt(lx, lx)
1722 real(kind=rp), intent(in) :: dzt(lx, lx)
1723 real(kind=rp) :: ur(lx, lx, lx)
1724 real(kind=rp) :: us(lx, lx, lx)
1725 real(kind=rp) :: ut(lx, lx, lx)
1726 real(kind=rp) :: wur(lx, lx, lx)
1727 real(kind=rp) :: wus(lx, lx, lx)
1728 real(kind=rp) :: wut(lx, lx, lx)
1729 integer :: e, i, j, k
1730
1731 !$omp do
1732 do e = 1, n
1733 do j = 1, lx * lx
1734 do i = 1, lx
1735 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1736 + dx(i,2) * u(2,j,1,e) &
1737 + dx(i,3) * u(3,j,1,e) &
1738 + dx(i,4) * u(4,j,1,e)
1739 end do
1740 end do
1741
1742 do k = 1, lx
1743 do j = 1, lx
1744 do i = 1, lx
1745 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1746 + dy(j,2) * u(i,2,k,e) &
1747 + dy(j,3) * u(i,3,k,e) &
1748 + dy(j,4) * u(i,4,k,e)
1749 end do
1750 end do
1751 end do
1752
1753 do k = 1, lx
1754 do i = 1, lx*lx
1755 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1756 + dz(k,2) * u(i,1,2,e) &
1757 + dz(k,3) * u(i,1,3,e) &
1758 + dz(k,4) * u(i,1,4,e)
1759 end do
1760 end do
1761
1762 do i = 1, lx*lx*lx
1763 ur(i,1,1) = h1(i,1,1,e) &
1764 * ( g11(i,1,1,e) * wur(i,1,1) &
1765 + g12(i,1,1,e) * wus(i,1,1) &
1766 + g13(i,1,1,e) * wut(i,1,1) )
1767 us(i,1,1) = h1(i,1,1,e) &
1768 * ( g12(i,1,1,e) * wur(i,1,1) &
1769 + g22(i,1,1,e) * wus(i,1,1) &
1770 + g23(i,1,1,e) * wut(i,1,1) )
1771 ut(i,1,1) = h1(i,1,1,e) &
1772 * ( g13(i,1,1,e) * wur(i,1,1) &
1773 + g23(i,1,1,e) * wus(i,1,1) &
1774 + g33(i,1,1,e) * wut(i,1,1) )
1775 end do
1776
1777 do j = 1, lx*lx
1778 do i = 1, lx
1779 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1780 + dxt(i,2) * ur(2,j,1) &
1781 + dxt(i,3) * ur(3,j,1) &
1782 + dxt(i,4) * ur(4,j,1)
1783 end do
1784 end do
1785
1786 do k = 1, lx
1787 do j = 1, lx
1788 do i = 1, lx
1789 w(i,j,k,e) = w(i,j,k,e) &
1790 + dyt(j,1) * us(i,1,k) &
1791 + dyt(j,2) * us(i,2,k) &
1792 + dyt(j,3) * us(i,3,k) &
1793 + dyt(j,4) * us(i,4,k)
1794 end do
1795 end do
1796 end do
1797
1798 do k = 1, lx
1799 do i = 1, lx*lx
1800 w(i,1,k,e) = w(i,1,k,e) &
1801 + dzt(k,1) * ut(i,1,1) &
1802 + dzt(k,2) * ut(i,1,2) &
1803 + dzt(k,3) * ut(i,1,3) &
1804 + dzt(k,4) * ut(i,1,4)
1805 end do
1806 end do
1807
1808 end do
1809 !$omp end do
1810 end subroutine ax_helm_lx4
1811
1812 subroutine ax_helm_lx3(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1813 h1, G11, G22, G33, G12, G13, G23, n)
1814 integer, parameter :: lx = 3
1815 integer, intent(in) :: n
1816 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1817 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1818 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1819 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1820 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1821 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1822 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1823 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1824 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1825 real(kind=rp), intent(in) :: dx(lx, lx)
1826 real(kind=rp), intent(in) :: dy(lx, lx)
1827 real(kind=rp), intent(in) :: dz(lx, lx)
1828 real(kind=rp), intent(in) :: dxt(lx, lx)
1829 real(kind=rp), intent(in) :: dyt(lx, lx)
1830 real(kind=rp), intent(in) :: dzt(lx, lx)
1831 real(kind=rp) :: ur(lx, lx, lx)
1832 real(kind=rp) :: us(lx, lx, lx)
1833 real(kind=rp) :: ut(lx, lx, lx)
1834 real(kind=rp) :: wur(lx, lx, lx)
1835 real(kind=rp) :: wus(lx, lx, lx)
1836 real(kind=rp) :: wut(lx, lx, lx)
1837 integer :: e, i, j, k
1838
1839 !$omp do
1840 do e = 1, n
1841 do j = 1, lx * lx
1842 do i = 1, lx
1843 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1844 + dx(i,2) * u(2,j,1,e) &
1845 + dx(i,3) * u(3,j,1,e)
1846 end do
1847 end do
1848
1849 do k = 1, lx
1850 do j = 1, lx
1851 do i = 1, lx
1852 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1853 + dy(j,2) * u(i,2,k,e) &
1854 + dy(j,3) * u(i,3,k,e)
1855 end do
1856 end do
1857 end do
1858
1859 do k = 1, lx
1860 do i = 1, lx*lx
1861 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1862 + dz(k,2) * u(i,1,2,e) &
1863 + dz(k,3) * u(i,1,3,e)
1864 end do
1865 end do
1866
1867 do i = 1, lx*lx*lx
1868 ur(i,1,1) = h1(i,1,1,e) &
1869 * ( g11(i,1,1,e) * wur(i,1,1) &
1870 + g12(i,1,1,e) * wus(i,1,1) &
1871 + g13(i,1,1,e) * wut(i,1,1) )
1872 us(i,1,1) = h1(i,1,1,e) &
1873 * ( g12(i,1,1,e) * wur(i,1,1) &
1874 + g22(i,1,1,e) * wus(i,1,1) &
1875 + g23(i,1,1,e) * wut(i,1,1) )
1876 ut(i,1,1) = h1(i,1,1,e) &
1877 * ( g13(i,1,1,e) * wur(i,1,1) &
1878 + g23(i,1,1,e) * wus(i,1,1) &
1879 + g33(i,1,1,e) * wut(i,1,1) )
1880 end do
1881
1882 do j = 1, lx*lx
1883 do i = 1, lx
1884 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1885 + dxt(i,2) * ur(2,j,1) &
1886 + dxt(i,3) * ur(3,j,1)
1887 end do
1888 end do
1889
1890 do k = 1, lx
1891 do j = 1, lx
1892 do i = 1, lx
1893 w(i,j,k,e) = w(i,j,k,e) &
1894 + dyt(j,1) * us(i,1,k) &
1895 + dyt(j,2) * us(i,2,k) &
1896 + dyt(j,3) * us(i,3,k)
1897 end do
1898 end do
1899 end do
1900
1901 do k = 1, lx
1902 do i = 1, lx*lx
1903 w(i,1,k,e) = w(i,1,k,e) &
1904 + dzt(k,1) * ut(i,1,1) &
1905 + dzt(k,2) * ut(i,1,2) &
1906 + dzt(k,3) * ut(i,1,3)
1907 end do
1908 end do
1909
1910 end do
1911 !$omp end do
1912 end subroutine ax_helm_lx3
1913
1914 subroutine ax_helm_lx2(w, u, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1915 h1, G11, G22, G33, G12, G13, G23, n)
1916 integer, parameter :: lx = 2
1917 integer, intent(in) :: n
1918 real(kind=rp), intent(inout) :: w(lx, lx, lx, n)
1919 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1920 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1921 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1922 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1923 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1924 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1925 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1926 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1927 real(kind=rp), intent(in) :: dx(lx, lx)
1928 real(kind=rp), intent(in) :: dy(lx, lx)
1929 real(kind=rp), intent(in) :: dz(lx, lx)
1930 real(kind=rp), intent(in) :: dxt(lx, lx)
1931 real(kind=rp), intent(in) :: dyt(lx, lx)
1932 real(kind=rp), intent(in) :: dzt(lx, lx)
1933 real(kind=rp) :: ur(lx, lx, lx)
1934 real(kind=rp) :: us(lx, lx, lx)
1935 real(kind=rp) :: ut(lx, lx, lx)
1936 real(kind=rp) :: wur(lx, lx, lx)
1937 real(kind=rp) :: wus(lx, lx, lx)
1938 real(kind=rp) :: wut(lx, lx, lx)
1939 integer :: e, i, j, k
1940
1941 !$omp do
1942 do e = 1, n
1943 do j = 1, lx * lx
1944 do i = 1, lx
1945 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1946 + dx(i,2) * u(2,j,1,e)
1947 end do
1948 end do
1949
1950 do k = 1, lx
1951 do j = 1, lx
1952 do i = 1, lx
1953 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1954 + dy(j,2) * u(i,2,k,e)
1955 end do
1956 end do
1957 end do
1958
1959 do k = 1, lx
1960 do i = 1, lx*lx
1961 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1962 + dz(k,2) * u(i,1,2,e)
1963 end do
1964 end do
1965
1966 do i = 1, lx*lx*lx
1967 ur(i,1,1) = h1(i,1,1,e) &
1968 * ( g11(i,1,1,e) * wur(i,1,1) &
1969 + g12(i,1,1,e) * wus(i,1,1) &
1970 + g13(i,1,1,e) * wut(i,1,1) )
1971 us(i,1,1) = h1(i,1,1,e) &
1972 * ( g12(i,1,1,e) * wur(i,1,1) &
1973 + g22(i,1,1,e) * wus(i,1,1) &
1974 + g23(i,1,1,e) * wut(i,1,1) )
1975 ut(i,1,1) = h1(i,1,1,e) &
1976 * ( g13(i,1,1,e) * wur(i,1,1) &
1977 + g23(i,1,1,e) * wus(i,1,1) &
1978 + g33(i,1,1,e) * wut(i,1,1) )
1979 end do
1980
1981 do j = 1, lx*lx
1982 do i = 1, lx
1983 w(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1984 + dxt(i,2) * ur(2,j,1)
1985 end do
1986 end do
1987
1988 do k = 1, lx
1989 do j = 1, lx
1990 do i = 1, lx
1991 w(i,j,k,e) = w(i,j,k,e) &
1992 + dyt(j,1) * us(i,1,k) &
1993 + dyt(j,2) * us(i,2,k)
1994 end do
1995 end do
1996 end do
1997
1998 do k = 1, lx
1999 do i = 1, lx*lx
2000 w(i,1,k,e) = w(i,1,k,e) &
2001 + dzt(k,1) * ut(i,1,1) &
2002 + dzt(k,2) * ut(i,1,2)
2003 end do
2004 end do
2005
2006 end do
2007 !$omp end do
2008 end subroutine ax_helm_lx2
2009
2010end module ax_helm_cpu
ax_helm_cpu
Definition ax_helm_cpu.f90:33
ax_helm_cpu::ax_helm_compute
subroutine ax_helm_compute(this, w, u, coef, msh, xh)
Compute product, taking three components of a vector field in an uncoupled manner.
Definition ax_helm_cpu.f90:89
ax_helm
Definition ax_helm.f90:33
ax_helm::ax_helm_compute_vector
subroutine ax_helm_compute_vector(this, au, av, aw, u, v, w, coef, msh, xh)
Definition ax_helm.f90:63
coefs
Coefficients.
Definition coef.f90:34
mesh
Defines a mesh.
Definition mesh.f90:34
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:14
space
Defines a function space.
Definition space.f90:34
ax_helm::ax_helm_t
Matrix-vector product for a Helmholtz problem.
Definition ax_helm.f90:44
ax_helm_cpu::ax_helm_cpu_t
CPU matrix-vector product for a Helmholtz problem.
Definition ax_helm_cpu.f90:43
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:135
mesh::mesh_t
Definition mesh.f90:73
space::space_t
The function space for the SEM solution fields.
Definition space.f90:64