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