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