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