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