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