Neko 0.9.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
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!
33module 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
50contains
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
1945end module ax_helm_cpu
ax_helm_cpu
Definition ax_helm_cpu.f90:33
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_compute
subroutine ax_helm_compute(w, u, coef, msh, xh)
Compute the product.
Definition ax_helm_cpu.f90:61
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_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_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_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_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_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_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_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_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_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_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_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_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
Definition ax_helm.f90:33
coefs
Coefficients.
Definition coef.f90:34
math
Definition math.f90:60
math::addcol4
subroutine, public addcol4(a, b, c, d, n)
Returns .
Definition math.f90:815
mesh
Defines a mesh.
Definition mesh.f90:34
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
space
Defines a function space.
Definition space.f90:34
ax_helm::ax_helm_t
Matrix-vector product for a Helmholtz problem.
Definition ax_helm.f90:44
ax_helm_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