Neko  0.9.99
A portable framework for high-order spectral element flow simulations
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 !
33 module 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 
48 contains
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 
1850 end module ax_helm_sx
ax_helm_sx
Definition: ax_helm_sx.f90:33
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::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_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_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_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::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_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_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_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_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_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_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_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_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_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
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:815
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