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