Neko 1.99.6
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
ax_helm_vector_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 2026, 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!
33submodule(ax_helm_cpu) ax_helm_vector_cpu
34 implicit none
35
36contains
37
51 module subroutine ax_helm_compute_vector(this, au, av, aw, &
52 u, v, w, coef, msh, xh)
53 class(ax_helm_cpu_t), intent(in) :: this
54 type(mesh_t), intent(in) :: msh
55 type(space_t), intent(in) :: Xh
56 type(coef_t), intent(in) :: coef
57 real(kind=rp), intent(inout) :: au(xh%lx, xh%ly, xh%lz, msh%nelv)
58 real(kind=rp), intent(inout) :: av(xh%lx, xh%ly, xh%lz, msh%nelv)
59 real(kind=rp), intent(inout) :: aw(xh%lx, xh%ly, xh%lz, msh%nelv)
60 real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, msh%nelv)
61 real(kind=rp), intent(in) :: v(xh%lx, xh%ly, xh%lz, msh%nelv)
62 real(kind=rp), intent(in) :: w(xh%lx, xh%ly, xh%lz, msh%nelv)
63 integer :: i
64
65 !$omp parallel
66 select case(xh%lx)
67 case (14)
68 call ax_helm_vector_lx14(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
69 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
70 coef%G12, coef%G13, coef%G23, msh%nelv)
71 case (13)
72 call ax_helm_vector_lx13(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
73 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
74 coef%G12, coef%G13, coef%G23, msh%nelv)
75 case (12)
76 call ax_helm_vector_lx12(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
77 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
78 coef%G12, coef%G13, coef%G23, msh%nelv)
79 case (11)
80 call ax_helm_vector_lx11(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
81 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
82 coef%G12, coef%G13, coef%G23, msh%nelv)
83 case (10)
84 call ax_helm_vector_lx10(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
85 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
86 coef%G12, coef%G13, coef%G23, msh%nelv)
87 case (9)
88 call ax_helm_vector_lx9(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
89 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
90 coef%G12, coef%G13, coef%G23, msh%nelv)
91 case (8)
92 call ax_helm_vector_lx8(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
93 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
94 coef%G12, coef%G13, coef%G23, msh%nelv)
95 case (7)
96 call ax_helm_vector_lx7(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
97 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
98 coef%G12, coef%G13, coef%G23, msh%nelv)
99 case (6)
100 call ax_helm_vector_lx6(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
101 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
102 coef%G12, coef%G13, coef%G23, msh%nelv)
103 case (5)
104 call ax_helm_vector_lx5(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
105 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
106 coef%G12, coef%G13, coef%G23, msh%nelv)
107 case (4)
108 call ax_helm_vector_lx4(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
109 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
110 coef%G12, coef%G13, coef%G23, msh%nelv)
111 case (3)
112 call ax_helm_vector_lx3(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
113 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
114 coef%G12, coef%G13, coef%G23, msh%nelv)
115 case (2)
116 call ax_helm_vector_lx2(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
117 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
118 coef%G12, coef%G13, coef%G23, msh%nelv)
119 case default
120 call ax_helm_vector_lx(au, av, aw, u, v, w, xh%dx, xh%dy, xh%dz, &
121 xh%dxt, xh%dyt, xh%dzt, coef%h1, coef%G11, coef%G22, coef%G33, &
122 coef%G12, coef%G13, coef%G23, msh%nelv, xh%lx)
123 end select
124
125 if (coef%ifh2) then
126 !$omp do private(i)
127 do i = 1, coef%dof%size()
128 au(i,1,1,1) = au(i,1,1,1) + &
129 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * u(i,1,1,1)
130 av(i,1,1,1) = av(i,1,1,1) + &
131 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * v(i,1,1,1)
132 aw(i,1,1,1) = aw(i,1,1,1) + &
133 coef%h2(i,1,1,1) * coef%B(i,1,1,1) * w(i,1,1,1)
134 end do
135 !$omp end do
136 end if
137 !$omp end parallel
138
139 end subroutine ax_helm_compute_vector
140
153 subroutine ax_helm_vector_lx(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
154 h1, G11, G22, G33, G12, G13, G23, n, lx)
155 integer, intent(in) :: n, lx
156 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
157 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
158 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
159 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
160 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
161 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
162 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
163 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
164 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
165 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
166 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
167 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
168 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
169 real(kind=rp), intent(in) :: dx(lx, lx)
170 real(kind=rp), intent(in) :: dy(lx, lx)
171 real(kind=rp), intent(in) :: dz(lx, lx)
172 real(kind=rp), intent(in) :: dxt(lx, lx)
173 real(kind=rp), intent(in) :: dyt(lx, lx)
174 real(kind=rp), intent(in) :: dzt(lx, lx)
175 real(kind=rp) :: ur(lx, lx, lx)
176 real(kind=rp) :: us(lx, lx, lx)
177 real(kind=rp) :: ut(lx, lx, lx)
178 real(kind=rp) :: vr(lx, lx, lx)
179 real(kind=rp) :: vs(lx, lx, lx)
180 real(kind=rp) :: vt(lx, lx, lx)
181 real(kind=rp) :: wr(lx, lx, lx)
182 real(kind=rp) :: ws(lx, lx, lx)
183 real(kind=rp) :: wt(lx, lx, lx)
184 real(kind=rp) :: wur(lx, lx, lx)
185 real(kind=rp) :: wvr(lx, lx, lx)
186 real(kind=rp) :: wwr(lx, lx, lx)
187 real(kind=rp) :: wus(lx, lx, lx)
188 real(kind=rp) :: wvs(lx, lx, lx)
189 real(kind=rp) :: wws(lx, lx, lx)
190 real(kind=rp) :: wut(lx, lx, lx)
191 real(kind=rp) :: wvt(lx, lx, lx)
192 real(kind=rp) :: wwt(lx, lx, lx)
193 real(kind=rp) :: t1, t2, t3
194 integer :: e, i, j, k, l
195
196 !$omp do
197 do e = 1, n
198 do j = 1, lx * lx
199 do i = 1, lx
200 t1 = 0.0_rp
201 t2 = 0.0_rp
202 t3 = 0.0_rp
203 do k = 1, lx
204 t1 = t1 + dx(i,k) * u(k,j,1,e)
205 t2 = t2 + dx(i,k) * v(k,j,1,e)
206 t3 = t3 + dx(i,k) * w(k,j,1,e)
207 end do
208 wur(i,j,1) = t1
209 wvr(i,j,1) = t2
210 wwr(i,j,1) = t3
211 end do
212 end do
213
214 do k = 1, lx
215 do j = 1, lx
216 do i = 1, lx
217 t1 = 0.0_rp
218 t2 = 0.0_rp
219 t3 = 0.0_rp
220 do l = 1, lx
221 t1 = t1 + dy(j,l) * u(i,l,k,e)
222 t2 = t2 + dy(j,l) * v(i,l,k,e)
223 t3 = t3 + dy(j,l) * w(i,l,k,e)
224 end do
225 wus(i,j,k) = t1
226 wvs(i,j,k) = t2
227 wws(i,j,k) = t3
228 end do
229 end do
230 end do
231
232 do k = 1, lx
233 do i = 1, lx*lx
234 t1 = 0.0_rp
235 t2 = 0.0_rp
236 t3 = 0.0_rp
237 do l = 1, lx
238 t1 = t1 + dz(k,l) * u(i,1,l,e)
239 t2 = t2 + dz(k,l) * v(i,1,l,e)
240 t3 = t3 + dz(k,l) * w(i,1,l,e)
241 end do
242 wut(i,1,k) = t1
243 wvt(i,1,k) = t2
244 wwt(i,1,k) = t3
245 end do
246 end do
247
248 do i = 1, lx*lx*lx
249 ur(i,1,1) = h1(i,1,1,e) &
250 * ( g11(i,1,1,e) * wur(i,1,1) &
251 + g12(i,1,1,e) * wus(i,1,1) &
252 + g13(i,1,1,e) * wut(i,1,1) )
253 us(i,1,1) = h1(i,1,1,e) &
254 * ( g12(i,1,1,e) * wur(i,1,1) &
255 + g22(i,1,1,e) * wus(i,1,1) &
256 + g23(i,1,1,e) * wut(i,1,1) )
257 ut(i,1,1) = h1(i,1,1,e) &
258 * ( g13(i,1,1,e) * wur(i,1,1) &
259 + g23(i,1,1,e) * wus(i,1,1) &
260 + g33(i,1,1,e) * wut(i,1,1) )
261
262 vr(i,1,1) = h1(i,1,1,e) &
263 * ( g11(i,1,1,e) * wvr(i,1,1) &
264 + g12(i,1,1,e) * wvs(i,1,1) &
265 + g13(i,1,1,e) * wvt(i,1,1) )
266 vs(i,1,1) = h1(i,1,1,e) &
267 * ( g12(i,1,1,e) * wvr(i,1,1) &
268 + g22(i,1,1,e) * wvs(i,1,1) &
269 + g23(i,1,1,e) * wvt(i,1,1) )
270 vt(i,1,1) = h1(i,1,1,e) &
271 * ( g13(i,1,1,e) * wvr(i,1,1) &
272 + g23(i,1,1,e) * wvs(i,1,1) &
273 + g33(i,1,1,e) * wvt(i,1,1) )
274
275 wr(i,1,1) = h1(i,1,1,e) &
276 * ( g11(i,1,1,e) * wwr(i,1,1) &
277 + g12(i,1,1,e) * wws(i,1,1) &
278 + g13(i,1,1,e) * wwt(i,1,1) )
279 ws(i,1,1) = h1(i,1,1,e) &
280 * ( g12(i,1,1,e) * wwr(i,1,1) &
281 + g22(i,1,1,e) * wws(i,1,1) &
282 + g23(i,1,1,e) * wwt(i,1,1) )
283 wt(i,1,1) = h1(i,1,1,e) &
284 * ( g13(i,1,1,e) * wwr(i,1,1) &
285 + g23(i,1,1,e) * wws(i,1,1) &
286 + g33(i,1,1,e) * wwt(i,1,1) )
287 end do
288
289 do j = 1, lx*lx
290 do i = 1, lx
291 t1 = 0.0_rp
292 t2 = 0.0_rp
293 t3 = 0.0_rp
294 do k = 1, lx
295 t1 = t1 + dxt(i,k) * ur(k,j,1)
296 t2 = t2 + dxt(i,k) * vr(k,j,1)
297 t3 = t3 + dxt(i,k) * wr(k,j,1)
298 end do
299 au(i,j,1,e) = t1
300 av(i,j,1,e) = t2
301 aw(i,j,1,e) = t3
302 end do
303 end do
304
305 do k = 1, lx
306 do j = 1, lx
307 do i = 1, lx
308 t1 = 0.0_rp
309 t2 = 0.0_rp
310 t3 = 0.0_rp
311 do l = 1, lx
312 t1 = t1 + dyt(j,l) * us(i,l,k)
313 t2 = t2 + dyt(j,l) * vs(i,l,k)
314 t3 = t3 + dyt(j,l) * ws(i,l,k)
315 end do
316 au(i,j,k,e) = au(i,j,k,e) + t1
317 av(i,j,k,e) = av(i,j,k,e) + t2
318 aw(i,j,k,e) = aw(i,j,k,e) + t3
319 end do
320 end do
321 end do
322
323 do k = 1, lx
324 do i = 1, lx*lx
325 t1 = 0.0_rp
326 t2 = 0.0_rp
327 t3 = 0.0_rp
328 do l = 1, lx
329 t1 = t1 + dzt(k,l) * ut(i,1,l)
330 t2 = t2 + dzt(k,l) * vt(i,1,l)
331 t3 = t3 + dzt(k,l) * wt(i,1,l)
332 end do
333 au(i,1,k,e) = au(i,1,k,e) + t1
334 av(i,1,k,e) = av(i,1,k,e) + t2
335 aw(i,1,k,e) = aw(i,1,k,e) + t3
336 end do
337 end do
338
339 end do
340 !$omp end do
341 end subroutine ax_helm_vector_lx
342
343 subroutine ax_helm_vector_lx14(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
344 h1, G11, G22, G33, G12, G13, G23, n)
345 integer, parameter :: lx = 14
346 integer, intent(in) :: n
347 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
348 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
349 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
350 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
351 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
352 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
353 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
354 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
355 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
356 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
357 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
358 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
359 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
360 real(kind=rp), intent(in) :: dx(lx, lx)
361 real(kind=rp), intent(in) :: dy(lx, lx)
362 real(kind=rp), intent(in) :: dz(lx, lx)
363 real(kind=rp), intent(in) :: dxt(lx, lx)
364 real(kind=rp), intent(in) :: dyt(lx, lx)
365 real(kind=rp), intent(in) :: dzt(lx, lx)
366 real(kind=rp) :: ur(lx, lx, lx)
367 real(kind=rp) :: us(lx, lx, lx)
368 real(kind=rp) :: ut(lx, lx, lx)
369 real(kind=rp) :: vr(lx, lx, lx)
370 real(kind=rp) :: vs(lx, lx, lx)
371 real(kind=rp) :: vt(lx, lx, lx)
372 real(kind=rp) :: wr(lx, lx, lx)
373 real(kind=rp) :: ws(lx, lx, lx)
374 real(kind=rp) :: wt(lx, lx, lx)
375 real(kind=rp) :: wur(lx, lx, lx)
376 real(kind=rp) :: wus(lx, lx, lx)
377 real(kind=rp) :: wut(lx, lx, lx)
378 real(kind=rp) :: wvr(lx, lx, lx)
379 real(kind=rp) :: wvs(lx, lx, lx)
380 real(kind=rp) :: wvt(lx, lx, lx)
381 real(kind=rp) :: wwr(lx, lx, lx)
382 real(kind=rp) :: wws(lx, lx, lx)
383 real(kind=rp) :: wwt(lx, lx, lx)
384 integer :: e, i, j, k
385
386 !$omp do
387 do e = 1, n
388 do j = 1, lx * lx
389 do i = 1, lx
390 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
391 + dx(i,2) * u(2,j,1,e) &
392 + dx(i,3) * u(3,j,1,e) &
393 + dx(i,4) * u(4,j,1,e) &
394 + dx(i,5) * u(5,j,1,e) &
395 + dx(i,6) * u(6,j,1,e) &
396 + dx(i,7) * u(7,j,1,e) &
397 + dx(i,8) * u(8,j,1,e) &
398 + dx(i,9) * u(9,j,1,e) &
399 + dx(i,10) * u(10,j,1,e) &
400 + dx(i,11) * u(11,j,1,e) &
401 + dx(i,12) * u(12,j,1,e) &
402 + dx(i,13) * u(13,j,1,e) &
403 + dx(i,14) * u(14,j,1,e)
404
405 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
406 + dx(i,2) * v(2,j,1,e) &
407 + dx(i,3) * v(3,j,1,e) &
408 + dx(i,4) * v(4,j,1,e) &
409 + dx(i,5) * v(5,j,1,e) &
410 + dx(i,6) * v(6,j,1,e) &
411 + dx(i,7) * v(7,j,1,e) &
412 + dx(i,8) * v(8,j,1,e) &
413 + dx(i,9) * v(9,j,1,e) &
414 + dx(i,10) * v(10,j,1,e) &
415 + dx(i,11) * v(11,j,1,e) &
416 + dx(i,12) * v(12,j,1,e) &
417 + dx(i,13) * v(13,j,1,e) &
418 + dx(i,14) * v(14,j,1,e)
419
420 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
421 + dx(i,2) * w(2,j,1,e) &
422 + dx(i,3) * w(3,j,1,e) &
423 + dx(i,4) * w(4,j,1,e) &
424 + dx(i,5) * w(5,j,1,e) &
425 + dx(i,6) * w(6,j,1,e) &
426 + dx(i,7) * w(7,j,1,e) &
427 + dx(i,8) * w(8,j,1,e) &
428 + dx(i,9) * w(9,j,1,e) &
429 + dx(i,10) * w(10,j,1,e) &
430 + dx(i,11) * w(11,j,1,e) &
431 + dx(i,12) * w(12,j,1,e) &
432 + dx(i,13) * w(13,j,1,e) &
433 + dx(i,14) * w(14,j,1,e)
434 end do
435 end do
436
437 do k = 1, lx
438 do j = 1, lx
439 do i = 1, lx
440 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
441 + dy(j,2) * u(i,2,k,e) &
442 + dy(j,3) * u(i,3,k,e) &
443 + dy(j,4) * u(i,4,k,e) &
444 + dy(j,5) * u(i,5,k,e) &
445 + dy(j,6) * u(i,6,k,e) &
446 + dy(j,7) * u(i,7,k,e) &
447 + dy(j,8) * u(i,8,k,e) &
448 + dy(j,9) * u(i,9,k,e) &
449 + dy(j,10) * u(i,10,k,e) &
450 + dy(j,11) * u(i,11,k,e) &
451 + dy(j,12) * u(i,12,k,e) &
452 + dy(j,13) * u(i,13,k,e) &
453 + dy(j,14) * u(i,14,k,e)
454
455 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
456 + dy(j,2) * v(i,2,k,e) &
457 + dy(j,3) * v(i,3,k,e) &
458 + dy(j,4) * v(i,4,k,e) &
459 + dy(j,5) * v(i,5,k,e) &
460 + dy(j,6) * v(i,6,k,e) &
461 + dy(j,7) * v(i,7,k,e) &
462 + dy(j,8) * v(i,8,k,e) &
463 + dy(j,9) * v(i,9,k,e) &
464 + dy(j,10) * v(i,10,k,e) &
465 + dy(j,11) * v(i,11,k,e) &
466 + dy(j,12) * v(i,12,k,e) &
467 + dy(j,13) * v(i,13,k,e) &
468 + dy(j,14) * v(i,14,k,e)
469
470 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
471 + dy(j,2) * w(i,2,k,e) &
472 + dy(j,3) * w(i,3,k,e) &
473 + dy(j,4) * w(i,4,k,e) &
474 + dy(j,5) * w(i,5,k,e) &
475 + dy(j,6) * w(i,6,k,e) &
476 + dy(j,7) * w(i,7,k,e) &
477 + dy(j,8) * w(i,8,k,e) &
478 + dy(j,9) * w(i,9,k,e) &
479 + dy(j,10) * w(i,10,k,e) &
480 + dy(j,11) * w(i,11,k,e) &
481 + dy(j,12) * w(i,12,k,e) &
482 + dy(j,13) * w(i,13,k,e) &
483 + dy(j,14) * w(i,14,k,e)
484 end do
485 end do
486 end do
487
488 do k = 1, lx
489 do i = 1, lx*lx
490 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
491 + dz(k,2) * u(i,1,2,e) &
492 + dz(k,3) * u(i,1,3,e) &
493 + dz(k,4) * u(i,1,4,e) &
494 + dz(k,5) * u(i,1,5,e) &
495 + dz(k,6) * u(i,1,6,e) &
496 + dz(k,7) * u(i,1,7,e) &
497 + dz(k,8) * u(i,1,8,e) &
498 + dz(k,9) * u(i,1,9,e) &
499 + dz(k,10) * u(i,1,10,e) &
500 + dz(k,11) * u(i,1,11,e) &
501 + dz(k,12) * u(i,1,12,e) &
502 + dz(k,13) * u(i,1,13,e) &
503 + dz(k,14) * u(i,1,14,e)
504
505 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
506 + dz(k,2) * v(i,1,2,e) &
507 + dz(k,3) * v(i,1,3,e) &
508 + dz(k,4) * v(i,1,4,e) &
509 + dz(k,5) * v(i,1,5,e) &
510 + dz(k,6) * v(i,1,6,e) &
511 + dz(k,7) * v(i,1,7,e) &
512 + dz(k,8) * v(i,1,8,e) &
513 + dz(k,9) * v(i,1,9,e) &
514 + dz(k,10) * v(i,1,10,e) &
515 + dz(k,11) * v(i,1,11,e) &
516 + dz(k,12) * v(i,1,12,e) &
517 + dz(k,13) * v(i,1,13,e) &
518 + dz(k,14) * v(i,1,14,e)
519
520 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
521 + dz(k,2) * w(i,1,2,e) &
522 + dz(k,3) * w(i,1,3,e) &
523 + dz(k,4) * w(i,1,4,e) &
524 + dz(k,5) * w(i,1,5,e) &
525 + dz(k,6) * w(i,1,6,e) &
526 + dz(k,7) * w(i,1,7,e) &
527 + dz(k,8) * w(i,1,8,e) &
528 + dz(k,9) * w(i,1,9,e) &
529 + dz(k,10) * w(i,1,10,e) &
530 + dz(k,11) * w(i,1,11,e) &
531 + dz(k,12) * w(i,1,12,e) &
532 + dz(k,13) * w(i,1,13,e) &
533 + dz(k,14) * w(i,1,14,e)
534 end do
535 end do
536
537 do i = 1, lx*lx*lx
538 ur(i,1,1) = h1(i,1,1,e) &
539 * ( g11(i,1,1,e) * wur(i,1,1) &
540 + g12(i,1,1,e) * wus(i,1,1) &
541 + g13(i,1,1,e) * wut(i,1,1) )
542 us(i,1,1) = h1(i,1,1,e) &
543 * ( g12(i,1,1,e) * wur(i,1,1) &
544 + g22(i,1,1,e) * wus(i,1,1) &
545 + g23(i,1,1,e) * wut(i,1,1) )
546 ut(i,1,1) = h1(i,1,1,e) &
547 * ( g13(i,1,1,e) * wur(i,1,1) &
548 + g23(i,1,1,e) * wus(i,1,1) &
549 + g33(i,1,1,e) * wut(i,1,1) )
550
551 vr(i,1,1) = h1(i,1,1,e) &
552 * ( g11(i,1,1,e) * wvr(i,1,1) &
553 + g12(i,1,1,e) * wvs(i,1,1) &
554 + g13(i,1,1,e) * wvt(i,1,1) )
555 vs(i,1,1) = h1(i,1,1,e) &
556 * ( g12(i,1,1,e) * wvr(i,1,1) &
557 + g22(i,1,1,e) * wvs(i,1,1) &
558 + g23(i,1,1,e) * wvt(i,1,1) )
559 vt(i,1,1) = h1(i,1,1,e) &
560 * ( g13(i,1,1,e) * wvr(i,1,1) &
561 + g23(i,1,1,e) * wvs(i,1,1) &
562 + g33(i,1,1,e) * wvt(i,1,1) )
563
564 wr(i,1,1) = h1(i,1,1,e) &
565 * ( g11(i,1,1,e) * wwr(i,1,1) &
566 + g12(i,1,1,e) * wws(i,1,1) &
567 + g13(i,1,1,e) * wwt(i,1,1) )
568 ws(i,1,1) = h1(i,1,1,e) &
569 * ( g12(i,1,1,e) * wwr(i,1,1) &
570 + g22(i,1,1,e) * wws(i,1,1) &
571 + g23(i,1,1,e) * wwt(i,1,1) )
572 wt(i,1,1) = h1(i,1,1,e) &
573 * ( g13(i,1,1,e) * wwr(i,1,1) &
574 + g23(i,1,1,e) * wws(i,1,1) &
575 + g33(i,1,1,e) * wwt(i,1,1) )
576 end do
577
578 do j = 1, lx*lx
579 do i = 1, lx
580 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
581 + dxt(i,2) * ur(2,j,1) &
582 + dxt(i,3) * ur(3,j,1) &
583 + dxt(i,4) * ur(4,j,1) &
584 + dxt(i,5) * ur(5,j,1) &
585 + dxt(i,6) * ur(6,j,1) &
586 + dxt(i,7) * ur(7,j,1) &
587 + dxt(i,8) * ur(8,j,1) &
588 + dxt(i,9) * ur(9,j,1) &
589 + dxt(i,10) * ur(10,j,1) &
590 + dxt(i,11) * ur(11,j,1) &
591 + dxt(i,12) * ur(12,j,1) &
592 + dxt(i,13) * ur(13,j,1) &
593 + dxt(i,14) * ur(14,j,1)
594
595 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
596 + dxt(i,2) * vr(2,j,1) &
597 + dxt(i,3) * vr(3,j,1) &
598 + dxt(i,4) * vr(4,j,1) &
599 + dxt(i,5) * vr(5,j,1) &
600 + dxt(i,6) * vr(6,j,1) &
601 + dxt(i,7) * vr(7,j,1) &
602 + dxt(i,8) * vr(8,j,1) &
603 + dxt(i,9) * vr(9,j,1) &
604 + dxt(i,10) * vr(10,j,1) &
605 + dxt(i,11) * vr(11,j,1) &
606 + dxt(i,12) * vr(12,j,1) &
607 + dxt(i,13) * vr(13,j,1) &
608 + dxt(i,14) * vr(14,j,1)
609
610 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
611 + dxt(i,2) * wr(2,j,1) &
612 + dxt(i,3) * wr(3,j,1) &
613 + dxt(i,4) * wr(4,j,1) &
614 + dxt(i,5) * wr(5,j,1) &
615 + dxt(i,6) * wr(6,j,1) &
616 + dxt(i,7) * wr(7,j,1) &
617 + dxt(i,8) * wr(8,j,1) &
618 + dxt(i,9) * wr(9,j,1) &
619 + dxt(i,10) * wr(10,j,1) &
620 + dxt(i,11) * wr(11,j,1) &
621 + dxt(i,12) * wr(12,j,1) &
622 + dxt(i,13) * wr(13,j,1) &
623 + dxt(i,14) * wr(14,j,1)
624 end do
625 end do
626
627 do k = 1, lx
628 do j = 1, lx
629 do i = 1, lx
630 au(i,j,k,e) = au(i,j,k,e) &
631 + dyt(j,1) * us(i,1,k) &
632 + dyt(j,2) * us(i,2,k) &
633 + dyt(j,3) * us(i,3,k) &
634 + dyt(j,4) * us(i,4,k) &
635 + dyt(j,5) * us(i,5,k) &
636 + dyt(j,6) * us(i,6,k) &
637 + dyt(j,7) * us(i,7,k) &
638 + dyt(j,8) * us(i,8,k) &
639 + dyt(j,9) * us(i,9,k) &
640 + dyt(j,10) * us(i,10,k) &
641 + dyt(j,11) * us(i,11,k) &
642 + dyt(j,12) * us(i,12,k) &
643 + dyt(j,13) * us(i,13,k) &
644 + dyt(j,14) * us(i,14,k)
645
646 av(i,j,k,e) = av(i,j,k,e) &
647 + dyt(j,1) * vs(i,1,k) &
648 + dyt(j,2) * vs(i,2,k) &
649 + dyt(j,3) * vs(i,3,k) &
650 + dyt(j,4) * vs(i,4,k) &
651 + dyt(j,5) * vs(i,5,k) &
652 + dyt(j,6) * vs(i,6,k) &
653 + dyt(j,7) * vs(i,7,k) &
654 + dyt(j,8) * vs(i,8,k) &
655 + dyt(j,9) * vs(i,9,k) &
656 + dyt(j,10) * vs(i,10,k) &
657 + dyt(j,11) * vs(i,11,k) &
658 + dyt(j,12) * vs(i,12,k) &
659 + dyt(j,13) * vs(i,13,k) &
660 + dyt(j,14) * vs(i,14,k)
661
662 aw(i,j,k,e) = aw(i,j,k,e) &
663 + dyt(j,1) * ws(i,1,k) &
664 + dyt(j,2) * ws(i,2,k) &
665 + dyt(j,3) * ws(i,3,k) &
666 + dyt(j,4) * ws(i,4,k) &
667 + dyt(j,5) * ws(i,5,k) &
668 + dyt(j,6) * ws(i,6,k) &
669 + dyt(j,7) * ws(i,7,k) &
670 + dyt(j,8) * ws(i,8,k) &
671 + dyt(j,9) * ws(i,9,k) &
672 + dyt(j,10) * ws(i,10,k) &
673 + dyt(j,11) * ws(i,11,k) &
674 + dyt(j,12) * ws(i,12,k) &
675 + dyt(j,13) * ws(i,13,k) &
676 + dyt(j,14) * ws(i,14,k)
677 end do
678 end do
679 end do
680
681 do k = 1, lx
682 do i = 1, lx*lx
683 au(i,1,k,e) = au(i,1,k,e) &
684 + dzt(k,1) * ut(i,1,1) &
685 + dzt(k,2) * ut(i,1,2) &
686 + dzt(k,3) * ut(i,1,3) &
687 + dzt(k,4) * ut(i,1,4) &
688 + dzt(k,5) * ut(i,1,5) &
689 + dzt(k,6) * ut(i,1,6) &
690 + dzt(k,7) * ut(i,1,7) &
691 + dzt(k,8) * ut(i,1,8) &
692 + dzt(k,9) * ut(i,1,9) &
693 + dzt(k,10) * ut(i,1,10) &
694 + dzt(k,11) * ut(i,1,11) &
695 + dzt(k,12) * ut(i,1,12) &
696 + dzt(k,13) * ut(i,1,13) &
697 + dzt(k,14) * ut(i,1,14)
698
699 av(i,1,k,e) = av(i,1,k,e) &
700 + dzt(k,1) * vt(i,1,1) &
701 + dzt(k,2) * vt(i,1,2) &
702 + dzt(k,3) * vt(i,1,3) &
703 + dzt(k,4) * vt(i,1,4) &
704 + dzt(k,5) * vt(i,1,5) &
705 + dzt(k,6) * vt(i,1,6) &
706 + dzt(k,7) * vt(i,1,7) &
707 + dzt(k,8) * vt(i,1,8) &
708 + dzt(k,9) * vt(i,1,9) &
709 + dzt(k,10) * vt(i,1,10) &
710 + dzt(k,11) * vt(i,1,11) &
711 + dzt(k,12) * vt(i,1,12) &
712 + dzt(k,13) * vt(i,1,13) &
713 + dzt(k,14) * vt(i,1,14)
714
715 aw(i,1,k,e) = aw(i,1,k,e) &
716 + dzt(k,1) * wt(i,1,1) &
717 + dzt(k,2) * wt(i,1,2) &
718 + dzt(k,3) * wt(i,1,3) &
719 + dzt(k,4) * wt(i,1,4) &
720 + dzt(k,5) * wt(i,1,5) &
721 + dzt(k,6) * wt(i,1,6) &
722 + dzt(k,7) * wt(i,1,7) &
723 + dzt(k,8) * wt(i,1,8) &
724 + dzt(k,9) * wt(i,1,9) &
725 + dzt(k,10) * wt(i,1,10) &
726 + dzt(k,11) * wt(i,1,11) &
727 + dzt(k,12) * wt(i,1,12) &
728 + dzt(k,13) * wt(i,1,13) &
729 + dzt(k,14) * wt(i,1,14)
730 end do
731 end do
732
733 end do
734 !$omp end do
735 end subroutine ax_helm_vector_lx14
736
737 subroutine ax_helm_vector_lx13(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
738 h1, G11, G22, G33, G12, G13, G23, n)
739 integer, parameter :: lx = 13
740 integer, intent(in) :: n
741 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
742 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
743 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
744 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
745 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
746 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
747 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
748 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
749 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
750 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
751 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
752 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
753 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
754 real(kind=rp), intent(in) :: dx(lx, lx)
755 real(kind=rp), intent(in) :: dy(lx, lx)
756 real(kind=rp), intent(in) :: dz(lx, lx)
757 real(kind=rp), intent(in) :: dxt(lx, lx)
758 real(kind=rp), intent(in) :: dyt(lx, lx)
759 real(kind=rp), intent(in) :: dzt(lx, lx)
760 real(kind=rp) :: ur(lx, lx, lx)
761 real(kind=rp) :: us(lx, lx, lx)
762 real(kind=rp) :: ut(lx, lx, lx)
763 real(kind=rp) :: vr(lx, lx, lx)
764 real(kind=rp) :: vs(lx, lx, lx)
765 real(kind=rp) :: vt(lx, lx, lx)
766 real(kind=rp) :: wr(lx, lx, lx)
767 real(kind=rp) :: ws(lx, lx, lx)
768 real(kind=rp) :: wt(lx, lx, lx)
769 real(kind=rp) :: wur(lx, lx, lx)
770 real(kind=rp) :: wus(lx, lx, lx)
771 real(kind=rp) :: wut(lx, lx, lx)
772 real(kind=rp) :: wvr(lx, lx, lx)
773 real(kind=rp) :: wvs(lx, lx, lx)
774 real(kind=rp) :: wvt(lx, lx, lx)
775 real(kind=rp) :: wwr(lx, lx, lx)
776 real(kind=rp) :: wws(lx, lx, lx)
777 real(kind=rp) :: wwt(lx, lx, lx)
778 integer :: e, i, j, k
779
780 !$omp do
781 do e = 1, n
782 do j = 1, lx * lx
783 do i = 1, lx
784 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
785 + dx(i,2) * u(2,j,1,e) &
786 + dx(i,3) * u(3,j,1,e) &
787 + dx(i,4) * u(4,j,1,e) &
788 + dx(i,5) * u(5,j,1,e) &
789 + dx(i,6) * u(6,j,1,e) &
790 + dx(i,7) * u(7,j,1,e) &
791 + dx(i,8) * u(8,j,1,e) &
792 + dx(i,9) * u(9,j,1,e) &
793 + dx(i,10) * u(10,j,1,e) &
794 + dx(i,11) * u(11,j,1,e) &
795 + dx(i,12) * u(12,j,1,e) &
796 + dx(i,13) * u(13,j,1,e)
797
798 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
799 + dx(i,2) * v(2,j,1,e) &
800 + dx(i,3) * v(3,j,1,e) &
801 + dx(i,4) * v(4,j,1,e) &
802 + dx(i,5) * v(5,j,1,e) &
803 + dx(i,6) * v(6,j,1,e) &
804 + dx(i,7) * v(7,j,1,e) &
805 + dx(i,8) * v(8,j,1,e) &
806 + dx(i,9) * v(9,j,1,e) &
807 + dx(i,10) * v(10,j,1,e) &
808 + dx(i,11) * v(11,j,1,e) &
809 + dx(i,12) * v(12,j,1,e) &
810 + dx(i,13) * v(13,j,1,e)
811
812 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
813 + dx(i,2) * w(2,j,1,e) &
814 + dx(i,3) * w(3,j,1,e) &
815 + dx(i,4) * w(4,j,1,e) &
816 + dx(i,5) * w(5,j,1,e) &
817 + dx(i,6) * w(6,j,1,e) &
818 + dx(i,7) * w(7,j,1,e) &
819 + dx(i,8) * w(8,j,1,e) &
820 + dx(i,9) * w(9,j,1,e) &
821 + dx(i,10) * w(10,j,1,e) &
822 + dx(i,11) * w(11,j,1,e) &
823 + dx(i,12) * w(12,j,1,e) &
824 + dx(i,13) * w(13,j,1,e)
825 end do
826 end do
827
828 do k = 1, lx
829 do j = 1, lx
830 do i = 1, lx
831 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
832 + dy(j,2) * u(i,2,k,e) &
833 + dy(j,3) * u(i,3,k,e) &
834 + dy(j,4) * u(i,4,k,e) &
835 + dy(j,5) * u(i,5,k,e) &
836 + dy(j,6) * u(i,6,k,e) &
837 + dy(j,7) * u(i,7,k,e) &
838 + dy(j,8) * u(i,8,k,e) &
839 + dy(j,9) * u(i,9,k,e) &
840 + dy(j,10) * u(i,10,k,e) &
841 + dy(j,11) * u(i,11,k,e) &
842 + dy(j,12) * u(i,12,k,e) &
843 + dy(j,13) * u(i,13,k,e)
844
845 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
846 + dy(j,2) * v(i,2,k,e) &
847 + dy(j,3) * v(i,3,k,e) &
848 + dy(j,4) * v(i,4,k,e) &
849 + dy(j,5) * v(i,5,k,e) &
850 + dy(j,6) * v(i,6,k,e) &
851 + dy(j,7) * v(i,7,k,e) &
852 + dy(j,8) * v(i,8,k,e) &
853 + dy(j,9) * v(i,9,k,e) &
854 + dy(j,10) * v(i,10,k,e) &
855 + dy(j,11) * v(i,11,k,e) &
856 + dy(j,12) * v(i,12,k,e) &
857 + dy(j,13) * v(i,13,k,e)
858
859 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
860 + dy(j,2) * w(i,2,k,e) &
861 + dy(j,3) * w(i,3,k,e) &
862 + dy(j,4) * w(i,4,k,e) &
863 + dy(j,5) * w(i,5,k,e) &
864 + dy(j,6) * w(i,6,k,e) &
865 + dy(j,7) * w(i,7,k,e) &
866 + dy(j,8) * w(i,8,k,e) &
867 + dy(j,9) * w(i,9,k,e) &
868 + dy(j,10) * w(i,10,k,e) &
869 + dy(j,11) * w(i,11,k,e) &
870 + dy(j,12) * w(i,12,k,e) &
871 + dy(j,13) * w(i,13,k,e)
872 end do
873 end do
874 end do
875
876 do k = 1, lx
877 do i = 1, lx*lx
878 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
879 + dz(k,2) * u(i,1,2,e) &
880 + dz(k,3) * u(i,1,3,e) &
881 + dz(k,4) * u(i,1,4,e) &
882 + dz(k,5) * u(i,1,5,e) &
883 + dz(k,6) * u(i,1,6,e) &
884 + dz(k,7) * u(i,1,7,e) &
885 + dz(k,8) * u(i,1,8,e) &
886 + dz(k,9) * u(i,1,9,e) &
887 + dz(k,10) * u(i,1,10,e) &
888 + dz(k,11) * u(i,1,11,e) &
889 + dz(k,12) * u(i,1,12,e) &
890 + dz(k,13) * u(i,1,13,e)
891
892 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
893 + dz(k,2) * v(i,1,2,e) &
894 + dz(k,3) * v(i,1,3,e) &
895 + dz(k,4) * v(i,1,4,e) &
896 + dz(k,5) * v(i,1,5,e) &
897 + dz(k,6) * v(i,1,6,e) &
898 + dz(k,7) * v(i,1,7,e) &
899 + dz(k,8) * v(i,1,8,e) &
900 + dz(k,9) * v(i,1,9,e) &
901 + dz(k,10) * v(i,1,10,e) &
902 + dz(k,11) * v(i,1,11,e) &
903 + dz(k,12) * v(i,1,12,e) &
904 + dz(k,13) * v(i,1,13,e)
905
906 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
907 + dz(k,2) * w(i,1,2,e) &
908 + dz(k,3) * w(i,1,3,e) &
909 + dz(k,4) * w(i,1,4,e) &
910 + dz(k,5) * w(i,1,5,e) &
911 + dz(k,6) * w(i,1,6,e) &
912 + dz(k,7) * w(i,1,7,e) &
913 + dz(k,8) * w(i,1,8,e) &
914 + dz(k,9) * w(i,1,9,e) &
915 + dz(k,10) * w(i,1,10,e) &
916 + dz(k,11) * w(i,1,11,e) &
917 + dz(k,12) * w(i,1,12,e) &
918 + dz(k,13) * w(i,1,13,e)
919 end do
920 end do
921
922 do i = 1, lx*lx*lx
923 ur(i,1,1) = h1(i,1,1,e) &
924 * ( g11(i,1,1,e) * wur(i,1,1) &
925 + g12(i,1,1,e) * wus(i,1,1) &
926 + g13(i,1,1,e) * wut(i,1,1) )
927 us(i,1,1) = h1(i,1,1,e) &
928 * ( g12(i,1,1,e) * wur(i,1,1) &
929 + g22(i,1,1,e) * wus(i,1,1) &
930 + g23(i,1,1,e) * wut(i,1,1) )
931 ut(i,1,1) = h1(i,1,1,e) &
932 * ( g13(i,1,1,e) * wur(i,1,1) &
933 + g23(i,1,1,e) * wus(i,1,1) &
934 + g33(i,1,1,e) * wut(i,1,1) )
935
936 vr(i,1,1) = h1(i,1,1,e) &
937 * ( g11(i,1,1,e) * wvr(i,1,1) &
938 + g12(i,1,1,e) * wvs(i,1,1) &
939 + g13(i,1,1,e) * wvt(i,1,1) )
940 vs(i,1,1) = h1(i,1,1,e) &
941 * ( g12(i,1,1,e) * wvr(i,1,1) &
942 + g22(i,1,1,e) * wvs(i,1,1) &
943 + g23(i,1,1,e) * wvt(i,1,1) )
944 vt(i,1,1) = h1(i,1,1,e) &
945 * ( g13(i,1,1,e) * wvr(i,1,1) &
946 + g23(i,1,1,e) * wvs(i,1,1) &
947 + g33(i,1,1,e) * wvt(i,1,1) )
948
949 wr(i,1,1) = h1(i,1,1,e) &
950 * ( g11(i,1,1,e) * wwr(i,1,1) &
951 + g12(i,1,1,e) * wws(i,1,1) &
952 + g13(i,1,1,e) * wwt(i,1,1) )
953 ws(i,1,1) = h1(i,1,1,e) &
954 * ( g12(i,1,1,e) * wwr(i,1,1) &
955 + g22(i,1,1,e) * wws(i,1,1) &
956 + g23(i,1,1,e) * wwt(i,1,1) )
957 wt(i,1,1) = h1(i,1,1,e) &
958 * ( g13(i,1,1,e) * wwr(i,1,1) &
959 + g23(i,1,1,e) * wws(i,1,1) &
960 + g33(i,1,1,e) * wwt(i,1,1) )
961 end do
962
963 do j = 1, lx*lx
964 do i = 1, lx
965 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
966 + dxt(i,2) * ur(2,j,1) &
967 + dxt(i,3) * ur(3,j,1) &
968 + dxt(i,4) * ur(4,j,1) &
969 + dxt(i,5) * ur(5,j,1) &
970 + dxt(i,6) * ur(6,j,1) &
971 + dxt(i,7) * ur(7,j,1) &
972 + dxt(i,8) * ur(8,j,1) &
973 + dxt(i,9) * ur(9,j,1) &
974 + dxt(i,10) * ur(10,j,1) &
975 + dxt(i,11) * ur(11,j,1) &
976 + dxt(i,12) * ur(12,j,1) &
977 + dxt(i,13) * ur(13,j,1)
978
979 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
980 + dxt(i,2) * vr(2,j,1) &
981 + dxt(i,3) * vr(3,j,1) &
982 + dxt(i,4) * vr(4,j,1) &
983 + dxt(i,5) * vr(5,j,1) &
984 + dxt(i,6) * vr(6,j,1) &
985 + dxt(i,7) * vr(7,j,1) &
986 + dxt(i,8) * vr(8,j,1) &
987 + dxt(i,9) * vr(9,j,1) &
988 + dxt(i,10) * vr(10,j,1) &
989 + dxt(i,11) * vr(11,j,1) &
990 + dxt(i,12) * vr(12,j,1) &
991 + dxt(i,13) * vr(13,j,1)
992
993 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
994 + dxt(i,2) * wr(2,j,1) &
995 + dxt(i,3) * wr(3,j,1) &
996 + dxt(i,4) * wr(4,j,1) &
997 + dxt(i,5) * wr(5,j,1) &
998 + dxt(i,6) * wr(6,j,1) &
999 + dxt(i,7) * wr(7,j,1) &
1000 + dxt(i,8) * wr(8,j,1) &
1001 + dxt(i,9) * wr(9,j,1) &
1002 + dxt(i,10) * wr(10,j,1) &
1003 + dxt(i,11) * wr(11,j,1) &
1004 + dxt(i,12) * wr(12,j,1) &
1005 + dxt(i,13) * wr(13,j,1)
1006 end do
1007 end do
1008
1009 do k = 1, lx
1010 do j = 1, lx
1011 do i = 1, lx
1012 au(i,j,k,e) = au(i,j,k,e) &
1013 + dyt(j,1) * us(i,1,k) &
1014 + dyt(j,2) * us(i,2,k) &
1015 + dyt(j,3) * us(i,3,k) &
1016 + dyt(j,4) * us(i,4,k) &
1017 + dyt(j,5) * us(i,5,k) &
1018 + dyt(j,6) * us(i,6,k) &
1019 + dyt(j,7) * us(i,7,k) &
1020 + dyt(j,8) * us(i,8,k) &
1021 + dyt(j,9) * us(i,9,k) &
1022 + dyt(j,10) * us(i,10,k) &
1023 + dyt(j,11) * us(i,11,k) &
1024 + dyt(j,12) * us(i,12,k) &
1025 + dyt(j,13) * us(i,13,k)
1026
1027 av(i,j,k,e) = av(i,j,k,e) &
1028 + dyt(j,1) * vs(i,1,k) &
1029 + dyt(j,2) * vs(i,2,k) &
1030 + dyt(j,3) * vs(i,3,k) &
1031 + dyt(j,4) * vs(i,4,k) &
1032 + dyt(j,5) * vs(i,5,k) &
1033 + dyt(j,6) * vs(i,6,k) &
1034 + dyt(j,7) * vs(i,7,k) &
1035 + dyt(j,8) * vs(i,8,k) &
1036 + dyt(j,9) * vs(i,9,k) &
1037 + dyt(j,10) * vs(i,10,k) &
1038 + dyt(j,11) * vs(i,11,k) &
1039 + dyt(j,12) * vs(i,12,k) &
1040 + dyt(j,13) * vs(i,13,k)
1041
1042 aw(i,j,k,e) = aw(i,j,k,e) &
1043 + dyt(j,1) * ws(i,1,k) &
1044 + dyt(j,2) * ws(i,2,k) &
1045 + dyt(j,3) * ws(i,3,k) &
1046 + dyt(j,4) * ws(i,4,k) &
1047 + dyt(j,5) * ws(i,5,k) &
1048 + dyt(j,6) * ws(i,6,k) &
1049 + dyt(j,7) * ws(i,7,k) &
1050 + dyt(j,8) * ws(i,8,k) &
1051 + dyt(j,9) * ws(i,9,k) &
1052 + dyt(j,10) * ws(i,10,k) &
1053 + dyt(j,11) * ws(i,11,k) &
1054 + dyt(j,12) * ws(i,12,k) &
1055 + dyt(j,13) * ws(i,13,k)
1056 end do
1057 end do
1058 end do
1059
1060 do k = 1, lx
1061 do i = 1, lx*lx
1062 au(i,1,k,e) = au(i,1,k,e) &
1063 + dzt(k,1) * ut(i,1,1) &
1064 + dzt(k,2) * ut(i,1,2) &
1065 + dzt(k,3) * ut(i,1,3) &
1066 + dzt(k,4) * ut(i,1,4) &
1067 + dzt(k,5) * ut(i,1,5) &
1068 + dzt(k,6) * ut(i,1,6) &
1069 + dzt(k,7) * ut(i,1,7) &
1070 + dzt(k,8) * ut(i,1,8) &
1071 + dzt(k,9) * ut(i,1,9) &
1072 + dzt(k,10) * ut(i,1,10) &
1073 + dzt(k,11) * ut(i,1,11) &
1074 + dzt(k,12) * ut(i,1,12) &
1075 + dzt(k,13) * ut(i,1,13)
1076
1077 av(i,1,k,e) = av(i,1,k,e) &
1078 + dzt(k,1) * vt(i,1,1) &
1079 + dzt(k,2) * vt(i,1,2) &
1080 + dzt(k,3) * vt(i,1,3) &
1081 + dzt(k,4) * vt(i,1,4) &
1082 + dzt(k,5) * vt(i,1,5) &
1083 + dzt(k,6) * vt(i,1,6) &
1084 + dzt(k,7) * vt(i,1,7) &
1085 + dzt(k,8) * vt(i,1,8) &
1086 + dzt(k,9) * vt(i,1,9) &
1087 + dzt(k,10) * vt(i,1,10) &
1088 + dzt(k,11) * vt(i,1,11) &
1089 + dzt(k,12) * vt(i,1,12) &
1090 + dzt(k,13) * vt(i,1,13)
1091
1092 aw(i,1,k,e) = aw(i,1,k,e) &
1093 + dzt(k,1) * wt(i,1,1) &
1094 + dzt(k,2) * wt(i,1,2) &
1095 + dzt(k,3) * wt(i,1,3) &
1096 + dzt(k,4) * wt(i,1,4) &
1097 + dzt(k,5) * wt(i,1,5) &
1098 + dzt(k,6) * wt(i,1,6) &
1099 + dzt(k,7) * wt(i,1,7) &
1100 + dzt(k,8) * wt(i,1,8) &
1101 + dzt(k,9) * wt(i,1,9) &
1102 + dzt(k,10) * wt(i,1,10) &
1103 + dzt(k,11) * wt(i,1,11) &
1104 + dzt(k,12) * wt(i,1,12) &
1105 + dzt(k,13) * wt(i,1,13)
1106 end do
1107 end do
1108
1109 end do
1110 !$omp end do
1111 end subroutine ax_helm_vector_lx13
1112
1113 subroutine ax_helm_vector_lx12(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1114 h1, G11, G22, G33, G12, G13, G23, n)
1115 integer, parameter :: lx = 12
1116 integer, intent(in) :: n
1117 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
1118 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
1119 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
1120 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1121 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
1122 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
1123 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1124 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1125 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1126 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1127 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1128 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1129 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1130 real(kind=rp), intent(in) :: dx(lx, lx)
1131 real(kind=rp), intent(in) :: dy(lx, lx)
1132 real(kind=rp), intent(in) :: dz(lx, lx)
1133 real(kind=rp), intent(in) :: dxt(lx, lx)
1134 real(kind=rp), intent(in) :: dyt(lx, lx)
1135 real(kind=rp), intent(in) :: dzt(lx, lx)
1136 real(kind=rp) :: ur(lx, lx, lx)
1137 real(kind=rp) :: us(lx, lx, lx)
1138 real(kind=rp) :: ut(lx, lx, lx)
1139 real(kind=rp) :: vr(lx, lx, lx)
1140 real(kind=rp) :: vs(lx, lx, lx)
1141 real(kind=rp) :: vt(lx, lx, lx)
1142 real(kind=rp) :: wr(lx, lx, lx)
1143 real(kind=rp) :: ws(lx, lx, lx)
1144 real(kind=rp) :: wt(lx, lx, lx)
1145 real(kind=rp) :: wur(lx, lx, lx)
1146 real(kind=rp) :: wus(lx, lx, lx)
1147 real(kind=rp) :: wut(lx, lx, lx)
1148 real(kind=rp) :: wvr(lx, lx, lx)
1149 real(kind=rp) :: wvs(lx, lx, lx)
1150 real(kind=rp) :: wvt(lx, lx, lx)
1151 real(kind=rp) :: wwr(lx, lx, lx)
1152 real(kind=rp) :: wws(lx, lx, lx)
1153 real(kind=rp) :: wwt(lx, lx, lx)
1154 integer :: e, i, j, k
1155
1156 !$omp do
1157 do e = 1, n
1158 do j = 1, lx * lx
1159 do i = 1, lx
1160 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1161 + dx(i,2) * u(2,j,1,e) &
1162 + dx(i,3) * u(3,j,1,e) &
1163 + dx(i,4) * u(4,j,1,e) &
1164 + dx(i,5) * u(5,j,1,e) &
1165 + dx(i,6) * u(6,j,1,e) &
1166 + dx(i,7) * u(7,j,1,e) &
1167 + dx(i,8) * u(8,j,1,e) &
1168 + dx(i,9) * u(9,j,1,e) &
1169 + dx(i,10) * u(10,j,1,e) &
1170 + dx(i,11) * u(11,j,1,e) &
1171 + dx(i,12) * u(12,j,1,e)
1172
1173 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
1174 + dx(i,2) * v(2,j,1,e) &
1175 + dx(i,3) * v(3,j,1,e) &
1176 + dx(i,4) * v(4,j,1,e) &
1177 + dx(i,5) * v(5,j,1,e) &
1178 + dx(i,6) * v(6,j,1,e) &
1179 + dx(i,7) * v(7,j,1,e) &
1180 + dx(i,8) * v(8,j,1,e) &
1181 + dx(i,9) * v(9,j,1,e) &
1182 + dx(i,10) * v(10,j,1,e) &
1183 + dx(i,11) * v(11,j,1,e) &
1184 + dx(i,12) * v(12,j,1,e)
1185
1186 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
1187 + dx(i,2) * w(2,j,1,e) &
1188 + dx(i,3) * w(3,j,1,e) &
1189 + dx(i,4) * w(4,j,1,e) &
1190 + dx(i,5) * w(5,j,1,e) &
1191 + dx(i,6) * w(6,j,1,e) &
1192 + dx(i,7) * w(7,j,1,e) &
1193 + dx(i,8) * w(8,j,1,e) &
1194 + dx(i,9) * w(9,j,1,e) &
1195 + dx(i,10) * w(10,j,1,e) &
1196 + dx(i,11) * w(11,j,1,e) &
1197 + dx(i,12) * w(12,j,1,e)
1198 end do
1199 end do
1200
1201 do k = 1, lx
1202 do j = 1, lx
1203 do i = 1, lx
1204 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1205 + dy(j,2) * u(i,2,k,e) &
1206 + dy(j,3) * u(i,3,k,e) &
1207 + dy(j,4) * u(i,4,k,e) &
1208 + dy(j,5) * u(i,5,k,e) &
1209 + dy(j,6) * u(i,6,k,e) &
1210 + dy(j,7) * u(i,7,k,e) &
1211 + dy(j,8) * u(i,8,k,e) &
1212 + dy(j,9) * u(i,9,k,e) &
1213 + dy(j,10) * u(i,10,k,e) &
1214 + dy(j,11) * u(i,11,k,e) &
1215 + dy(j,12) * u(i,12,k,e)
1216
1217 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
1218 + dy(j,2) * v(i,2,k,e) &
1219 + dy(j,3) * v(i,3,k,e) &
1220 + dy(j,4) * v(i,4,k,e) &
1221 + dy(j,5) * v(i,5,k,e) &
1222 + dy(j,6) * v(i,6,k,e) &
1223 + dy(j,7) * v(i,7,k,e) &
1224 + dy(j,8) * v(i,8,k,e) &
1225 + dy(j,9) * v(i,9,k,e) &
1226 + dy(j,10) * v(i,10,k,e) &
1227 + dy(j,11) * v(i,11,k,e) &
1228 + dy(j,12) * v(i,12,k,e)
1229
1230 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
1231 + dy(j,2) * w(i,2,k,e) &
1232 + dy(j,3) * w(i,3,k,e) &
1233 + dy(j,4) * w(i,4,k,e) &
1234 + dy(j,5) * w(i,5,k,e) &
1235 + dy(j,6) * w(i,6,k,e) &
1236 + dy(j,7) * w(i,7,k,e) &
1237 + dy(j,8) * w(i,8,k,e) &
1238 + dy(j,9) * w(i,9,k,e) &
1239 + dy(j,10) * w(i,10,k,e) &
1240 + dy(j,11) * w(i,11,k,e) &
1241 + dy(j,12) * w(i,12,k,e)
1242 end do
1243 end do
1244 end do
1245
1246 do k = 1, lx
1247 do i = 1, lx*lx
1248 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1249 + dz(k,2) * u(i,1,2,e) &
1250 + dz(k,3) * u(i,1,3,e) &
1251 + dz(k,4) * u(i,1,4,e) &
1252 + dz(k,5) * u(i,1,5,e) &
1253 + dz(k,6) * u(i,1,6,e) &
1254 + dz(k,7) * u(i,1,7,e) &
1255 + dz(k,8) * u(i,1,8,e) &
1256 + dz(k,9) * u(i,1,9,e) &
1257 + dz(k,10) * u(i,1,10,e) &
1258 + dz(k,11) * u(i,1,11,e) &
1259 + dz(k,12) * u(i,1,12,e)
1260
1261 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
1262 + dz(k,2) * v(i,1,2,e) &
1263 + dz(k,3) * v(i,1,3,e) &
1264 + dz(k,4) * v(i,1,4,e) &
1265 + dz(k,5) * v(i,1,5,e) &
1266 + dz(k,6) * v(i,1,6,e) &
1267 + dz(k,7) * v(i,1,7,e) &
1268 + dz(k,8) * v(i,1,8,e) &
1269 + dz(k,9) * v(i,1,9,e) &
1270 + dz(k,10) * v(i,1,10,e) &
1271 + dz(k,11) * v(i,1,11,e) &
1272 + dz(k,12) * v(i,1,12,e)
1273
1274 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
1275 + dz(k,2) * w(i,1,2,e) &
1276 + dz(k,3) * w(i,1,3,e) &
1277 + dz(k,4) * w(i,1,4,e) &
1278 + dz(k,5) * w(i,1,5,e) &
1279 + dz(k,6) * w(i,1,6,e) &
1280 + dz(k,7) * w(i,1,7,e) &
1281 + dz(k,8) * w(i,1,8,e) &
1282 + dz(k,9) * w(i,1,9,e) &
1283 + dz(k,10) * w(i,1,10,e) &
1284 + dz(k,11) * w(i,1,11,e) &
1285 + dz(k,12) * w(i,1,12,e)
1286 end do
1287 end do
1288
1289 do i = 1, lx*lx*lx
1290 ur(i,1,1) = h1(i,1,1,e) &
1291 * ( g11(i,1,1,e) * wur(i,1,1) &
1292 + g12(i,1,1,e) * wus(i,1,1) &
1293 + g13(i,1,1,e) * wut(i,1,1) )
1294 us(i,1,1) = h1(i,1,1,e) &
1295 * ( g12(i,1,1,e) * wur(i,1,1) &
1296 + g22(i,1,1,e) * wus(i,1,1) &
1297 + g23(i,1,1,e) * wut(i,1,1) )
1298 ut(i,1,1) = h1(i,1,1,e) &
1299 * ( g13(i,1,1,e) * wur(i,1,1) &
1300 + g23(i,1,1,e) * wus(i,1,1) &
1301 + g33(i,1,1,e) * wut(i,1,1) )
1302
1303 vr(i,1,1) = h1(i,1,1,e) &
1304 * ( g11(i,1,1,e) * wvr(i,1,1) &
1305 + g12(i,1,1,e) * wvs(i,1,1) &
1306 + g13(i,1,1,e) * wvt(i,1,1) )
1307 vs(i,1,1) = h1(i,1,1,e) &
1308 * ( g12(i,1,1,e) * wvr(i,1,1) &
1309 + g22(i,1,1,e) * wvs(i,1,1) &
1310 + g23(i,1,1,e) * wvt(i,1,1) )
1311 vt(i,1,1) = h1(i,1,1,e) &
1312 * ( g13(i,1,1,e) * wvr(i,1,1) &
1313 + g23(i,1,1,e) * wvs(i,1,1) &
1314 + g33(i,1,1,e) * wvt(i,1,1) )
1315
1316 wr(i,1,1) = h1(i,1,1,e) &
1317 * ( g11(i,1,1,e) * wwr(i,1,1) &
1318 + g12(i,1,1,e) * wws(i,1,1) &
1319 + g13(i,1,1,e) * wwt(i,1,1) )
1320 ws(i,1,1) = h1(i,1,1,e) &
1321 * ( g12(i,1,1,e) * wwr(i,1,1) &
1322 + g22(i,1,1,e) * wws(i,1,1) &
1323 + g23(i,1,1,e) * wwt(i,1,1) )
1324 wt(i,1,1) = h1(i,1,1,e) &
1325 * ( g13(i,1,1,e) * wwr(i,1,1) &
1326 + g23(i,1,1,e) * wws(i,1,1) &
1327 + g33(i,1,1,e) * wwt(i,1,1) )
1328 end do
1329
1330 do j = 1, lx*lx
1331 do i = 1, lx
1332 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1333 + dxt(i,2) * ur(2,j,1) &
1334 + dxt(i,3) * ur(3,j,1) &
1335 + dxt(i,4) * ur(4,j,1) &
1336 + dxt(i,5) * ur(5,j,1) &
1337 + dxt(i,6) * ur(6,j,1) &
1338 + dxt(i,7) * ur(7,j,1) &
1339 + dxt(i,8) * ur(8,j,1) &
1340 + dxt(i,9) * ur(9,j,1) &
1341 + dxt(i,10) * ur(10,j,1) &
1342 + dxt(i,11) * ur(11,j,1) &
1343 + dxt(i,12) * ur(12,j,1)
1344
1345 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
1346 + dxt(i,2) * vr(2,j,1) &
1347 + dxt(i,3) * vr(3,j,1) &
1348 + dxt(i,4) * vr(4,j,1) &
1349 + dxt(i,5) * vr(5,j,1) &
1350 + dxt(i,6) * vr(6,j,1) &
1351 + dxt(i,7) * vr(7,j,1) &
1352 + dxt(i,8) * vr(8,j,1) &
1353 + dxt(i,9) * vr(9,j,1) &
1354 + dxt(i,10) * vr(10,j,1) &
1355 + dxt(i,11) * vr(11,j,1) &
1356 + dxt(i,12) * vr(12,j,1)
1357
1358 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
1359 + dxt(i,2) * wr(2,j,1) &
1360 + dxt(i,3) * wr(3,j,1) &
1361 + dxt(i,4) * wr(4,j,1) &
1362 + dxt(i,5) * wr(5,j,1) &
1363 + dxt(i,6) * wr(6,j,1) &
1364 + dxt(i,7) * wr(7,j,1) &
1365 + dxt(i,8) * wr(8,j,1) &
1366 + dxt(i,9) * wr(9,j,1) &
1367 + dxt(i,10) * wr(10,j,1) &
1368 + dxt(i,11) * wr(11,j,1) &
1369 + dxt(i,12) * wr(12,j,1)
1370 end do
1371 end do
1372
1373 do k = 1, lx
1374 do j = 1, lx
1375 do i = 1, lx
1376 au(i,j,k,e) = au(i,j,k,e) &
1377 + dyt(j,1) * us(i,1,k) &
1378 + dyt(j,2) * us(i,2,k) &
1379 + dyt(j,3) * us(i,3,k) &
1380 + dyt(j,4) * us(i,4,k) &
1381 + dyt(j,5) * us(i,5,k) &
1382 + dyt(j,6) * us(i,6,k) &
1383 + dyt(j,7) * us(i,7,k) &
1384 + dyt(j,8) * us(i,8,k) &
1385 + dyt(j,9) * us(i,9,k) &
1386 + dyt(j,10) * us(i,10,k) &
1387 + dyt(j,11) * us(i,11,k) &
1388 + dyt(j,12) * us(i,12,k)
1389
1390 av(i,j,k,e) = av(i,j,k,e) &
1391 + dyt(j,1) * vs(i,1,k) &
1392 + dyt(j,2) * vs(i,2,k) &
1393 + dyt(j,3) * vs(i,3,k) &
1394 + dyt(j,4) * vs(i,4,k) &
1395 + dyt(j,5) * vs(i,5,k) &
1396 + dyt(j,6) * vs(i,6,k) &
1397 + dyt(j,7) * vs(i,7,k) &
1398 + dyt(j,8) * vs(i,8,k) &
1399 + dyt(j,9) * vs(i,9,k) &
1400 + dyt(j,10) * vs(i,10,k) &
1401 + dyt(j,11) * vs(i,11,k) &
1402 + dyt(j,12) * vs(i,12,k)
1403
1404 aw(i,j,k,e) = aw(i,j,k,e) &
1405 + dyt(j,1) * ws(i,1,k) &
1406 + dyt(j,2) * ws(i,2,k) &
1407 + dyt(j,3) * ws(i,3,k) &
1408 + dyt(j,4) * ws(i,4,k) &
1409 + dyt(j,5) * ws(i,5,k) &
1410 + dyt(j,6) * ws(i,6,k) &
1411 + dyt(j,7) * ws(i,7,k) &
1412 + dyt(j,8) * ws(i,8,k) &
1413 + dyt(j,9) * ws(i,9,k) &
1414 + dyt(j,10) * ws(i,10,k) &
1415 + dyt(j,11) * ws(i,11,k) &
1416 + dyt(j,12) * ws(i,12,k)
1417 end do
1418 end do
1419 end do
1420
1421 do k = 1, lx
1422 do i = 1, lx*lx
1423 au(i,1,k,e) = au(i,1,k,e) &
1424 + dzt(k,1) * ut(i,1,1) &
1425 + dzt(k,2) * ut(i,1,2) &
1426 + dzt(k,3) * ut(i,1,3) &
1427 + dzt(k,4) * ut(i,1,4) &
1428 + dzt(k,5) * ut(i,1,5) &
1429 + dzt(k,6) * ut(i,1,6) &
1430 + dzt(k,7) * ut(i,1,7) &
1431 + dzt(k,8) * ut(i,1,8) &
1432 + dzt(k,9) * ut(i,1,9) &
1433 + dzt(k,10) * ut(i,1,10) &
1434 + dzt(k,11) * ut(i,1,11) &
1435 + dzt(k,12) * ut(i,1,12)
1436
1437 av(i,1,k,e) = av(i,1,k,e) &
1438 + dzt(k,1) * vt(i,1,1) &
1439 + dzt(k,2) * vt(i,1,2) &
1440 + dzt(k,3) * vt(i,1,3) &
1441 + dzt(k,4) * vt(i,1,4) &
1442 + dzt(k,5) * vt(i,1,5) &
1443 + dzt(k,6) * vt(i,1,6) &
1444 + dzt(k,7) * vt(i,1,7) &
1445 + dzt(k,8) * vt(i,1,8) &
1446 + dzt(k,9) * vt(i,1,9) &
1447 + dzt(k,10) * vt(i,1,10) &
1448 + dzt(k,11) * vt(i,1,11) &
1449 + dzt(k,12) * vt(i,1,12)
1450
1451 aw(i,1,k,e) = aw(i,1,k,e) &
1452 + dzt(k,1) * wt(i,1,1) &
1453 + dzt(k,2) * wt(i,1,2) &
1454 + dzt(k,3) * wt(i,1,3) &
1455 + dzt(k,4) * wt(i,1,4) &
1456 + dzt(k,5) * wt(i,1,5) &
1457 + dzt(k,6) * wt(i,1,6) &
1458 + dzt(k,7) * wt(i,1,7) &
1459 + dzt(k,8) * wt(i,1,8) &
1460 + dzt(k,9) * wt(i,1,9) &
1461 + dzt(k,10) * wt(i,1,10) &
1462 + dzt(k,11) * wt(i,1,11) &
1463 + dzt(k,12) * wt(i,1,12)
1464 end do
1465 end do
1466
1467 end do
1468 !$omp end do
1469 end subroutine ax_helm_vector_lx12
1470
1471 subroutine ax_helm_vector_lx11(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1472 h1, G11, G22, G33, G12, G13, G23, n)
1473 integer, parameter :: lx = 11
1474 integer, intent(in) :: n
1475 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
1476 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
1477 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
1478 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1479 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
1480 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
1481 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1482 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1483 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1484 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1485 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1486 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1487 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1488 real(kind=rp), intent(in) :: dx(lx, lx)
1489 real(kind=rp), intent(in) :: dy(lx, lx)
1490 real(kind=rp), intent(in) :: dz(lx, lx)
1491 real(kind=rp), intent(in) :: dxt(lx, lx)
1492 real(kind=rp), intent(in) :: dyt(lx, lx)
1493 real(kind=rp), intent(in) :: dzt(lx, lx)
1494 real(kind=rp) :: ur(lx, lx, lx)
1495 real(kind=rp) :: us(lx, lx, lx)
1496 real(kind=rp) :: ut(lx, lx, lx)
1497 real(kind=rp) :: vr(lx, lx, lx)
1498 real(kind=rp) :: vs(lx, lx, lx)
1499 real(kind=rp) :: vt(lx, lx, lx)
1500 real(kind=rp) :: wr(lx, lx, lx)
1501 real(kind=rp) :: ws(lx, lx, lx)
1502 real(kind=rp) :: wt(lx, lx, lx)
1503 real(kind=rp) :: wur(lx, lx, lx)
1504 real(kind=rp) :: wus(lx, lx, lx)
1505 real(kind=rp) :: wut(lx, lx, lx)
1506 real(kind=rp) :: wvr(lx, lx, lx)
1507 real(kind=rp) :: wvs(lx, lx, lx)
1508 real(kind=rp) :: wvt(lx, lx, lx)
1509 real(kind=rp) :: wwr(lx, lx, lx)
1510 real(kind=rp) :: wws(lx, lx, lx)
1511 real(kind=rp) :: wwt(lx, lx, lx)
1512 integer :: e, i, j, k
1513
1514 !$omp do
1515 do e = 1, n
1516 do j = 1, lx * lx
1517 do i = 1, lx
1518 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1519 + dx(i,2) * u(2,j,1,e) &
1520 + dx(i,3) * u(3,j,1,e) &
1521 + dx(i,4) * u(4,j,1,e) &
1522 + dx(i,5) * u(5,j,1,e) &
1523 + dx(i,6) * u(6,j,1,e) &
1524 + dx(i,7) * u(7,j,1,e) &
1525 + dx(i,8) * u(8,j,1,e) &
1526 + dx(i,9) * u(9,j,1,e) &
1527 + dx(i,10) * u(10,j,1,e) &
1528 + dx(i,11) * u(11,j,1,e)
1529
1530 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
1531 + dx(i,2) * v(2,j,1,e) &
1532 + dx(i,3) * v(3,j,1,e) &
1533 + dx(i,4) * v(4,j,1,e) &
1534 + dx(i,5) * v(5,j,1,e) &
1535 + dx(i,6) * v(6,j,1,e) &
1536 + dx(i,7) * v(7,j,1,e) &
1537 + dx(i,8) * v(8,j,1,e) &
1538 + dx(i,9) * v(9,j,1,e) &
1539 + dx(i,10) * v(10,j,1,e) &
1540 + dx(i,11) * v(11,j,1,e)
1541
1542 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
1543 + dx(i,2) * w(2,j,1,e) &
1544 + dx(i,3) * w(3,j,1,e) &
1545 + dx(i,4) * w(4,j,1,e) &
1546 + dx(i,5) * w(5,j,1,e) &
1547 + dx(i,6) * w(6,j,1,e) &
1548 + dx(i,7) * w(7,j,1,e) &
1549 + dx(i,8) * w(8,j,1,e) &
1550 + dx(i,9) * w(9,j,1,e) &
1551 + dx(i,10) * w(10,j,1,e) &
1552 + dx(i,11) * w(11,j,1,e)
1553 end do
1554 end do
1555
1556 do k = 1, lx
1557 do j = 1, lx
1558 do i = 1, lx
1559 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1560 + dy(j,2) * u(i,2,k,e) &
1561 + dy(j,3) * u(i,3,k,e) &
1562 + dy(j,4) * u(i,4,k,e) &
1563 + dy(j,5) * u(i,5,k,e) &
1564 + dy(j,6) * u(i,6,k,e) &
1565 + dy(j,7) * u(i,7,k,e) &
1566 + dy(j,8) * u(i,8,k,e) &
1567 + dy(j,9) * u(i,9,k,e) &
1568 + dy(j,10) * u(i,10,k,e) &
1569 + dy(j,11) * u(i,11,k,e)
1570
1571 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
1572 + dy(j,2) * v(i,2,k,e) &
1573 + dy(j,3) * v(i,3,k,e) &
1574 + dy(j,4) * v(i,4,k,e) &
1575 + dy(j,5) * v(i,5,k,e) &
1576 + dy(j,6) * v(i,6,k,e) &
1577 + dy(j,7) * v(i,7,k,e) &
1578 + dy(j,8) * v(i,8,k,e) &
1579 + dy(j,9) * v(i,9,k,e) &
1580 + dy(j,10) * v(i,10,k,e) &
1581 + dy(j,11) * v(i,11,k,e)
1582
1583 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
1584 + dy(j,2) * w(i,2,k,e) &
1585 + dy(j,3) * w(i,3,k,e) &
1586 + dy(j,4) * w(i,4,k,e) &
1587 + dy(j,5) * w(i,5,k,e) &
1588 + dy(j,6) * w(i,6,k,e) &
1589 + dy(j,7) * w(i,7,k,e) &
1590 + dy(j,8) * w(i,8,k,e) &
1591 + dy(j,9) * w(i,9,k,e) &
1592 + dy(j,10) * w(i,10,k,e) &
1593 + dy(j,11) * w(i,11,k,e)
1594 end do
1595 end do
1596 end do
1597
1598 do k = 1, lx
1599 do i = 1, lx*lx
1600 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1601 + dz(k,2) * u(i,1,2,e) &
1602 + dz(k,3) * u(i,1,3,e) &
1603 + dz(k,4) * u(i,1,4,e) &
1604 + dz(k,5) * u(i,1,5,e) &
1605 + dz(k,6) * u(i,1,6,e) &
1606 + dz(k,7) * u(i,1,7,e) &
1607 + dz(k,8) * u(i,1,8,e) &
1608 + dz(k,9) * u(i,1,9,e) &
1609 + dz(k,10) * u(i,1,10,e) &
1610 + dz(k,11) * u(i,1,11,e)
1611
1612 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
1613 + dz(k,2) * v(i,1,2,e) &
1614 + dz(k,3) * v(i,1,3,e) &
1615 + dz(k,4) * v(i,1,4,e) &
1616 + dz(k,5) * v(i,1,5,e) &
1617 + dz(k,6) * v(i,1,6,e) &
1618 + dz(k,7) * v(i,1,7,e) &
1619 + dz(k,8) * v(i,1,8,e) &
1620 + dz(k,9) * v(i,1,9,e) &
1621 + dz(k,10) * v(i,1,10,e) &
1622 + dz(k,11) * v(i,1,11,e)
1623
1624 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
1625 + dz(k,2) * w(i,1,2,e) &
1626 + dz(k,3) * w(i,1,3,e) &
1627 + dz(k,4) * w(i,1,4,e) &
1628 + dz(k,5) * w(i,1,5,e) &
1629 + dz(k,6) * w(i,1,6,e) &
1630 + dz(k,7) * w(i,1,7,e) &
1631 + dz(k,8) * w(i,1,8,e) &
1632 + dz(k,9) * w(i,1,9,e) &
1633 + dz(k,10) * w(i,1,10,e) &
1634 + dz(k,11) * w(i,1,11,e)
1635 end do
1636 end do
1637
1638 do i = 1, lx*lx*lx
1639 ur(i,1,1) = h1(i,1,1,e) &
1640 * ( g11(i,1,1,e) * wur(i,1,1) &
1641 + g12(i,1,1,e) * wus(i,1,1) &
1642 + g13(i,1,1,e) * wut(i,1,1) )
1643 us(i,1,1) = h1(i,1,1,e) &
1644 * ( g12(i,1,1,e) * wur(i,1,1) &
1645 + g22(i,1,1,e) * wus(i,1,1) &
1646 + g23(i,1,1,e) * wut(i,1,1) )
1647 ut(i,1,1) = h1(i,1,1,e) &
1648 * ( g13(i,1,1,e) * wur(i,1,1) &
1649 + g23(i,1,1,e) * wus(i,1,1) &
1650 + g33(i,1,1,e) * wut(i,1,1) )
1651
1652 vr(i,1,1) = h1(i,1,1,e) &
1653 * ( g11(i,1,1,e) * wvr(i,1,1) &
1654 + g12(i,1,1,e) * wvs(i,1,1) &
1655 + g13(i,1,1,e) * wvt(i,1,1) )
1656 vs(i,1,1) = h1(i,1,1,e) &
1657 * ( g12(i,1,1,e) * wvr(i,1,1) &
1658 + g22(i,1,1,e) * wvs(i,1,1) &
1659 + g23(i,1,1,e) * wvt(i,1,1) )
1660 vt(i,1,1) = h1(i,1,1,e) &
1661 * ( g13(i,1,1,e) * wvr(i,1,1) &
1662 + g23(i,1,1,e) * wvs(i,1,1) &
1663 + g33(i,1,1,e) * wvt(i,1,1) )
1664
1665 wr(i,1,1) = h1(i,1,1,e) &
1666 * ( g11(i,1,1,e) * wwr(i,1,1) &
1667 + g12(i,1,1,e) * wws(i,1,1) &
1668 + g13(i,1,1,e) * wwt(i,1,1) )
1669 ws(i,1,1) = h1(i,1,1,e) &
1670 * ( g12(i,1,1,e) * wwr(i,1,1) &
1671 + g22(i,1,1,e) * wws(i,1,1) &
1672 + g23(i,1,1,e) * wwt(i,1,1) )
1673 wt(i,1,1) = h1(i,1,1,e) &
1674 * ( g13(i,1,1,e) * wwr(i,1,1) &
1675 + g23(i,1,1,e) * wws(i,1,1) &
1676 + g33(i,1,1,e) * wwt(i,1,1) )
1677 end do
1678
1679 do j = 1, lx*lx
1680 do i = 1, lx
1681 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
1682 + dxt(i,2) * ur(2,j,1) &
1683 + dxt(i,3) * ur(3,j,1) &
1684 + dxt(i,4) * ur(4,j,1) &
1685 + dxt(i,5) * ur(5,j,1) &
1686 + dxt(i,6) * ur(6,j,1) &
1687 + dxt(i,7) * ur(7,j,1) &
1688 + dxt(i,8) * ur(8,j,1) &
1689 + dxt(i,9) * ur(9,j,1) &
1690 + dxt(i,10) * ur(10,j,1) &
1691 + dxt(i,11) * ur(11,j,1)
1692
1693 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
1694 + dxt(i,2) * vr(2,j,1) &
1695 + dxt(i,3) * vr(3,j,1) &
1696 + dxt(i,4) * vr(4,j,1) &
1697 + dxt(i,5) * vr(5,j,1) &
1698 + dxt(i,6) * vr(6,j,1) &
1699 + dxt(i,7) * vr(7,j,1) &
1700 + dxt(i,8) * vr(8,j,1) &
1701 + dxt(i,9) * vr(9,j,1) &
1702 + dxt(i,10) * vr(10,j,1) &
1703 + dxt(i,11) * vr(11,j,1)
1704
1705 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
1706 + dxt(i,2) * wr(2,j,1) &
1707 + dxt(i,3) * wr(3,j,1) &
1708 + dxt(i,4) * wr(4,j,1) &
1709 + dxt(i,5) * wr(5,j,1) &
1710 + dxt(i,6) * wr(6,j,1) &
1711 + dxt(i,7) * wr(7,j,1) &
1712 + dxt(i,8) * wr(8,j,1) &
1713 + dxt(i,9) * wr(9,j,1) &
1714 + dxt(i,10) * wr(10,j,1) &
1715 + dxt(i,11) * wr(11,j,1)
1716 end do
1717 end do
1718
1719 do k = 1, lx
1720 do j = 1, lx
1721 do i = 1, lx
1722 au(i,j,k,e) = au(i,j,k,e) &
1723 + dyt(j,1) * us(i,1,k) &
1724 + dyt(j,2) * us(i,2,k) &
1725 + dyt(j,3) * us(i,3,k) &
1726 + dyt(j,4) * us(i,4,k) &
1727 + dyt(j,5) * us(i,5,k) &
1728 + dyt(j,6) * us(i,6,k) &
1729 + dyt(j,7) * us(i,7,k) &
1730 + dyt(j,8) * us(i,8,k) &
1731 + dyt(j,9) * us(i,9,k) &
1732 + dyt(j,10) * us(i,10,k) &
1733 + dyt(j,11) * us(i,11,k)
1734
1735 av(i,j,k,e) = av(i,j,k,e) &
1736 + dyt(j,1) * vs(i,1,k) &
1737 + dyt(j,2) * vs(i,2,k) &
1738 + dyt(j,3) * vs(i,3,k) &
1739 + dyt(j,4) * vs(i,4,k) &
1740 + dyt(j,5) * vs(i,5,k) &
1741 + dyt(j,6) * vs(i,6,k) &
1742 + dyt(j,7) * vs(i,7,k) &
1743 + dyt(j,8) * vs(i,8,k) &
1744 + dyt(j,9) * vs(i,9,k) &
1745 + dyt(j,10) * vs(i,10,k) &
1746 + dyt(j,11) * vs(i,11,k)
1747
1748 aw(i,j,k,e) = aw(i,j,k,e) &
1749 + dyt(j,1) * ws(i,1,k) &
1750 + dyt(j,2) * ws(i,2,k) &
1751 + dyt(j,3) * ws(i,3,k) &
1752 + dyt(j,4) * ws(i,4,k) &
1753 + dyt(j,5) * ws(i,5,k) &
1754 + dyt(j,6) * ws(i,6,k) &
1755 + dyt(j,7) * ws(i,7,k) &
1756 + dyt(j,8) * ws(i,8,k) &
1757 + dyt(j,9) * ws(i,9,k) &
1758 + dyt(j,10) * ws(i,10,k) &
1759 + dyt(j,11) * ws(i,11,k)
1760 end do
1761 end do
1762 end do
1763
1764 do k = 1, lx
1765 do i = 1, lx*lx
1766 au(i,1,k,e) = au(i,1,k,e) &
1767 + dzt(k,1) * ut(i,1,1) &
1768 + dzt(k,2) * ut(i,1,2) &
1769 + dzt(k,3) * ut(i,1,3) &
1770 + dzt(k,4) * ut(i,1,4) &
1771 + dzt(k,5) * ut(i,1,5) &
1772 + dzt(k,6) * ut(i,1,6) &
1773 + dzt(k,7) * ut(i,1,7) &
1774 + dzt(k,8) * ut(i,1,8) &
1775 + dzt(k,9) * ut(i,1,9) &
1776 + dzt(k,10) * ut(i,1,10) &
1777 + dzt(k,11) * ut(i,1,11)
1778
1779 av(i,1,k,e) = av(i,1,k,e) &
1780 + dzt(k,1) * vt(i,1,1) &
1781 + dzt(k,2) * vt(i,1,2) &
1782 + dzt(k,3) * vt(i,1,3) &
1783 + dzt(k,4) * vt(i,1,4) &
1784 + dzt(k,5) * vt(i,1,5) &
1785 + dzt(k,6) * vt(i,1,6) &
1786 + dzt(k,7) * vt(i,1,7) &
1787 + dzt(k,8) * vt(i,1,8) &
1788 + dzt(k,9) * vt(i,1,9) &
1789 + dzt(k,10) * vt(i,1,10) &
1790 + dzt(k,11) * vt(i,1,11)
1791
1792 aw(i,1,k,e) = aw(i,1,k,e) &
1793 + dzt(k,1) * wt(i,1,1) &
1794 + dzt(k,2) * wt(i,1,2) &
1795 + dzt(k,3) * wt(i,1,3) &
1796 + dzt(k,4) * wt(i,1,4) &
1797 + dzt(k,5) * wt(i,1,5) &
1798 + dzt(k,6) * wt(i,1,6) &
1799 + dzt(k,7) * wt(i,1,7) &
1800 + dzt(k,8) * wt(i,1,8) &
1801 + dzt(k,9) * wt(i,1,9) &
1802 + dzt(k,10) * wt(i,1,10) &
1803 + dzt(k,11) * wt(i,1,11)
1804 end do
1805 end do
1806
1807 end do
1808 !$omp end do
1809 end subroutine ax_helm_vector_lx11
1810
1811 subroutine ax_helm_vector_lx10(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
1812 h1, G11, G22, G33, G12, G13, G23, n)
1813 integer, parameter :: lx = 10
1814 integer, intent(in) :: n
1815 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
1816 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
1817 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
1818 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
1819 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
1820 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
1821 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
1822 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1823 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1824 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1825 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1826 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1827 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1828 real(kind=rp), intent(in) :: dx(lx, lx)
1829 real(kind=rp), intent(in) :: dy(lx, lx)
1830 real(kind=rp), intent(in) :: dz(lx, lx)
1831 real(kind=rp), intent(in) :: dxt(lx, lx)
1832 real(kind=rp), intent(in) :: dyt(lx, lx)
1833 real(kind=rp), intent(in) :: dzt(lx, lx)
1834 real(kind=rp) :: ur(lx, lx, lx)
1835 real(kind=rp) :: us(lx, lx, lx)
1836 real(kind=rp) :: ut(lx, lx, lx)
1837 real(kind=rp) :: vr(lx, lx, lx)
1838 real(kind=rp) :: vs(lx, lx, lx)
1839 real(kind=rp) :: vt(lx, lx, lx)
1840 real(kind=rp) :: wr(lx, lx, lx)
1841 real(kind=rp) :: ws(lx, lx, lx)
1842 real(kind=rp) :: wt(lx, lx, lx)
1843 real(kind=rp) :: wur(lx, lx, lx)
1844 real(kind=rp) :: wus(lx, lx, lx)
1845 real(kind=rp) :: wut(lx, lx, lx)
1846 real(kind=rp) :: wvr(lx, lx, lx)
1847 real(kind=rp) :: wvs(lx, lx, lx)
1848 real(kind=rp) :: wvt(lx, lx, lx)
1849 real(kind=rp) :: wwr(lx, lx, lx)
1850 real(kind=rp) :: wws(lx, lx, lx)
1851 real(kind=rp) :: wwt(lx, lx, lx)
1852 integer :: e, i, j, k
1853
1854 !$omp do
1855 do e = 1, n
1856 do j = 1, lx * lx
1857 do i = 1, lx
1858 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1859 + dx(i,2) * u(2,j,1,e) &
1860 + dx(i,3) * u(3,j,1,e) &
1861 + dx(i,4) * u(4,j,1,e) &
1862 + dx(i,5) * u(5,j,1,e) &
1863 + dx(i,6) * u(6,j,1,e) &
1864 + dx(i,7) * u(7,j,1,e) &
1865 + dx(i,8) * u(8,j,1,e) &
1866 + dx(i,9) * u(9,j,1,e) &
1867 + dx(i,10) * u(10,j,1,e)
1868
1869 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
1870 + dx(i,2) * v(2,j,1,e) &
1871 + dx(i,3) * v(3,j,1,e) &
1872 + dx(i,4) * v(4,j,1,e) &
1873 + dx(i,5) * v(5,j,1,e) &
1874 + dx(i,6) * v(6,j,1,e) &
1875 + dx(i,7) * v(7,j,1,e) &
1876 + dx(i,8) * v(8,j,1,e) &
1877 + dx(i,9) * v(9,j,1,e) &
1878 + dx(i,10) * v(10,j,1,e)
1879
1880 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
1881 + dx(i,2) * w(2,j,1,e) &
1882 + dx(i,3) * w(3,j,1,e) &
1883 + dx(i,4) * w(4,j,1,e) &
1884 + dx(i,5) * w(5,j,1,e) &
1885 + dx(i,6) * w(6,j,1,e) &
1886 + dx(i,7) * w(7,j,1,e) &
1887 + dx(i,8) * w(8,j,1,e) &
1888 + dx(i,9) * w(9,j,1,e) &
1889 + dx(i,10) * w(10,j,1,e)
1890 end do
1891 end do
1892
1893 do k = 1, lx
1894 do j = 1, lx
1895 do i = 1, lx
1896 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
1897 + dy(j,2) * u(i,2,k,e) &
1898 + dy(j,3) * u(i,3,k,e) &
1899 + dy(j,4) * u(i,4,k,e) &
1900 + dy(j,5) * u(i,5,k,e) &
1901 + dy(j,6) * u(i,6,k,e) &
1902 + dy(j,7) * u(i,7,k,e) &
1903 + dy(j,8) * u(i,8,k,e) &
1904 + dy(j,9) * u(i,9,k,e) &
1905 + dy(j,10) * u(i,10,k,e)
1906
1907 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
1908 + dy(j,2) * v(i,2,k,e) &
1909 + dy(j,3) * v(i,3,k,e) &
1910 + dy(j,4) * v(i,4,k,e) &
1911 + dy(j,5) * v(i,5,k,e) &
1912 + dy(j,6) * v(i,6,k,e) &
1913 + dy(j,7) * v(i,7,k,e) &
1914 + dy(j,8) * v(i,8,k,e) &
1915 + dy(j,9) * v(i,9,k,e) &
1916 + dy(j,10) * v(i,10,k,e)
1917
1918 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
1919 + dy(j,2) * w(i,2,k,e) &
1920 + dy(j,3) * w(i,3,k,e) &
1921 + dy(j,4) * w(i,4,k,e) &
1922 + dy(j,5) * w(i,5,k,e) &
1923 + dy(j,6) * w(i,6,k,e) &
1924 + dy(j,7) * w(i,7,k,e) &
1925 + dy(j,8) * w(i,8,k,e) &
1926 + dy(j,9) * w(i,9,k,e) &
1927 + dy(j,10) * w(i,10,k,e)
1928 end do
1929 end do
1930 end do
1931
1932 do k = 1, lx
1933 do i = 1, lx*lx
1934 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1935 + dz(k,2) * u(i,1,2,e) &
1936 + dz(k,3) * u(i,1,3,e) &
1937 + dz(k,4) * u(i,1,4,e) &
1938 + dz(k,5) * u(i,1,5,e) &
1939 + dz(k,6) * u(i,1,6,e) &
1940 + dz(k,7) * u(i,1,7,e) &
1941 + dz(k,8) * u(i,1,8,e) &
1942 + dz(k,9) * u(i,1,9,e) &
1943 + dz(k,10) * u(i,1,10,e)
1944
1945 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
1946 + dz(k,2) * v(i,1,2,e) &
1947 + dz(k,3) * v(i,1,3,e) &
1948 + dz(k,4) * v(i,1,4,e) &
1949 + dz(k,5) * v(i,1,5,e) &
1950 + dz(k,6) * v(i,1,6,e) &
1951 + dz(k,7) * v(i,1,7,e) &
1952 + dz(k,8) * v(i,1,8,e) &
1953 + dz(k,9) * v(i,1,9,e) &
1954 + dz(k,10) * v(i,1,10,e)
1955
1956 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
1957 + dz(k,2) * w(i,1,2,e) &
1958 + dz(k,3) * w(i,1,3,e) &
1959 + dz(k,4) * w(i,1,4,e) &
1960 + dz(k,5) * w(i,1,5,e) &
1961 + dz(k,6) * w(i,1,6,e) &
1962 + dz(k,7) * w(i,1,7,e) &
1963 + dz(k,8) * w(i,1,8,e) &
1964 + dz(k,9) * w(i,1,9,e) &
1965 + dz(k,10) * w(i,1,10,e)
1966 end do
1967 end do
1968
1969 do i = 1, lx*lx*lx
1970 ur(i,1,1) = h1(i,1,1,e) &
1971 * ( g11(i,1,1,e) * wur(i,1,1) &
1972 + g12(i,1,1,e) * wus(i,1,1) &
1973 + g13(i,1,1,e) * wut(i,1,1) )
1974 us(i,1,1) = h1(i,1,1,e) &
1975 * ( g12(i,1,1,e) * wur(i,1,1) &
1976 + g22(i,1,1,e) * wus(i,1,1) &
1977 + g23(i,1,1,e) * wut(i,1,1) )
1978 ut(i,1,1) = h1(i,1,1,e) &
1979 * ( g13(i,1,1,e) * wur(i,1,1) &
1980 + g23(i,1,1,e) * wus(i,1,1) &
1981 + g33(i,1,1,e) * wut(i,1,1) )
1982
1983 vr(i,1,1) = h1(i,1,1,e) &
1984 * ( g11(i,1,1,e) * wvr(i,1,1) &
1985 + g12(i,1,1,e) * wvs(i,1,1) &
1986 + g13(i,1,1,e) * wvt(i,1,1) )
1987 vs(i,1,1) = h1(i,1,1,e) &
1988 * ( g12(i,1,1,e) * wvr(i,1,1) &
1989 + g22(i,1,1,e) * wvs(i,1,1) &
1990 + g23(i,1,1,e) * wvt(i,1,1) )
1991 vt(i,1,1) = h1(i,1,1,e) &
1992 * ( g13(i,1,1,e) * wvr(i,1,1) &
1993 + g23(i,1,1,e) * wvs(i,1,1) &
1994 + g33(i,1,1,e) * wvt(i,1,1) )
1995
1996 wr(i,1,1) = h1(i,1,1,e) &
1997 * ( g11(i,1,1,e) * wwr(i,1,1) &
1998 + g12(i,1,1,e) * wws(i,1,1) &
1999 + g13(i,1,1,e) * wwt(i,1,1) )
2000 ws(i,1,1) = h1(i,1,1,e) &
2001 * ( g12(i,1,1,e) * wwr(i,1,1) &
2002 + g22(i,1,1,e) * wws(i,1,1) &
2003 + g23(i,1,1,e) * wwt(i,1,1) )
2004 wt(i,1,1) = h1(i,1,1,e) &
2005 * ( g13(i,1,1,e) * wwr(i,1,1) &
2006 + g23(i,1,1,e) * wws(i,1,1) &
2007 + g33(i,1,1,e) * wwt(i,1,1) )
2008 end do
2009
2010 do j = 1, lx*lx
2011 do i = 1, lx
2012 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
2013 + dxt(i,2) * ur(2,j,1) &
2014 + dxt(i,3) * ur(3,j,1) &
2015 + dxt(i,4) * ur(4,j,1) &
2016 + dxt(i,5) * ur(5,j,1) &
2017 + dxt(i,6) * ur(6,j,1) &
2018 + dxt(i,7) * ur(7,j,1) &
2019 + dxt(i,8) * ur(8,j,1) &
2020 + dxt(i,9) * ur(9,j,1) &
2021 + dxt(i,10) * ur(10,j,1)
2022
2023 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
2024 + dxt(i,2) * vr(2,j,1) &
2025 + dxt(i,3) * vr(3,j,1) &
2026 + dxt(i,4) * vr(4,j,1) &
2027 + dxt(i,5) * vr(5,j,1) &
2028 + dxt(i,6) * vr(6,j,1) &
2029 + dxt(i,7) * vr(7,j,1) &
2030 + dxt(i,8) * vr(8,j,1) &
2031 + dxt(i,9) * vr(9,j,1) &
2032 + dxt(i,10) * vr(10,j,1)
2033
2034 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
2035 + dxt(i,2) * wr(2,j,1) &
2036 + dxt(i,3) * wr(3,j,1) &
2037 + dxt(i,4) * wr(4,j,1) &
2038 + dxt(i,5) * wr(5,j,1) &
2039 + dxt(i,6) * wr(6,j,1) &
2040 + dxt(i,7) * wr(7,j,1) &
2041 + dxt(i,8) * wr(8,j,1) &
2042 + dxt(i,9) * wr(9,j,1) &
2043 + dxt(i,10) * wr(10,j,1)
2044 end do
2045 end do
2046
2047 do k = 1, lx
2048 do j = 1, lx
2049 do i = 1, lx
2050 au(i,j,k,e) = au(i,j,k,e) &
2051 + dyt(j,1) * us(i,1,k) &
2052 + dyt(j,2) * us(i,2,k) &
2053 + dyt(j,3) * us(i,3,k) &
2054 + dyt(j,4) * us(i,4,k) &
2055 + dyt(j,5) * us(i,5,k) &
2056 + dyt(j,6) * us(i,6,k) &
2057 + dyt(j,7) * us(i,7,k) &
2058 + dyt(j,8) * us(i,8,k) &
2059 + dyt(j,9) * us(i,9,k) &
2060 + dyt(j,10) * us(i,10,k)
2061
2062 av(i,j,k,e) = av(i,j,k,e) &
2063 + dyt(j,1) * vs(i,1,k) &
2064 + dyt(j,2) * vs(i,2,k) &
2065 + dyt(j,3) * vs(i,3,k) &
2066 + dyt(j,4) * vs(i,4,k) &
2067 + dyt(j,5) * vs(i,5,k) &
2068 + dyt(j,6) * vs(i,6,k) &
2069 + dyt(j,7) * vs(i,7,k) &
2070 + dyt(j,8) * vs(i,8,k) &
2071 + dyt(j,9) * vs(i,9,k) &
2072 + dyt(j,10) * vs(i,10,k)
2073
2074 aw(i,j,k,e) = aw(i,j,k,e) &
2075 + dyt(j,1) * ws(i,1,k) &
2076 + dyt(j,2) * ws(i,2,k) &
2077 + dyt(j,3) * ws(i,3,k) &
2078 + dyt(j,4) * ws(i,4,k) &
2079 + dyt(j,5) * ws(i,5,k) &
2080 + dyt(j,6) * ws(i,6,k) &
2081 + dyt(j,7) * ws(i,7,k) &
2082 + dyt(j,8) * ws(i,8,k) &
2083 + dyt(j,9) * ws(i,9,k) &
2084 + dyt(j,10) * ws(i,10,k)
2085 end do
2086 end do
2087 end do
2088
2089 do k = 1, lx
2090 do i = 1, lx*lx
2091 au(i,1,k,e) = au(i,1,k,e) &
2092 + dzt(k,1) * ut(i,1,1) &
2093 + dzt(k,2) * ut(i,1,2) &
2094 + dzt(k,3) * ut(i,1,3) &
2095 + dzt(k,4) * ut(i,1,4) &
2096 + dzt(k,5) * ut(i,1,5) &
2097 + dzt(k,6) * ut(i,1,6) &
2098 + dzt(k,7) * ut(i,1,7) &
2099 + dzt(k,8) * ut(i,1,8) &
2100 + dzt(k,9) * ut(i,1,9) &
2101 + dzt(k,10) * ut(i,1,10)
2102
2103 av(i,1,k,e) = av(i,1,k,e) &
2104 + dzt(k,1) * vt(i,1,1) &
2105 + dzt(k,2) * vt(i,1,2) &
2106 + dzt(k,3) * vt(i,1,3) &
2107 + dzt(k,4) * vt(i,1,4) &
2108 + dzt(k,5) * vt(i,1,5) &
2109 + dzt(k,6) * vt(i,1,6) &
2110 + dzt(k,7) * vt(i,1,7) &
2111 + dzt(k,8) * vt(i,1,8) &
2112 + dzt(k,9) * vt(i,1,9) &
2113 + dzt(k,10) * vt(i,1,10)
2114
2115 aw(i,1,k,e) = aw(i,1,k,e) &
2116 + dzt(k,1) * wt(i,1,1) &
2117 + dzt(k,2) * wt(i,1,2) &
2118 + dzt(k,3) * wt(i,1,3) &
2119 + dzt(k,4) * wt(i,1,4) &
2120 + dzt(k,5) * wt(i,1,5) &
2121 + dzt(k,6) * wt(i,1,6) &
2122 + dzt(k,7) * wt(i,1,7) &
2123 + dzt(k,8) * wt(i,1,8) &
2124 + dzt(k,9) * wt(i,1,9) &
2125 + dzt(k,10) * wt(i,1,10)
2126 end do
2127 end do
2128
2129 end do
2130 !$omp end do
2131 end subroutine ax_helm_vector_lx10
2132
2133 subroutine ax_helm_vector_lx9(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
2134 h1, G11, G22, G33, G12, G13, G23, n)
2135 integer, parameter :: lx = 9
2136 integer, intent(in) :: n
2137 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2138 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2139 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2140 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2141 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
2142 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
2143 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
2144 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
2145 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
2146 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
2147 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
2148 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
2149 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
2150 real(kind=rp), intent(in) :: dx(lx, lx)
2151 real(kind=rp), intent(in) :: dy(lx, lx)
2152 real(kind=rp), intent(in) :: dz(lx, lx)
2153 real(kind=rp), intent(in) :: dxt(lx, lx)
2154 real(kind=rp), intent(in) :: dyt(lx, lx)
2155 real(kind=rp), intent(in) :: dzt(lx, lx)
2156 real(kind=rp) :: ur(lx, lx, lx)
2157 real(kind=rp) :: us(lx, lx, lx)
2158 real(kind=rp) :: ut(lx, lx, lx)
2159 real(kind=rp) :: vr(lx, lx, lx)
2160 real(kind=rp) :: vs(lx, lx, lx)
2161 real(kind=rp) :: vt(lx, lx, lx)
2162 real(kind=rp) :: wr(lx, lx, lx)
2163 real(kind=rp) :: ws(lx, lx, lx)
2164 real(kind=rp) :: wt(lx, lx, lx)
2165 real(kind=rp) :: wur(lx, lx, lx)
2166 real(kind=rp) :: wus(lx, lx, lx)
2167 real(kind=rp) :: wut(lx, lx, lx)
2168 real(kind=rp) :: wvr(lx, lx, lx)
2169 real(kind=rp) :: wvs(lx, lx, lx)
2170 real(kind=rp) :: wvt(lx, lx, lx)
2171 real(kind=rp) :: wwr(lx, lx, lx)
2172 real(kind=rp) :: wws(lx, lx, lx)
2173 real(kind=rp) :: wwt(lx, lx, lx)
2174 integer :: e, i, j, k
2175
2176 !$omp do
2177 do e = 1, n
2178 do j = 1, lx * lx
2179 do i = 1, lx
2180 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
2181 + dx(i,2) * u(2,j,1,e) &
2182 + dx(i,3) * u(3,j,1,e) &
2183 + dx(i,4) * u(4,j,1,e) &
2184 + dx(i,5) * u(5,j,1,e) &
2185 + dx(i,6) * u(6,j,1,e) &
2186 + dx(i,7) * u(7,j,1,e) &
2187 + dx(i,8) * u(8,j,1,e) &
2188 + dx(i,9) * u(9,j,1,e)
2189
2190 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
2191 + dx(i,2) * v(2,j,1,e) &
2192 + dx(i,3) * v(3,j,1,e) &
2193 + dx(i,4) * v(4,j,1,e) &
2194 + dx(i,5) * v(5,j,1,e) &
2195 + dx(i,6) * v(6,j,1,e) &
2196 + dx(i,7) * v(7,j,1,e) &
2197 + dx(i,8) * v(8,j,1,e) &
2198 + dx(i,9) * v(9,j,1,e)
2199
2200 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
2201 + dx(i,2) * w(2,j,1,e) &
2202 + dx(i,3) * w(3,j,1,e) &
2203 + dx(i,4) * w(4,j,1,e) &
2204 + dx(i,5) * w(5,j,1,e) &
2205 + dx(i,6) * w(6,j,1,e) &
2206 + dx(i,7) * w(7,j,1,e) &
2207 + dx(i,8) * w(8,j,1,e) &
2208 + dx(i,9) * w(9,j,1,e)
2209 end do
2210 end do
2211
2212 do k = 1, lx
2213 do j = 1, lx
2214 do i = 1, lx
2215 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
2216 + dy(j,2) * u(i,2,k,e) &
2217 + dy(j,3) * u(i,3,k,e) &
2218 + dy(j,4) * u(i,4,k,e) &
2219 + dy(j,5) * u(i,5,k,e) &
2220 + dy(j,6) * u(i,6,k,e) &
2221 + dy(j,7) * u(i,7,k,e) &
2222 + dy(j,8) * u(i,8,k,e) &
2223 + dy(j,9) * u(i,9,k,e)
2224
2225 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
2226 + dy(j,2) * v(i,2,k,e) &
2227 + dy(j,3) * v(i,3,k,e) &
2228 + dy(j,4) * v(i,4,k,e) &
2229 + dy(j,5) * v(i,5,k,e) &
2230 + dy(j,6) * v(i,6,k,e) &
2231 + dy(j,7) * v(i,7,k,e) &
2232 + dy(j,8) * v(i,8,k,e) &
2233 + dy(j,9) * v(i,9,k,e)
2234
2235 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
2236 + dy(j,2) * w(i,2,k,e) &
2237 + dy(j,3) * w(i,3,k,e) &
2238 + dy(j,4) * w(i,4,k,e) &
2239 + dy(j,5) * w(i,5,k,e) &
2240 + dy(j,6) * w(i,6,k,e) &
2241 + dy(j,7) * w(i,7,k,e) &
2242 + dy(j,8) * w(i,8,k,e) &
2243 + dy(j,9) * w(i,9,k,e)
2244 end do
2245 end do
2246 end do
2247
2248 do k = 1, lx
2249 do i = 1, lx*lx
2250 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2251 + dz(k,2) * u(i,1,2,e) &
2252 + dz(k,3) * u(i,1,3,e) &
2253 + dz(k,4) * u(i,1,4,e) &
2254 + dz(k,5) * u(i,1,5,e) &
2255 + dz(k,6) * u(i,1,6,e) &
2256 + dz(k,7) * u(i,1,7,e) &
2257 + dz(k,8) * u(i,1,8,e) &
2258 + dz(k,9) * u(i,1,9,e)
2259
2260 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2261 + dz(k,2) * v(i,1,2,e) &
2262 + dz(k,3) * v(i,1,3,e) &
2263 + dz(k,4) * v(i,1,4,e) &
2264 + dz(k,5) * v(i,1,5,e) &
2265 + dz(k,6) * v(i,1,6,e) &
2266 + dz(k,7) * v(i,1,7,e) &
2267 + dz(k,8) * v(i,1,8,e) &
2268 + dz(k,9) * v(i,1,9,e)
2269
2270 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2271 + dz(k,2) * w(i,1,2,e) &
2272 + dz(k,3) * w(i,1,3,e) &
2273 + dz(k,4) * w(i,1,4,e) &
2274 + dz(k,5) * w(i,1,5,e) &
2275 + dz(k,6) * w(i,1,6,e) &
2276 + dz(k,7) * w(i,1,7,e) &
2277 + dz(k,8) * w(i,1,8,e) &
2278 + dz(k,9) * w(i,1,9,e)
2279 end do
2280 end do
2281
2282 do i = 1, lx*lx*lx
2283 ur(i,1,1) = h1(i,1,1,e) &
2284 * ( g11(i,1,1,e) * wur(i,1,1) &
2285 + g12(i,1,1,e) * wus(i,1,1) &
2286 + g13(i,1,1,e) * wut(i,1,1) )
2287 us(i,1,1) = h1(i,1,1,e) &
2288 * ( g12(i,1,1,e) * wur(i,1,1) &
2289 + g22(i,1,1,e) * wus(i,1,1) &
2290 + g23(i,1,1,e) * wut(i,1,1) )
2291 ut(i,1,1) = h1(i,1,1,e) &
2292 * ( g13(i,1,1,e) * wur(i,1,1) &
2293 + g23(i,1,1,e) * wus(i,1,1) &
2294 + g33(i,1,1,e) * wut(i,1,1) )
2295
2296 vr(i,1,1) = h1(i,1,1,e) &
2297 * ( g11(i,1,1,e) * wvr(i,1,1) &
2298 + g12(i,1,1,e) * wvs(i,1,1) &
2299 + g13(i,1,1,e) * wvt(i,1,1) )
2300 vs(i,1,1) = h1(i,1,1,e) &
2301 * ( g12(i,1,1,e) * wvr(i,1,1) &
2302 + g22(i,1,1,e) * wvs(i,1,1) &
2303 + g23(i,1,1,e) * wvt(i,1,1) )
2304 vt(i,1,1) = h1(i,1,1,e) &
2305 * ( g13(i,1,1,e) * wvr(i,1,1) &
2306 + g23(i,1,1,e) * wvs(i,1,1) &
2307 + g33(i,1,1,e) * wvt(i,1,1) )
2308
2309 wr(i,1,1) = h1(i,1,1,e) &
2310 * ( g11(i,1,1,e) * wwr(i,1,1) &
2311 + g12(i,1,1,e) * wws(i,1,1) &
2312 + g13(i,1,1,e) * wwt(i,1,1) )
2313 ws(i,1,1) = h1(i,1,1,e) &
2314 * ( g12(i,1,1,e) * wwr(i,1,1) &
2315 + g22(i,1,1,e) * wws(i,1,1) &
2316 + g23(i,1,1,e) * wwt(i,1,1) )
2317 wt(i,1,1) = h1(i,1,1,e) &
2318 * ( g13(i,1,1,e) * wwr(i,1,1) &
2319 + g23(i,1,1,e) * wws(i,1,1) &
2320 + g33(i,1,1,e) * wwt(i,1,1) )
2321 end do
2322
2323 do j = 1, lx*lx
2324 do i = 1, lx
2325 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
2326 + dxt(i,2) * ur(2,j,1) &
2327 + dxt(i,3) * ur(3,j,1) &
2328 + dxt(i,4) * ur(4,j,1) &
2329 + dxt(i,5) * ur(5,j,1) &
2330 + dxt(i,6) * ur(6,j,1) &
2331 + dxt(i,7) * ur(7,j,1) &
2332 + dxt(i,8) * ur(8,j,1) &
2333 + dxt(i,9) * ur(9,j,1)
2334
2335 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
2336 + dxt(i,2) * vr(2,j,1) &
2337 + dxt(i,3) * vr(3,j,1) &
2338 + dxt(i,4) * vr(4,j,1) &
2339 + dxt(i,5) * vr(5,j,1) &
2340 + dxt(i,6) * vr(6,j,1) &
2341 + dxt(i,7) * vr(7,j,1) &
2342 + dxt(i,8) * vr(8,j,1) &
2343 + dxt(i,9) * vr(9,j,1)
2344
2345 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
2346 + dxt(i,2) * wr(2,j,1) &
2347 + dxt(i,3) * wr(3,j,1) &
2348 + dxt(i,4) * wr(4,j,1) &
2349 + dxt(i,5) * wr(5,j,1) &
2350 + dxt(i,6) * wr(6,j,1) &
2351 + dxt(i,7) * wr(7,j,1) &
2352 + dxt(i,8) * wr(8,j,1) &
2353 + dxt(i,9) * wr(9,j,1)
2354 end do
2355 end do
2356
2357 do k = 1, lx
2358 do j = 1, lx
2359 do i = 1, lx
2360 au(i,j,k,e) = au(i,j,k,e) &
2361 + dyt(j,1) * us(i,1,k) &
2362 + dyt(j,2) * us(i,2,k) &
2363 + dyt(j,3) * us(i,3,k) &
2364 + dyt(j,4) * us(i,4,k) &
2365 + dyt(j,5) * us(i,5,k) &
2366 + dyt(j,6) * us(i,6,k) &
2367 + dyt(j,7) * us(i,7,k) &
2368 + dyt(j,8) * us(i,8,k) &
2369 + dyt(j,9) * us(i,9,k)
2370
2371 av(i,j,k,e) = av(i,j,k,e) &
2372 + dyt(j,1) * vs(i,1,k) &
2373 + dyt(j,2) * vs(i,2,k) &
2374 + dyt(j,3) * vs(i,3,k) &
2375 + dyt(j,4) * vs(i,4,k) &
2376 + dyt(j,5) * vs(i,5,k) &
2377 + dyt(j,6) * vs(i,6,k) &
2378 + dyt(j,7) * vs(i,7,k) &
2379 + dyt(j,8) * vs(i,8,k) &
2380 + dyt(j,9) * vs(i,9,k)
2381
2382 aw(i,j,k,e) = aw(i,j,k,e) &
2383 + dyt(j,1) * ws(i,1,k) &
2384 + dyt(j,2) * ws(i,2,k) &
2385 + dyt(j,3) * ws(i,3,k) &
2386 + dyt(j,4) * ws(i,4,k) &
2387 + dyt(j,5) * ws(i,5,k) &
2388 + dyt(j,6) * ws(i,6,k) &
2389 + dyt(j,7) * ws(i,7,k) &
2390 + dyt(j,8) * ws(i,8,k) &
2391 + dyt(j,9) * ws(i,9,k)
2392 end do
2393 end do
2394 end do
2395
2396 do k = 1, lx
2397 do i = 1, lx*lx
2398 au(i,1,k,e) = au(i,1,k,e) &
2399 + dzt(k,1) * ut(i,1,1) &
2400 + dzt(k,2) * ut(i,1,2) &
2401 + dzt(k,3) * ut(i,1,3) &
2402 + dzt(k,4) * ut(i,1,4) &
2403 + dzt(k,5) * ut(i,1,5) &
2404 + dzt(k,6) * ut(i,1,6) &
2405 + dzt(k,7) * ut(i,1,7) &
2406 + dzt(k,8) * ut(i,1,8) &
2407 + dzt(k,9) * ut(i,1,9)
2408
2409 av(i,1,k,e) = av(i,1,k,e) &
2410 + dzt(k,1) * vt(i,1,1) &
2411 + dzt(k,2) * vt(i,1,2) &
2412 + dzt(k,3) * vt(i,1,3) &
2413 + dzt(k,4) * vt(i,1,4) &
2414 + dzt(k,5) * vt(i,1,5) &
2415 + dzt(k,6) * vt(i,1,6) &
2416 + dzt(k,7) * vt(i,1,7) &
2417 + dzt(k,8) * vt(i,1,8) &
2418 + dzt(k,9) * vt(i,1,9)
2419
2420 aw(i,1,k,e) = aw(i,1,k,e) &
2421 + dzt(k,1) * wt(i,1,1) &
2422 + dzt(k,2) * wt(i,1,2) &
2423 + dzt(k,3) * wt(i,1,3) &
2424 + dzt(k,4) * wt(i,1,4) &
2425 + dzt(k,5) * wt(i,1,5) &
2426 + dzt(k,6) * wt(i,1,6) &
2427 + dzt(k,7) * wt(i,1,7) &
2428 + dzt(k,8) * wt(i,1,8) &
2429 + dzt(k,9) * wt(i,1,9)
2430 end do
2431 end do
2432
2433 end do
2434 !$omp end do
2435 end subroutine ax_helm_vector_lx9
2436
2437 subroutine ax_helm_vector_lx8(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
2438 h1, G11, G22, G33, G12, G13, G23, n)
2439 integer, parameter :: lx = 8
2440 integer, intent(in) :: n
2441 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2442 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2443 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2444 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2445 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
2446 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
2447 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
2448 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
2449 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
2450 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
2451 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
2452 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
2453 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
2454 real(kind=rp), intent(in) :: dx(lx, lx)
2455 real(kind=rp), intent(in) :: dy(lx, lx)
2456 real(kind=rp), intent(in) :: dz(lx, lx)
2457 real(kind=rp), intent(in) :: dxt(lx, lx)
2458 real(kind=rp), intent(in) :: dyt(lx, lx)
2459 real(kind=rp), intent(in) :: dzt(lx, lx)
2460 real(kind=rp) :: ur(lx, lx, lx)
2461 real(kind=rp) :: us(lx, lx, lx)
2462 real(kind=rp) :: ut(lx, lx, lx)
2463 real(kind=rp) :: vr(lx, lx, lx)
2464 real(kind=rp) :: vs(lx, lx, lx)
2465 real(kind=rp) :: vt(lx, lx, lx)
2466 real(kind=rp) :: wr(lx, lx, lx)
2467 real(kind=rp) :: ws(lx, lx, lx)
2468 real(kind=rp) :: wt(lx, lx, lx)
2469 real(kind=rp) :: wur(lx, lx, lx)
2470 real(kind=rp) :: wus(lx, lx, lx)
2471 real(kind=rp) :: wut(lx, lx, lx)
2472 real(kind=rp) :: wvr(lx, lx, lx)
2473 real(kind=rp) :: wvs(lx, lx, lx)
2474 real(kind=rp) :: wvt(lx, lx, lx)
2475 real(kind=rp) :: wwr(lx, lx, lx)
2476 real(kind=rp) :: wws(lx, lx, lx)
2477 real(kind=rp) :: wwt(lx, lx, lx)
2478 integer :: e, i, j, k
2479
2480 !$omp do
2481 do e = 1, n
2482 do j = 1, lx * lx
2483 do i = 1, lx
2484 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
2485 + dx(i,2) * u(2,j,1,e) &
2486 + dx(i,3) * u(3,j,1,e) &
2487 + dx(i,4) * u(4,j,1,e) &
2488 + dx(i,5) * u(5,j,1,e) &
2489 + dx(i,6) * u(6,j,1,e) &
2490 + dx(i,7) * u(7,j,1,e) &
2491 + dx(i,8) * u(8,j,1,e)
2492
2493 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
2494 + dx(i,2) * v(2,j,1,e) &
2495 + dx(i,3) * v(3,j,1,e) &
2496 + dx(i,4) * v(4,j,1,e) &
2497 + dx(i,5) * v(5,j,1,e) &
2498 + dx(i,6) * v(6,j,1,e) &
2499 + dx(i,7) * v(7,j,1,e) &
2500 + dx(i,8) * v(8,j,1,e)
2501
2502 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
2503 + dx(i,2) * w(2,j,1,e) &
2504 + dx(i,3) * w(3,j,1,e) &
2505 + dx(i,4) * w(4,j,1,e) &
2506 + dx(i,5) * w(5,j,1,e) &
2507 + dx(i,6) * w(6,j,1,e) &
2508 + dx(i,7) * w(7,j,1,e) &
2509 + dx(i,8) * w(8,j,1,e)
2510 end do
2511 end do
2512
2513 do k = 1, lx
2514 do j = 1, lx
2515 do i = 1, lx
2516 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
2517 + dy(j,2) * u(i,2,k,e) &
2518 + dy(j,3) * u(i,3,k,e) &
2519 + dy(j,4) * u(i,4,k,e) &
2520 + dy(j,5) * u(i,5,k,e) &
2521 + dy(j,6) * u(i,6,k,e) &
2522 + dy(j,7) * u(i,7,k,e) &
2523 + dy(j,8) * u(i,8,k,e)
2524
2525 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
2526 + dy(j,2) * v(i,2,k,e) &
2527 + dy(j,3) * v(i,3,k,e) &
2528 + dy(j,4) * v(i,4,k,e) &
2529 + dy(j,5) * v(i,5,k,e) &
2530 + dy(j,6) * v(i,6,k,e) &
2531 + dy(j,7) * v(i,7,k,e) &
2532 + dy(j,8) * v(i,8,k,e)
2533
2534 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
2535 + dy(j,2) * w(i,2,k,e) &
2536 + dy(j,3) * w(i,3,k,e) &
2537 + dy(j,4) * w(i,4,k,e) &
2538 + dy(j,5) * w(i,5,k,e) &
2539 + dy(j,6) * w(i,6,k,e) &
2540 + dy(j,7) * w(i,7,k,e) &
2541 + dy(j,8) * w(i,8,k,e)
2542 end do
2543 end do
2544 end do
2545
2546 do k = 1, lx
2547 do i = 1, lx*lx
2548 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2549 + dz(k,2) * u(i,1,2,e) &
2550 + dz(k,3) * u(i,1,3,e) &
2551 + dz(k,4) * u(i,1,4,e) &
2552 + dz(k,5) * u(i,1,5,e) &
2553 + dz(k,6) * u(i,1,6,e) &
2554 + dz(k,7) * u(i,1,7,e) &
2555 + dz(k,8) * u(i,1,8,e)
2556
2557 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2558 + dz(k,2) * v(i,1,2,e) &
2559 + dz(k,3) * v(i,1,3,e) &
2560 + dz(k,4) * v(i,1,4,e) &
2561 + dz(k,5) * v(i,1,5,e) &
2562 + dz(k,6) * v(i,1,6,e) &
2563 + dz(k,7) * v(i,1,7,e) &
2564 + dz(k,8) * v(i,1,8,e)
2565
2566 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2567 + dz(k,2) * w(i,1,2,e) &
2568 + dz(k,3) * w(i,1,3,e) &
2569 + dz(k,4) * w(i,1,4,e) &
2570 + dz(k,5) * w(i,1,5,e) &
2571 + dz(k,6) * w(i,1,6,e) &
2572 + dz(k,7) * w(i,1,7,e) &
2573 + dz(k,8) * w(i,1,8,e)
2574 end do
2575 end do
2576
2577 do i = 1, lx*lx*lx
2578 ur(i,1,1) = h1(i,1,1,e) &
2579 * ( g11(i,1,1,e) * wur(i,1,1) &
2580 + g12(i,1,1,e) * wus(i,1,1) &
2581 + g13(i,1,1,e) * wut(i,1,1) )
2582 us(i,1,1) = h1(i,1,1,e) &
2583 * ( g12(i,1,1,e) * wur(i,1,1) &
2584 + g22(i,1,1,e) * wus(i,1,1) &
2585 + g23(i,1,1,e) * wut(i,1,1) )
2586 ut(i,1,1) = h1(i,1,1,e) &
2587 * ( g13(i,1,1,e) * wur(i,1,1) &
2588 + g23(i,1,1,e) * wus(i,1,1) &
2589 + g33(i,1,1,e) * wut(i,1,1) )
2590
2591 vr(i,1,1) = h1(i,1,1,e) &
2592 * ( g11(i,1,1,e) * wvr(i,1,1) &
2593 + g12(i,1,1,e) * wvs(i,1,1) &
2594 + g13(i,1,1,e) * wvt(i,1,1) )
2595 vs(i,1,1) = h1(i,1,1,e) &
2596 * ( g12(i,1,1,e) * wvr(i,1,1) &
2597 + g22(i,1,1,e) * wvs(i,1,1) &
2598 + g23(i,1,1,e) * wvt(i,1,1) )
2599 vt(i,1,1) = h1(i,1,1,e) &
2600 * ( g13(i,1,1,e) * wvr(i,1,1) &
2601 + g23(i,1,1,e) * wvs(i,1,1) &
2602 + g33(i,1,1,e) * wvt(i,1,1) )
2603
2604 wr(i,1,1) = h1(i,1,1,e) &
2605 * ( g11(i,1,1,e) * wwr(i,1,1) &
2606 + g12(i,1,1,e) * wws(i,1,1) &
2607 + g13(i,1,1,e) * wwt(i,1,1) )
2608 ws(i,1,1) = h1(i,1,1,e) &
2609 * ( g12(i,1,1,e) * wwr(i,1,1) &
2610 + g22(i,1,1,e) * wws(i,1,1) &
2611 + g23(i,1,1,e) * wwt(i,1,1) )
2612 wt(i,1,1) = h1(i,1,1,e) &
2613 * ( g13(i,1,1,e) * wwr(i,1,1) &
2614 + g23(i,1,1,e) * wws(i,1,1) &
2615 + g33(i,1,1,e) * wwt(i,1,1) )
2616 end do
2617
2618 do j = 1, lx*lx
2619 do i = 1, lx
2620 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
2621 + dxt(i,2) * ur(2,j,1) &
2622 + dxt(i,3) * ur(3,j,1) &
2623 + dxt(i,4) * ur(4,j,1) &
2624 + dxt(i,5) * ur(5,j,1) &
2625 + dxt(i,6) * ur(6,j,1) &
2626 + dxt(i,7) * ur(7,j,1) &
2627 + dxt(i,8) * ur(8,j,1)
2628
2629 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
2630 + dxt(i,2) * vr(2,j,1) &
2631 + dxt(i,3) * vr(3,j,1) &
2632 + dxt(i,4) * vr(4,j,1) &
2633 + dxt(i,5) * vr(5,j,1) &
2634 + dxt(i,6) * vr(6,j,1) &
2635 + dxt(i,7) * vr(7,j,1) &
2636 + dxt(i,8) * vr(8,j,1)
2637
2638 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
2639 + dxt(i,2) * wr(2,j,1) &
2640 + dxt(i,3) * wr(3,j,1) &
2641 + dxt(i,4) * wr(4,j,1) &
2642 + dxt(i,5) * wr(5,j,1) &
2643 + dxt(i,6) * wr(6,j,1) &
2644 + dxt(i,7) * wr(7,j,1) &
2645 + dxt(i,8) * wr(8,j,1)
2646 end do
2647 end do
2648
2649 do k = 1, lx
2650 do j = 1, lx
2651 do i = 1, lx
2652 au(i,j,k,e) = au(i,j,k,e) &
2653 + dyt(j,1) * us(i,1,k) &
2654 + dyt(j,2) * us(i,2,k) &
2655 + dyt(j,3) * us(i,3,k) &
2656 + dyt(j,4) * us(i,4,k) &
2657 + dyt(j,5) * us(i,5,k) &
2658 + dyt(j,6) * us(i,6,k) &
2659 + dyt(j,7) * us(i,7,k) &
2660 + dyt(j,8) * us(i,8,k)
2661
2662 av(i,j,k,e) = av(i,j,k,e) &
2663 + dyt(j,1) * vs(i,1,k) &
2664 + dyt(j,2) * vs(i,2,k) &
2665 + dyt(j,3) * vs(i,3,k) &
2666 + dyt(j,4) * vs(i,4,k) &
2667 + dyt(j,5) * vs(i,5,k) &
2668 + dyt(j,6) * vs(i,6,k) &
2669 + dyt(j,7) * vs(i,7,k) &
2670 + dyt(j,8) * vs(i,8,k)
2671
2672 aw(i,j,k,e) = aw(i,j,k,e) &
2673 + dyt(j,1) * ws(i,1,k) &
2674 + dyt(j,2) * ws(i,2,k) &
2675 + dyt(j,3) * ws(i,3,k) &
2676 + dyt(j,4) * ws(i,4,k) &
2677 + dyt(j,5) * ws(i,5,k) &
2678 + dyt(j,6) * ws(i,6,k) &
2679 + dyt(j,7) * ws(i,7,k) &
2680 + dyt(j,8) * ws(i,8,k)
2681 end do
2682 end do
2683 end do
2684
2685 do k = 1, lx
2686 do i = 1, lx*lx
2687 au(i,1,k,e) = au(i,1,k,e) &
2688 + dzt(k,1) * ut(i,1,1) &
2689 + dzt(k,2) * ut(i,1,2) &
2690 + dzt(k,3) * ut(i,1,3) &
2691 + dzt(k,4) * ut(i,1,4) &
2692 + dzt(k,5) * ut(i,1,5) &
2693 + dzt(k,6) * ut(i,1,6) &
2694 + dzt(k,7) * ut(i,1,7) &
2695 + dzt(k,8) * ut(i,1,8)
2696
2697 av(i,1,k,e) = av(i,1,k,e) &
2698 + dzt(k,1) * vt(i,1,1) &
2699 + dzt(k,2) * vt(i,1,2) &
2700 + dzt(k,3) * vt(i,1,3) &
2701 + dzt(k,4) * vt(i,1,4) &
2702 + dzt(k,5) * vt(i,1,5) &
2703 + dzt(k,6) * vt(i,1,6) &
2704 + dzt(k,7) * vt(i,1,7) &
2705 + dzt(k,8) * vt(i,1,8)
2706
2707 au(i,1,k,e) = au(i,1,k,e) &
2708 + dzt(k,1) * wt(i,1,1) &
2709 + dzt(k,2) * wt(i,1,2) &
2710 + dzt(k,3) * wt(i,1,3) &
2711 + dzt(k,4) * wt(i,1,4) &
2712 + dzt(k,5) * wt(i,1,5) &
2713 + dzt(k,6) * wt(i,1,6) &
2714 + dzt(k,7) * wt(i,1,7) &
2715 + dzt(k,8) * wt(i,1,8)
2716 end do
2717 end do
2718
2719 end do
2720 !$omp end do
2721 end subroutine ax_helm_vector_lx8
2722
2723 subroutine ax_helm_vector_lx7(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
2724 h1, G11, G22, G33, G12, G13, G23, n)
2725 integer, parameter :: lx = 7
2726 integer, intent(in) :: n
2727 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2728 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2729 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2730 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2731 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
2732 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
2733 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
2734 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
2735 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
2736 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
2737 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
2738 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
2739 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
2740 real(kind=rp), intent(in) :: dx(lx, lx)
2741 real(kind=rp), intent(in) :: dy(lx, lx)
2742 real(kind=rp), intent(in) :: dz(lx, lx)
2743 real(kind=rp), intent(in) :: dxt(lx, lx)
2744 real(kind=rp), intent(in) :: dyt(lx, lx)
2745 real(kind=rp), intent(in) :: dzt(lx, lx)
2746 real(kind=rp) :: ur(lx, lx, lx)
2747 real(kind=rp) :: us(lx, lx, lx)
2748 real(kind=rp) :: ut(lx, lx, lx)
2749 real(kind=rp) :: vr(lx, lx, lx)
2750 real(kind=rp) :: vs(lx, lx, lx)
2751 real(kind=rp) :: vt(lx, lx, lx)
2752 real(kind=rp) :: wr(lx, lx, lx)
2753 real(kind=rp) :: ws(lx, lx, lx)
2754 real(kind=rp) :: wt(lx, lx, lx)
2755 real(kind=rp) :: wur(lx, lx, lx)
2756 real(kind=rp) :: wus(lx, lx, lx)
2757 real(kind=rp) :: wut(lx, lx, lx)
2758 real(kind=rp) :: wvr(lx, lx, lx)
2759 real(kind=rp) :: wvs(lx, lx, lx)
2760 real(kind=rp) :: wvt(lx, lx, lx)
2761 real(kind=rp) :: wwr(lx, lx, lx)
2762 real(kind=rp) :: wws(lx, lx, lx)
2763 real(kind=rp) :: wwt(lx, lx, lx)
2764 integer :: e, i, j, k
2765
2766 !$omp do
2767 do e = 1, n
2768 do j = 1, lx * lx
2769 do i = 1, lx
2770 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
2771 + dx(i,2) * u(2,j,1,e) &
2772 + dx(i,3) * u(3,j,1,e) &
2773 + dx(i,4) * u(4,j,1,e) &
2774 + dx(i,5) * u(5,j,1,e) &
2775 + dx(i,6) * u(6,j,1,e) &
2776 + dx(i,7) * u(7,j,1,e)
2777
2778 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
2779 + dx(i,2) * v(2,j,1,e) &
2780 + dx(i,3) * v(3,j,1,e) &
2781 + dx(i,4) * v(4,j,1,e) &
2782 + dx(i,5) * v(5,j,1,e) &
2783 + dx(i,6) * v(6,j,1,e) &
2784 + dx(i,7) * v(7,j,1,e)
2785
2786 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
2787 + dx(i,2) * w(2,j,1,e) &
2788 + dx(i,3) * w(3,j,1,e) &
2789 + dx(i,4) * w(4,j,1,e) &
2790 + dx(i,5) * w(5,j,1,e) &
2791 + dx(i,6) * w(6,j,1,e) &
2792 + dx(i,7) * w(7,j,1,e)
2793 end do
2794 end do
2795
2796 do k = 1, lx
2797 do j = 1, lx
2798 do i = 1, lx
2799 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
2800 + dy(j,2) * u(i,2,k,e) &
2801 + dy(j,3) * u(i,3,k,e) &
2802 + dy(j,4) * u(i,4,k,e) &
2803 + dy(j,5) * u(i,5,k,e) &
2804 + dy(j,6) * u(i,6,k,e) &
2805 + dy(j,7) * u(i,7,k,e)
2806
2807 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
2808 + dy(j,2) * v(i,2,k,e) &
2809 + dy(j,3) * v(i,3,k,e) &
2810 + dy(j,4) * v(i,4,k,e) &
2811 + dy(j,5) * v(i,5,k,e) &
2812 + dy(j,6) * v(i,6,k,e) &
2813 + dy(j,7) * v(i,7,k,e)
2814
2815 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
2816 + dy(j,2) * w(i,2,k,e) &
2817 + dy(j,3) * w(i,3,k,e) &
2818 + dy(j,4) * w(i,4,k,e) &
2819 + dy(j,5) * w(i,5,k,e) &
2820 + dy(j,6) * w(i,6,k,e) &
2821 + dy(j,7) * w(i,7,k,e)
2822 end do
2823 end do
2824 end do
2825
2826 do k = 1, lx
2827 do i = 1, lx*lx
2828 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
2829 + dz(k,2) * u(i,1,2,e) &
2830 + dz(k,3) * u(i,1,3,e) &
2831 + dz(k,4) * u(i,1,4,e) &
2832 + dz(k,5) * u(i,1,5,e) &
2833 + dz(k,6) * u(i,1,6,e) &
2834 + dz(k,7) * u(i,1,7,e)
2835
2836 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
2837 + dz(k,2) * v(i,1,2,e) &
2838 + dz(k,3) * v(i,1,3,e) &
2839 + dz(k,4) * v(i,1,4,e) &
2840 + dz(k,5) * v(i,1,5,e) &
2841 + dz(k,6) * v(i,1,6,e) &
2842 + dz(k,7) * v(i,1,7,e)
2843
2844 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
2845 + dz(k,2) * w(i,1,2,e) &
2846 + dz(k,3) * w(i,1,3,e) &
2847 + dz(k,4) * w(i,1,4,e) &
2848 + dz(k,5) * w(i,1,5,e) &
2849 + dz(k,6) * w(i,1,6,e) &
2850 + dz(k,7) * w(i,1,7,e)
2851 end do
2852 end do
2853
2854 do i = 1, lx*lx*lx
2855 ur(i,1,1) = h1(i,1,1,e) &
2856 * ( g11(i,1,1,e) * wur(i,1,1) &
2857 + g12(i,1,1,e) * wus(i,1,1) &
2858 + g13(i,1,1,e) * wut(i,1,1) )
2859 us(i,1,1) = h1(i,1,1,e) &
2860 * ( g12(i,1,1,e) * wur(i,1,1) &
2861 + g22(i,1,1,e) * wus(i,1,1) &
2862 + g23(i,1,1,e) * wut(i,1,1) )
2863 ut(i,1,1) = h1(i,1,1,e) &
2864 * ( g13(i,1,1,e) * wur(i,1,1) &
2865 + g23(i,1,1,e) * wus(i,1,1) &
2866 + g33(i,1,1,e) * wut(i,1,1) )
2867
2868 vr(i,1,1) = h1(i,1,1,e) &
2869 * ( g11(i,1,1,e) * wvr(i,1,1) &
2870 + g12(i,1,1,e) * wvs(i,1,1) &
2871 + g13(i,1,1,e) * wvt(i,1,1) )
2872 vs(i,1,1) = h1(i,1,1,e) &
2873 * ( g12(i,1,1,e) * wvr(i,1,1) &
2874 + g22(i,1,1,e) * wvs(i,1,1) &
2875 + g23(i,1,1,e) * wvt(i,1,1) )
2876 vt(i,1,1) = h1(i,1,1,e) &
2877 * ( g13(i,1,1,e) * wvr(i,1,1) &
2878 + g23(i,1,1,e) * wvs(i,1,1) &
2879 + g33(i,1,1,e) * wvt(i,1,1) )
2880
2881 wr(i,1,1) = h1(i,1,1,e) &
2882 * ( g11(i,1,1,e) * wwr(i,1,1) &
2883 + g12(i,1,1,e) * wws(i,1,1) &
2884 + g13(i,1,1,e) * wwt(i,1,1) )
2885 ws(i,1,1) = h1(i,1,1,e) &
2886 * ( g12(i,1,1,e) * wwr(i,1,1) &
2887 + g22(i,1,1,e) * wws(i,1,1) &
2888 + g23(i,1,1,e) * wwt(i,1,1) )
2889 wt(i,1,1) = h1(i,1,1,e) &
2890 * ( g13(i,1,1,e) * wwr(i,1,1) &
2891 + g23(i,1,1,e) * wws(i,1,1) &
2892 + g33(i,1,1,e) * wwt(i,1,1) )
2893 end do
2894
2895 do j = 1, lx*lx
2896 do i = 1, lx
2897 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
2898 + dxt(i,2) * ur(2,j,1) &
2899 + dxt(i,3) * ur(3,j,1) &
2900 + dxt(i,4) * ur(4,j,1) &
2901 + dxt(i,5) * ur(5,j,1) &
2902 + dxt(i,6) * ur(6,j,1) &
2903 + dxt(i,7) * ur(7,j,1)
2904
2905 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
2906 + dxt(i,2) * vr(2,j,1) &
2907 + dxt(i,3) * vr(3,j,1) &
2908 + dxt(i,4) * vr(4,j,1) &
2909 + dxt(i,5) * vr(5,j,1) &
2910 + dxt(i,6) * vr(6,j,1) &
2911 + dxt(i,7) * vr(7,j,1)
2912
2913 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
2914 + dxt(i,2) * wr(2,j,1) &
2915 + dxt(i,3) * wr(3,j,1) &
2916 + dxt(i,4) * wr(4,j,1) &
2917 + dxt(i,5) * wr(5,j,1) &
2918 + dxt(i,6) * wr(6,j,1) &
2919 + dxt(i,7) * wr(7,j,1)
2920 end do
2921 end do
2922
2923 do k = 1, lx
2924 do j = 1, lx
2925 do i = 1, lx
2926 au(i,j,k,e) = au(i,j,k,e) &
2927 + dyt(j,1) * us(i,1,k) &
2928 + dyt(j,2) * us(i,2,k) &
2929 + dyt(j,3) * us(i,3,k) &
2930 + dyt(j,4) * us(i,4,k) &
2931 + dyt(j,5) * us(i,5,k) &
2932 + dyt(j,6) * us(i,6,k) &
2933 + dyt(j,7) * us(i,7,k)
2934
2935 av(i,j,k,e) = av(i,j,k,e) &
2936 + dyt(j,1) * vs(i,1,k) &
2937 + dyt(j,2) * vs(i,2,k) &
2938 + dyt(j,3) * vs(i,3,k) &
2939 + dyt(j,4) * vs(i,4,k) &
2940 + dyt(j,5) * vs(i,5,k) &
2941 + dyt(j,6) * vs(i,6,k) &
2942 + dyt(j,7) * vs(i,7,k)
2943
2944 aw(i,j,k,e) = aw(i,j,k,e) &
2945 + dyt(j,1) * ws(i,1,k) &
2946 + dyt(j,2) * ws(i,2,k) &
2947 + dyt(j,3) * ws(i,3,k) &
2948 + dyt(j,4) * ws(i,4,k) &
2949 + dyt(j,5) * ws(i,5,k) &
2950 + dyt(j,6) * ws(i,6,k) &
2951 + dyt(j,7) * ws(i,7,k)
2952 end do
2953 end do
2954 end do
2955
2956 do k = 1, lx
2957 do i = 1, lx*lx
2958 au(i,1,k,e) = au(i,1,k,e) &
2959 + dzt(k,1) * ut(i,1,1) &
2960 + dzt(k,2) * ut(i,1,2) &
2961 + dzt(k,3) * ut(i,1,3) &
2962 + dzt(k,4) * ut(i,1,4) &
2963 + dzt(k,5) * ut(i,1,5) &
2964 + dzt(k,6) * ut(i,1,6) &
2965 + dzt(k,7) * ut(i,1,7)
2966
2967 av(i,1,k,e) = av(i,1,k,e) &
2968 + dzt(k,1) * vt(i,1,1) &
2969 + dzt(k,2) * vt(i,1,2) &
2970 + dzt(k,3) * vt(i,1,3) &
2971 + dzt(k,4) * vt(i,1,4) &
2972 + dzt(k,5) * vt(i,1,5) &
2973 + dzt(k,6) * vt(i,1,6) &
2974 + dzt(k,7) * vt(i,1,7)
2975
2976 aw(i,1,k,e) = aw(i,1,k,e) &
2977 + dzt(k,1) * wt(i,1,1) &
2978 + dzt(k,2) * wt(i,1,2) &
2979 + dzt(k,3) * wt(i,1,3) &
2980 + dzt(k,4) * wt(i,1,4) &
2981 + dzt(k,5) * wt(i,1,5) &
2982 + dzt(k,6) * wt(i,1,6) &
2983 + dzt(k,7) * wt(i,1,7)
2984 end do
2985 end do
2986
2987 end do
2988 !$omp end do
2989 end subroutine ax_helm_vector_lx7
2990
2991 subroutine ax_helm_vector_lx6(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
2992 h1, G11, G22, G33, G12, G13, G23, n)
2993 integer, parameter :: lx = 6
2994 integer, intent(in) :: n
2995 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
2996 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
2997 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
2998 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
2999 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3000 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3001 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3002 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
3003 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
3004 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
3005 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
3006 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
3007 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
3008 real(kind=rp), intent(in) :: dx(lx, lx)
3009 real(kind=rp), intent(in) :: dy(lx, lx)
3010 real(kind=rp), intent(in) :: dz(lx, lx)
3011 real(kind=rp), intent(in) :: dxt(lx, lx)
3012 real(kind=rp), intent(in) :: dyt(lx, lx)
3013 real(kind=rp), intent(in) :: dzt(lx, lx)
3014 real(kind=rp) :: ur(lx, lx, lx)
3015 real(kind=rp) :: us(lx, lx, lx)
3016 real(kind=rp) :: ut(lx, lx, lx)
3017 real(kind=rp) :: vr(lx, lx, lx)
3018 real(kind=rp) :: vs(lx, lx, lx)
3019 real(kind=rp) :: vt(lx, lx, lx)
3020 real(kind=rp) :: wr(lx, lx, lx)
3021 real(kind=rp) :: ws(lx, lx, lx)
3022 real(kind=rp) :: wt(lx, lx, lx)
3023 real(kind=rp) :: wur(lx, lx, lx)
3024 real(kind=rp) :: wus(lx, lx, lx)
3025 real(kind=rp) :: wut(lx, lx, lx)
3026 real(kind=rp) :: wvr(lx, lx, lx)
3027 real(kind=rp) :: wvs(lx, lx, lx)
3028 real(kind=rp) :: wvt(lx, lx, lx)
3029 real(kind=rp) :: wwr(lx, lx, lx)
3030 real(kind=rp) :: wws(lx, lx, lx)
3031 real(kind=rp) :: wwt(lx, lx, lx)
3032 integer :: e, i, j, k
3033
3034 !$omp do
3035 do e = 1, n
3036 do j = 1, lx * lx
3037 do i = 1, lx
3038 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3039 + dx(i,2) * u(2,j,1,e) &
3040 + dx(i,3) * u(3,j,1,e) &
3041 + dx(i,4) * u(4,j,1,e) &
3042 + dx(i,5) * u(5,j,1,e) &
3043 + dx(i,6) * u(6,j,1,e)
3044
3045 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3046 + dx(i,2) * v(2,j,1,e) &
3047 + dx(i,3) * v(3,j,1,e) &
3048 + dx(i,4) * v(4,j,1,e) &
3049 + dx(i,5) * v(5,j,1,e) &
3050 + dx(i,6) * v(6,j,1,e)
3051
3052 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3053 + dx(i,2) * w(2,j,1,e) &
3054 + dx(i,3) * w(3,j,1,e) &
3055 + dx(i,4) * w(4,j,1,e) &
3056 + dx(i,5) * w(5,j,1,e) &
3057 + dx(i,6) * w(6,j,1,e)
3058 end do
3059 end do
3060
3061 do k = 1, lx
3062 do j = 1, lx
3063 do i = 1, lx
3064 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3065 + dy(j,2) * u(i,2,k,e) &
3066 + dy(j,3) * u(i,3,k,e) &
3067 + dy(j,4) * u(i,4,k,e) &
3068 + dy(j,5) * u(i,5,k,e) &
3069 + dy(j,6) * u(i,6,k,e)
3070
3071 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3072 + dy(j,2) * v(i,2,k,e) &
3073 + dy(j,3) * v(i,3,k,e) &
3074 + dy(j,4) * v(i,4,k,e) &
3075 + dy(j,5) * v(i,5,k,e) &
3076 + dy(j,6) * v(i,6,k,e)
3077
3078 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3079 + dy(j,2) * w(i,2,k,e) &
3080 + dy(j,3) * w(i,3,k,e) &
3081 + dy(j,4) * w(i,4,k,e) &
3082 + dy(j,5) * w(i,5,k,e) &
3083 + dy(j,6) * w(i,6,k,e)
3084 end do
3085 end do
3086 end do
3087
3088 do k = 1, lx
3089 do i = 1, lx*lx
3090 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3091 + dz(k,2) * u(i,1,2,e) &
3092 + dz(k,3) * u(i,1,3,e) &
3093 + dz(k,4) * u(i,1,4,e) &
3094 + dz(k,5) * u(i,1,5,e) &
3095 + dz(k,6) * u(i,1,6,e)
3096
3097 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3098 + dz(k,2) * v(i,1,2,e) &
3099 + dz(k,3) * v(i,1,3,e) &
3100 + dz(k,4) * v(i,1,4,e) &
3101 + dz(k,5) * v(i,1,5,e) &
3102 + dz(k,6) * v(i,1,6,e)
3103
3104 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3105 + dz(k,2) * w(i,1,2,e) &
3106 + dz(k,3) * w(i,1,3,e) &
3107 + dz(k,4) * w(i,1,4,e) &
3108 + dz(k,5) * w(i,1,5,e) &
3109 + dz(k,6) * w(i,1,6,e)
3110 end do
3111 end do
3112
3113 do i = 1, lx*lx*lx
3114 ur(i,1,1) = h1(i,1,1,e) &
3115 * ( g11(i,1,1,e) * wur(i,1,1) &
3116 + g12(i,1,1,e) * wus(i,1,1) &
3117 + g13(i,1,1,e) * wut(i,1,1) )
3118 us(i,1,1) = h1(i,1,1,e) &
3119 * ( g12(i,1,1,e) * wur(i,1,1) &
3120 + g22(i,1,1,e) * wus(i,1,1) &
3121 + g23(i,1,1,e) * wut(i,1,1) )
3122 ut(i,1,1) = h1(i,1,1,e) &
3123 * ( g13(i,1,1,e) * wur(i,1,1) &
3124 + g23(i,1,1,e) * wus(i,1,1) &
3125 + g33(i,1,1,e) * wut(i,1,1) )
3126
3127 vr(i,1,1) = h1(i,1,1,e) &
3128 * ( g11(i,1,1,e) * wvr(i,1,1) &
3129 + g12(i,1,1,e) * wvs(i,1,1) &
3130 + g13(i,1,1,e) * wvt(i,1,1) )
3131 vs(i,1,1) = h1(i,1,1,e) &
3132 * ( g12(i,1,1,e) * wvr(i,1,1) &
3133 + g22(i,1,1,e) * wvs(i,1,1) &
3134 + g23(i,1,1,e) * wvt(i,1,1) )
3135 vt(i,1,1) = h1(i,1,1,e) &
3136 * ( g13(i,1,1,e) * wvr(i,1,1) &
3137 + g23(i,1,1,e) * wvs(i,1,1) &
3138 + g33(i,1,1,e) * wvt(i,1,1) )
3139
3140 wr(i,1,1) = h1(i,1,1,e) &
3141 * ( g11(i,1,1,e) * wwr(i,1,1) &
3142 + g12(i,1,1,e) * wws(i,1,1) &
3143 + g13(i,1,1,e) * wwt(i,1,1) )
3144 ws(i,1,1) = h1(i,1,1,e) &
3145 * ( g12(i,1,1,e) * wwr(i,1,1) &
3146 + g22(i,1,1,e) * wws(i,1,1) &
3147 + g23(i,1,1,e) * wwt(i,1,1) )
3148 wt(i,1,1) = h1(i,1,1,e) &
3149 * ( g13(i,1,1,e) * wwr(i,1,1) &
3150 + g23(i,1,1,e) * wws(i,1,1) &
3151 + g33(i,1,1,e) * wwt(i,1,1) )
3152 end do
3153
3154 do j = 1, lx*lx
3155 do i = 1, lx
3156 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
3157 + dxt(i,2) * ur(2,j,1) &
3158 + dxt(i,3) * ur(3,j,1) &
3159 + dxt(i,4) * ur(4,j,1) &
3160 + dxt(i,5) * ur(5,j,1) &
3161 + dxt(i,6) * ur(6,j,1)
3162
3163 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
3164 + dxt(i,2) * vr(2,j,1) &
3165 + dxt(i,3) * vr(3,j,1) &
3166 + dxt(i,4) * vr(4,j,1) &
3167 + dxt(i,5) * vr(5,j,1) &
3168 + dxt(i,6) * vr(6,j,1)
3169
3170 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
3171 + dxt(i,2) * wr(2,j,1) &
3172 + dxt(i,3) * wr(3,j,1) &
3173 + dxt(i,4) * wr(4,j,1) &
3174 + dxt(i,5) * wr(5,j,1) &
3175 + dxt(i,6) * wr(6,j,1)
3176 end do
3177 end do
3178
3179 do k = 1, lx
3180 do j = 1, lx
3181 do i = 1, lx
3182 au(i,j,k,e) = au(i,j,k,e) &
3183 + dyt(j,1) * us(i,1,k) &
3184 + dyt(j,2) * us(i,2,k) &
3185 + dyt(j,3) * us(i,3,k) &
3186 + dyt(j,4) * us(i,4,k) &
3187 + dyt(j,5) * us(i,5,k) &
3188 + dyt(j,6) * us(i,6,k)
3189
3190 av(i,j,k,e) = av(i,j,k,e) &
3191 + dyt(j,1) * vs(i,1,k) &
3192 + dyt(j,2) * vs(i,2,k) &
3193 + dyt(j,3) * vs(i,3,k) &
3194 + dyt(j,4) * vs(i,4,k) &
3195 + dyt(j,5) * vs(i,5,k) &
3196 + dyt(j,6) * vs(i,6,k)
3197
3198 aw(i,j,k,e) = aw(i,j,k,e) &
3199 + dyt(j,1) * ws(i,1,k) &
3200 + dyt(j,2) * ws(i,2,k) &
3201 + dyt(j,3) * ws(i,3,k) &
3202 + dyt(j,4) * ws(i,4,k) &
3203 + dyt(j,5) * ws(i,5,k) &
3204 + dyt(j,6) * ws(i,6,k)
3205 end do
3206 end do
3207 end do
3208
3209 do k = 1, lx
3210 do i = 1, lx*lx
3211 au(i,1,k,e) = au(i,1,k,e) &
3212 + dzt(k,1) * ut(i,1,1) &
3213 + dzt(k,2) * ut(i,1,2) &
3214 + dzt(k,3) * ut(i,1,3) &
3215 + dzt(k,4) * ut(i,1,4) &
3216 + dzt(k,5) * ut(i,1,5) &
3217 + dzt(k,6) * ut(i,1,6)
3218
3219 av(i,1,k,e) = av(i,1,k,e) &
3220 + dzt(k,1) * vt(i,1,1) &
3221 + dzt(k,2) * vt(i,1,2) &
3222 + dzt(k,3) * vt(i,1,3) &
3223 + dzt(k,4) * vt(i,1,4) &
3224 + dzt(k,5) * vt(i,1,5) &
3225 + dzt(k,6) * vt(i,1,6)
3226
3227 aw(i,1,k,e) = aw(i,1,k,e) &
3228 + dzt(k,1) * wt(i,1,1) &
3229 + dzt(k,2) * wt(i,1,2) &
3230 + dzt(k,3) * wt(i,1,3) &
3231 + dzt(k,4) * wt(i,1,4) &
3232 + dzt(k,5) * wt(i,1,5) &
3233 + dzt(k,6) * wt(i,1,6)
3234 end do
3235 end do
3236
3237 end do
3238 !$omp end do
3239 end subroutine ax_helm_vector_lx6
3240
3241 subroutine ax_helm_vector_lx5(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
3242 h1, G11, G22, G33, G12, G13, G23, n)
3243 integer, parameter :: lx = 5
3244 integer, intent(in) :: n
3245 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3246 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3247 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3248 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3249 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3250 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3251 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3252 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
3253 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
3254 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
3255 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
3256 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
3257 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
3258 real(kind=rp), intent(in) :: dx(lx, lx)
3259 real(kind=rp), intent(in) :: dy(lx, lx)
3260 real(kind=rp), intent(in) :: dz(lx, lx)
3261 real(kind=rp), intent(in) :: dxt(lx, lx)
3262 real(kind=rp), intent(in) :: dyt(lx, lx)
3263 real(kind=rp), intent(in) :: dzt(lx, lx)
3264 real(kind=rp) :: ur(lx, lx, lx)
3265 real(kind=rp) :: us(lx, lx, lx)
3266 real(kind=rp) :: ut(lx, lx, lx)
3267 real(kind=rp) :: vr(lx, lx, lx)
3268 real(kind=rp) :: vs(lx, lx, lx)
3269 real(kind=rp) :: vt(lx, lx, lx)
3270 real(kind=rp) :: wr(lx, lx, lx)
3271 real(kind=rp) :: ws(lx, lx, lx)
3272 real(kind=rp) :: wt(lx, lx, lx)
3273 real(kind=rp) :: wur(lx, lx, lx)
3274 real(kind=rp) :: wus(lx, lx, lx)
3275 real(kind=rp) :: wut(lx, lx, lx)
3276 real(kind=rp) :: wvr(lx, lx, lx)
3277 real(kind=rp) :: wvs(lx, lx, lx)
3278 real(kind=rp) :: wvt(lx, lx, lx)
3279 real(kind=rp) :: wwr(lx, lx, lx)
3280 real(kind=rp) :: wws(lx, lx, lx)
3281 real(kind=rp) :: wwt(lx, lx, lx)
3282 integer :: e, i, j, k
3283
3284 !$omp do
3285 do e = 1, n
3286 do j = 1, lx * lx
3287 do i = 1, lx
3288 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3289 + dx(i,2) * u(2,j,1,e) &
3290 + dx(i,3) * u(3,j,1,e) &
3291 + dx(i,4) * u(4,j,1,e) &
3292 + dx(i,5) * u(5,j,1,e)
3293
3294 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3295 + dx(i,2) * v(2,j,1,e) &
3296 + dx(i,3) * v(3,j,1,e) &
3297 + dx(i,4) * v(4,j,1,e) &
3298 + dx(i,5) * v(5,j,1,e)
3299
3300 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3301 + dx(i,2) * w(2,j,1,e) &
3302 + dx(i,3) * w(3,j,1,e) &
3303 + dx(i,4) * w(4,j,1,e) &
3304 + dx(i,5) * w(5,j,1,e)
3305 end do
3306 end do
3307
3308 do k = 1, lx
3309 do j = 1, lx
3310 do i = 1, lx
3311 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3312 + dy(j,2) * u(i,2,k,e) &
3313 + dy(j,3) * u(i,3,k,e) &
3314 + dy(j,4) * u(i,4,k,e) &
3315 + dy(j,5) * u(i,5,k,e)
3316
3317 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3318 + dy(j,2) * v(i,2,k,e) &
3319 + dy(j,3) * v(i,3,k,e) &
3320 + dy(j,4) * v(i,4,k,e) &
3321 + dy(j,5) * v(i,5,k,e)
3322
3323 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3324 + dy(j,2) * w(i,2,k,e) &
3325 + dy(j,3) * w(i,3,k,e) &
3326 + dy(j,4) * w(i,4,k,e) &
3327 + dy(j,5) * w(i,5,k,e)
3328 end do
3329 end do
3330 end do
3331
3332 do k = 1, lx
3333 do i = 1, lx*lx
3334 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3335 + dz(k,2) * u(i,1,2,e) &
3336 + dz(k,3) * u(i,1,3,e) &
3337 + dz(k,4) * u(i,1,4,e) &
3338 + dz(k,5) * u(i,1,5,e)
3339
3340 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3341 + dz(k,2) * v(i,1,2,e) &
3342 + dz(k,3) * v(i,1,3,e) &
3343 + dz(k,4) * v(i,1,4,e) &
3344 + dz(k,5) * v(i,1,5,e)
3345
3346 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3347 + dz(k,2) * w(i,1,2,e) &
3348 + dz(k,3) * w(i,1,3,e) &
3349 + dz(k,4) * w(i,1,4,e) &
3350 + dz(k,5) * w(i,1,5,e)
3351 end do
3352 end do
3353
3354 do i = 1, lx*lx*lx
3355 ur(i,1,1) = h1(i,1,1,e) &
3356 * ( g11(i,1,1,e) * wur(i,1,1) &
3357 + g12(i,1,1,e) * wus(i,1,1) &
3358 + g13(i,1,1,e) * wut(i,1,1) )
3359 us(i,1,1) = h1(i,1,1,e) &
3360 * ( g12(i,1,1,e) * wur(i,1,1) &
3361 + g22(i,1,1,e) * wus(i,1,1) &
3362 + g23(i,1,1,e) * wut(i,1,1) )
3363 ut(i,1,1) = h1(i,1,1,e) &
3364 * ( g13(i,1,1,e) * wur(i,1,1) &
3365 + g23(i,1,1,e) * wus(i,1,1) &
3366 + g33(i,1,1,e) * wut(i,1,1) )
3367
3368 vr(i,1,1) = h1(i,1,1,e) &
3369 * ( g11(i,1,1,e) * wvr(i,1,1) &
3370 + g12(i,1,1,e) * wvs(i,1,1) &
3371 + g13(i,1,1,e) * wvt(i,1,1) )
3372 vs(i,1,1) = h1(i,1,1,e) &
3373 * ( g12(i,1,1,e) * wvr(i,1,1) &
3374 + g22(i,1,1,e) * wvs(i,1,1) &
3375 + g23(i,1,1,e) * wvt(i,1,1) )
3376 vt(i,1,1) = h1(i,1,1,e) &
3377 * ( g13(i,1,1,e) * wvr(i,1,1) &
3378 + g23(i,1,1,e) * wvs(i,1,1) &
3379 + g33(i,1,1,e) * wvt(i,1,1) )
3380
3381 wr(i,1,1) = h1(i,1,1,e) &
3382 * ( g11(i,1,1,e) * wwr(i,1,1) &
3383 + g12(i,1,1,e) * wws(i,1,1) &
3384 + g13(i,1,1,e) * wwt(i,1,1) )
3385 ws(i,1,1) = h1(i,1,1,e) &
3386 * ( g12(i,1,1,e) * wwr(i,1,1) &
3387 + g22(i,1,1,e) * wws(i,1,1) &
3388 + g23(i,1,1,e) * wwt(i,1,1) )
3389 wt(i,1,1) = h1(i,1,1,e) &
3390 * ( g13(i,1,1,e) * wwr(i,1,1) &
3391 + g23(i,1,1,e) * wws(i,1,1) &
3392 + g33(i,1,1,e) * wwt(i,1,1) )
3393 end do
3394
3395 do j = 1, lx*lx
3396 do i = 1, lx
3397 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
3398 + dxt(i,2) * ur(2,j,1) &
3399 + dxt(i,3) * ur(3,j,1) &
3400 + dxt(i,4) * ur(4,j,1) &
3401 + dxt(i,5) * ur(5,j,1)
3402
3403 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
3404 + dxt(i,2) * vr(2,j,1) &
3405 + dxt(i,3) * vr(3,j,1) &
3406 + dxt(i,4) * vr(4,j,1) &
3407 + dxt(i,5) * vr(5,j,1)
3408
3409 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
3410 + dxt(i,2) * wr(2,j,1) &
3411 + dxt(i,3) * wr(3,j,1) &
3412 + dxt(i,4) * wr(4,j,1) &
3413 + dxt(i,5) * wr(5,j,1)
3414 end do
3415 end do
3416
3417 do k = 1, lx
3418 do j = 1, lx
3419 do i = 1, lx
3420 au(i,j,k,e) = au(i,j,k,e) &
3421 + dyt(j,1) * us(i,1,k) &
3422 + dyt(j,2) * us(i,2,k) &
3423 + dyt(j,3) * us(i,3,k) &
3424 + dyt(j,4) * us(i,4,k) &
3425 + dyt(j,5) * us(i,5,k)
3426
3427 av(i,j,k,e) = av(i,j,k,e) &
3428 + dyt(j,1) * vs(i,1,k) &
3429 + dyt(j,2) * vs(i,2,k) &
3430 + dyt(j,3) * vs(i,3,k) &
3431 + dyt(j,4) * vs(i,4,k) &
3432 + dyt(j,5) * vs(i,5,k)
3433
3434 aw(i,j,k,e) = aw(i,j,k,e) &
3435 + dyt(j,1) * ws(i,1,k) &
3436 + dyt(j,2) * ws(i,2,k) &
3437 + dyt(j,3) * ws(i,3,k) &
3438 + dyt(j,4) * ws(i,4,k) &
3439 + dyt(j,5) * ws(i,5,k)
3440 end do
3441 end do
3442 end do
3443
3444 do k = 1, lx
3445 do i = 1, lx*lx
3446 au(i,1,k,e) = au(i,1,k,e) &
3447 + dzt(k,1) * ut(i,1,1) &
3448 + dzt(k,2) * ut(i,1,2) &
3449 + dzt(k,3) * ut(i,1,3) &
3450 + dzt(k,4) * ut(i,1,4) &
3451 + dzt(k,5) * ut(i,1,5)
3452
3453 av(i,1,k,e) = av(i,1,k,e) &
3454 + dzt(k,1) * vt(i,1,1) &
3455 + dzt(k,2) * vt(i,1,2) &
3456 + dzt(k,3) * vt(i,1,3) &
3457 + dzt(k,4) * vt(i,1,4) &
3458 + dzt(k,5) * vt(i,1,5)
3459
3460 aw(i,1,k,e) = aw(i,1,k,e) &
3461 + dzt(k,1) * wt(i,1,1) &
3462 + dzt(k,2) * wt(i,1,2) &
3463 + dzt(k,3) * wt(i,1,3) &
3464 + dzt(k,4) * wt(i,1,4) &
3465 + dzt(k,5) * wt(i,1,5)
3466 end do
3467 end do
3468
3469 end do
3470 !$omp end do
3471 end subroutine ax_helm_vector_lx5
3472
3473 subroutine ax_helm_vector_lx4(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
3474 h1, G11, G22, G33, G12, G13, G23, n)
3475 integer, parameter :: lx = 4
3476 integer, intent(in) :: n
3477 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3478 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3479 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3480 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3481 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3482 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3483 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3484 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
3485 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
3486 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
3487 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
3488 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
3489 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
3490 real(kind=rp), intent(in) :: dx(lx, lx)
3491 real(kind=rp), intent(in) :: dy(lx, lx)
3492 real(kind=rp), intent(in) :: dz(lx, lx)
3493 real(kind=rp), intent(in) :: dxt(lx, lx)
3494 real(kind=rp), intent(in) :: dyt(lx, lx)
3495 real(kind=rp), intent(in) :: dzt(lx, lx)
3496 real(kind=rp) :: ur(lx, lx, lx)
3497 real(kind=rp) :: us(lx, lx, lx)
3498 real(kind=rp) :: ut(lx, lx, lx)
3499 real(kind=rp) :: vr(lx, lx, lx)
3500 real(kind=rp) :: vs(lx, lx, lx)
3501 real(kind=rp) :: vt(lx, lx, lx)
3502 real(kind=rp) :: wr(lx, lx, lx)
3503 real(kind=rp) :: ws(lx, lx, lx)
3504 real(kind=rp) :: wt(lx, lx, lx)
3505 real(kind=rp) :: wur(lx, lx, lx)
3506 real(kind=rp) :: wus(lx, lx, lx)
3507 real(kind=rp) :: wut(lx, lx, lx)
3508 real(kind=rp) :: wvr(lx, lx, lx)
3509 real(kind=rp) :: wvs(lx, lx, lx)
3510 real(kind=rp) :: wvt(lx, lx, lx)
3511 real(kind=rp) :: wwr(lx, lx, lx)
3512 real(kind=rp) :: wws(lx, lx, lx)
3513 real(kind=rp) :: wwt(lx, lx, lx)
3514 integer :: e, i, j, k
3515
3516 !$omp do
3517 do e = 1, n
3518 do j = 1, lx * lx
3519 do i = 1, lx
3520 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3521 + dx(i,2) * u(2,j,1,e) &
3522 + dx(i,3) * u(3,j,1,e) &
3523 + dx(i,4) * u(4,j,1,e)
3524
3525 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3526 + dx(i,2) * v(2,j,1,e) &
3527 + dx(i,3) * v(3,j,1,e) &
3528 + dx(i,4) * v(4,j,1,e)
3529
3530 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3531 + dx(i,2) * w(2,j,1,e) &
3532 + dx(i,3) * w(3,j,1,e) &
3533 + dx(i,4) * w(4,j,1,e)
3534 end do
3535 end do
3536
3537 do k = 1, lx
3538 do j = 1, lx
3539 do i = 1, lx
3540 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3541 + dy(j,2) * u(i,2,k,e) &
3542 + dy(j,3) * u(i,3,k,e) &
3543 + dy(j,4) * u(i,4,k,e)
3544
3545 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3546 + dy(j,2) * v(i,2,k,e) &
3547 + dy(j,3) * v(i,3,k,e) &
3548 + dy(j,4) * v(i,4,k,e)
3549
3550 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3551 + dy(j,2) * w(i,2,k,e) &
3552 + dy(j,3) * w(i,3,k,e) &
3553 + dy(j,4) * w(i,4,k,e)
3554 end do
3555 end do
3556 end do
3557
3558 do k = 1, lx
3559 do i = 1, lx*lx
3560 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3561 + dz(k,2) * u(i,1,2,e) &
3562 + dz(k,3) * u(i,1,3,e) &
3563 + dz(k,4) * u(i,1,4,e)
3564
3565 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3566 + dz(k,2) * v(i,1,2,e) &
3567 + dz(k,3) * v(i,1,3,e) &
3568 + dz(k,4) * v(i,1,4,e)
3569
3570 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3571 + dz(k,2) * w(i,1,2,e) &
3572 + dz(k,3) * w(i,1,3,e) &
3573 + dz(k,4) * w(i,1,4,e)
3574 end do
3575 end do
3576
3577 do i = 1, lx*lx*lx
3578 ur(i,1,1) = h1(i,1,1,e) &
3579 * ( g11(i,1,1,e) * wur(i,1,1) &
3580 + g12(i,1,1,e) * wus(i,1,1) &
3581 + g13(i,1,1,e) * wut(i,1,1) )
3582 us(i,1,1) = h1(i,1,1,e) &
3583 * ( g12(i,1,1,e) * wur(i,1,1) &
3584 + g22(i,1,1,e) * wus(i,1,1) &
3585 + g23(i,1,1,e) * wut(i,1,1) )
3586 ut(i,1,1) = h1(i,1,1,e) &
3587 * ( g13(i,1,1,e) * wur(i,1,1) &
3588 + g23(i,1,1,e) * wus(i,1,1) &
3589 + g33(i,1,1,e) * wut(i,1,1) )
3590
3591 vr(i,1,1) = h1(i,1,1,e) &
3592 * ( g11(i,1,1,e) * wvr(i,1,1) &
3593 + g12(i,1,1,e) * wvs(i,1,1) &
3594 + g13(i,1,1,e) * wvt(i,1,1) )
3595 vs(i,1,1) = h1(i,1,1,e) &
3596 * ( g12(i,1,1,e) * wvr(i,1,1) &
3597 + g22(i,1,1,e) * wvs(i,1,1) &
3598 + g23(i,1,1,e) * wvt(i,1,1) )
3599 vt(i,1,1) = h1(i,1,1,e) &
3600 * ( g13(i,1,1,e) * wvr(i,1,1) &
3601 + g23(i,1,1,e) * wvs(i,1,1) &
3602 + g33(i,1,1,e) * wvt(i,1,1) )
3603
3604 wr(i,1,1) = h1(i,1,1,e) &
3605 * ( g11(i,1,1,e) * wwr(i,1,1) &
3606 + g12(i,1,1,e) * wws(i,1,1) &
3607 + g13(i,1,1,e) * wwt(i,1,1) )
3608 ws(i,1,1) = h1(i,1,1,e) &
3609 * ( g12(i,1,1,e) * wwr(i,1,1) &
3610 + g22(i,1,1,e) * wws(i,1,1) &
3611 + g23(i,1,1,e) * wwt(i,1,1) )
3612 wt(i,1,1) = h1(i,1,1,e) &
3613 * ( g13(i,1,1,e) * wwr(i,1,1) &
3614 + g23(i,1,1,e) * wws(i,1,1) &
3615 + g33(i,1,1,e) * wwt(i,1,1) )
3616 end do
3617
3618 do j = 1, lx*lx
3619 do i = 1, lx
3620 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
3621 + dxt(i,2) * ur(2,j,1) &
3622 + dxt(i,3) * ur(3,j,1) &
3623 + dxt(i,4) * ur(4,j,1)
3624
3625 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
3626 + dxt(i,2) * vr(2,j,1) &
3627 + dxt(i,3) * vr(3,j,1) &
3628 + dxt(i,4) * vr(4,j,1)
3629
3630 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
3631 + dxt(i,2) * wr(2,j,1) &
3632 + dxt(i,3) * wr(3,j,1) &
3633 + dxt(i,4) * wr(4,j,1)
3634 end do
3635 end do
3636
3637 do k = 1, lx
3638 do j = 1, lx
3639 do i = 1, lx
3640 au(i,j,k,e) = au(i,j,k,e) &
3641 + dyt(j,1) * us(i,1,k) &
3642 + dyt(j,2) * us(i,2,k) &
3643 + dyt(j,3) * us(i,3,k) &
3644 + dyt(j,4) * us(i,4,k)
3645
3646 av(i,j,k,e) = av(i,j,k,e) &
3647 + dyt(j,1) * vs(i,1,k) &
3648 + dyt(j,2) * vs(i,2,k) &
3649 + dyt(j,3) * vs(i,3,k) &
3650 + dyt(j,4) * vs(i,4,k)
3651
3652 aw(i,j,k,e) = aw(i,j,k,e) &
3653 + dyt(j,1) * ws(i,1,k) &
3654 + dyt(j,2) * ws(i,2,k) &
3655 + dyt(j,3) * ws(i,3,k) &
3656 + dyt(j,4) * ws(i,4,k)
3657 end do
3658 end do
3659 end do
3660
3661 do k = 1, lx
3662 do i = 1, lx*lx
3663 au(i,1,k,e) = au(i,1,k,e) &
3664 + dzt(k,1) * ut(i,1,1) &
3665 + dzt(k,2) * ut(i,1,2) &
3666 + dzt(k,3) * ut(i,1,3) &
3667 + dzt(k,4) * ut(i,1,4)
3668
3669 av(i,1,k,e) = av(i,1,k,e) &
3670 + dzt(k,1) * vt(i,1,1) &
3671 + dzt(k,2) * vt(i,1,2) &
3672 + dzt(k,3) * vt(i,1,3) &
3673 + dzt(k,4) * vt(i,1,4)
3674
3675 aw(i,1,k,e) = aw(i,1,k,e) &
3676 + dzt(k,1) * wt(i,1,1) &
3677 + dzt(k,2) * wt(i,1,2) &
3678 + dzt(k,3) * wt(i,1,3) &
3679 + dzt(k,4) * wt(i,1,4)
3680 end do
3681 end do
3682
3683 end do
3684 !$omp end do
3685 end subroutine ax_helm_vector_lx4
3686
3687 subroutine ax_helm_vector_lx3(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
3688 h1, G11, G22, G33, G12, G13, G23, n)
3689 integer, parameter :: lx = 3
3690 integer, intent(in) :: n
3691 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3692 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3693 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3694 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3695 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3696 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3697 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3698 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
3699 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
3700 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
3701 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
3702 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
3703 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
3704 real(kind=rp), intent(in) :: dx(lx, lx)
3705 real(kind=rp), intent(in) :: dy(lx, lx)
3706 real(kind=rp), intent(in) :: dz(lx, lx)
3707 real(kind=rp), intent(in) :: dxt(lx, lx)
3708 real(kind=rp), intent(in) :: dyt(lx, lx)
3709 real(kind=rp), intent(in) :: dzt(lx, lx)
3710 real(kind=rp) :: ur(lx, lx, lx)
3711 real(kind=rp) :: us(lx, lx, lx)
3712 real(kind=rp) :: ut(lx, lx, lx)
3713 real(kind=rp) :: vr(lx, lx, lx)
3714 real(kind=rp) :: vs(lx, lx, lx)
3715 real(kind=rp) :: vt(lx, lx, lx)
3716 real(kind=rp) :: wr(lx, lx, lx)
3717 real(kind=rp) :: ws(lx, lx, lx)
3718 real(kind=rp) :: wt(lx, lx, lx)
3719 real(kind=rp) :: wur(lx, lx, lx)
3720 real(kind=rp) :: wus(lx, lx, lx)
3721 real(kind=rp) :: wut(lx, lx, lx)
3722 real(kind=rp) :: wvr(lx, lx, lx)
3723 real(kind=rp) :: wvs(lx, lx, lx)
3724 real(kind=rp) :: wvt(lx, lx, lx)
3725 real(kind=rp) :: wwr(lx, lx, lx)
3726 real(kind=rp) :: wws(lx, lx, lx)
3727 real(kind=rp) :: wwt(lx, lx, lx)
3728 integer :: e, i, j, k
3729
3730 !$omp do
3731 do e = 1, n
3732 do j = 1, lx * lx
3733 do i = 1, lx
3734 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3735 + dx(i,2) * u(2,j,1,e) &
3736 + dx(i,3) * u(3,j,1,e)
3737
3738 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3739 + dx(i,2) * v(2,j,1,e) &
3740 + dx(i,3) * v(3,j,1,e)
3741
3742 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3743 + dx(i,2) * w(2,j,1,e) &
3744 + dx(i,3) * w(3,j,1,e)
3745 end do
3746 end do
3747
3748 do k = 1, lx
3749 do j = 1, lx
3750 do i = 1, lx
3751 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3752 + dy(j,2) * u(i,2,k,e) &
3753 + dy(j,3) * u(i,3,k,e)
3754
3755 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3756 + dy(j,2) * v(i,2,k,e) &
3757 + dy(j,3) * v(i,3,k,e)
3758
3759 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3760 + dy(j,2) * w(i,2,k,e) &
3761 + dy(j,3) * w(i,3,k,e)
3762 end do
3763 end do
3764 end do
3765
3766 do k = 1, lx
3767 do i = 1, lx*lx
3768 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3769 + dz(k,2) * u(i,1,2,e) &
3770 + dz(k,3) * u(i,1,3,e)
3771
3772 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3773 + dz(k,2) * v(i,1,2,e) &
3774 + dz(k,3) * v(i,1,3,e)
3775
3776 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3777 + dz(k,2) * w(i,1,2,e) &
3778 + dz(k,3) * w(i,1,3,e)
3779 end do
3780 end do
3781
3782 do i = 1, lx*lx*lx
3783 ur(i,1,1) = h1(i,1,1,e) &
3784 * ( g11(i,1,1,e) * wur(i,1,1) &
3785 + g12(i,1,1,e) * wus(i,1,1) &
3786 + g13(i,1,1,e) * wut(i,1,1) )
3787 us(i,1,1) = h1(i,1,1,e) &
3788 * ( g12(i,1,1,e) * wur(i,1,1) &
3789 + g22(i,1,1,e) * wus(i,1,1) &
3790 + g23(i,1,1,e) * wut(i,1,1) )
3791 ut(i,1,1) = h1(i,1,1,e) &
3792 * ( g13(i,1,1,e) * wur(i,1,1) &
3793 + g23(i,1,1,e) * wus(i,1,1) &
3794 + g33(i,1,1,e) * wut(i,1,1) )
3795
3796 vr(i,1,1) = h1(i,1,1,e) &
3797 * ( g11(i,1,1,e) * wvr(i,1,1) &
3798 + g12(i,1,1,e) * wvs(i,1,1) &
3799 + g13(i,1,1,e) * wvt(i,1,1) )
3800 vs(i,1,1) = h1(i,1,1,e) &
3801 * ( g12(i,1,1,e) * wvr(i,1,1) &
3802 + g22(i,1,1,e) * wvs(i,1,1) &
3803 + g23(i,1,1,e) * wvt(i,1,1) )
3804 vt(i,1,1) = h1(i,1,1,e) &
3805 * ( g13(i,1,1,e) * wwr(i,1,1) &
3806 + g23(i,1,1,e) * wws(i,1,1) &
3807 + g33(i,1,1,e) * wwt(i,1,1) )
3808
3809 wr(i,1,1) = h1(i,1,1,e) &
3810 * ( g11(i,1,1,e) * wwr(i,1,1) &
3811 + g12(i,1,1,e) * wws(i,1,1) &
3812 + g13(i,1,1,e) * wwt(i,1,1) )
3813 ws(i,1,1) = h1(i,1,1,e) &
3814 * ( g12(i,1,1,e) * wwr(i,1,1) &
3815 + g22(i,1,1,e) * wws(i,1,1) &
3816 + g23(i,1,1,e) * wwt(i,1,1) )
3817 wt(i,1,1) = h1(i,1,1,e) &
3818 * ( g13(i,1,1,e) * wwr(i,1,1) &
3819 + g23(i,1,1,e) * wws(i,1,1) &
3820 + g33(i,1,1,e) * wwt(i,1,1) )
3821 end do
3822
3823 do j = 1, lx*lx
3824 do i = 1, lx
3825 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
3826 + dxt(i,2) * ur(2,j,1) &
3827 + dxt(i,3) * ur(3,j,1)
3828
3829 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
3830 + dxt(i,2) * vr(2,j,1) &
3831 + dxt(i,3) * vr(3,j,1)
3832
3833 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
3834 + dxt(i,2) * wr(2,j,1) &
3835 + dxt(i,3) * wr(3,j,1)
3836 end do
3837 end do
3838
3839 do k = 1, lx
3840 do j = 1, lx
3841 do i = 1, lx
3842 au(i,j,k,e) = au(i,j,k,e) &
3843 + dyt(j,1) * us(i,1,k) &
3844 + dyt(j,2) * us(i,2,k) &
3845 + dyt(j,3) * us(i,3,k)
3846
3847 av(i,j,k,e) = av(i,j,k,e) &
3848 + dyt(j,1) * vs(i,1,k) &
3849 + dyt(j,2) * vs(i,2,k) &
3850 + dyt(j,3) * vs(i,3,k)
3851
3852 aw(i,j,k,e) = aw(i,j,k,e) &
3853 + dyt(j,1) * ws(i,1,k) &
3854 + dyt(j,2) * ws(i,2,k) &
3855 + dyt(j,3) * ws(i,3,k)
3856 end do
3857 end do
3858 end do
3859
3860 do k = 1, lx
3861 do i = 1, lx*lx
3862 au(i,1,k,e) = au(i,1,k,e) &
3863 + dzt(k,1) * ut(i,1,1) &
3864 + dzt(k,2) * ut(i,1,2) &
3865 + dzt(k,3) * ut(i,1,3)
3866
3867 av(i,1,k,e) = av(i,1,k,e) &
3868 + dzt(k,1) * vt(i,1,1) &
3869 + dzt(k,2) * vt(i,1,2) &
3870 + dzt(k,3) * vt(i,1,3)
3871
3872 aw(i,1,k,e) = aw(i,1,k,e) &
3873 + dzt(k,1) * wt(i,1,1) &
3874 + dzt(k,2) * wt(i,1,2) &
3875 + dzt(k,3) * wt(i,1,3)
3876 end do
3877 end do
3878
3879 end do
3880 !$omp end do
3881 end subroutine ax_helm_vector_lx3
3882
3883 subroutine ax_helm_vector_lx2(au, av, aw, u, v, w, Dx, Dy, Dz, Dxt, Dyt, Dzt, &
3884 h1, G11, G22, G33, G12, G13, G23, n)
3885 integer, parameter :: lx = 2
3886 integer, intent(in) :: n
3887 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
3888 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
3889 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
3890 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
3891 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
3892 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
3893 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
3894 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
3895 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
3896 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
3897 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
3898 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
3899 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
3900 real(kind=rp), intent(in) :: dx(lx, lx)
3901 real(kind=rp), intent(in) :: dy(lx, lx)
3902 real(kind=rp), intent(in) :: dz(lx, lx)
3903 real(kind=rp), intent(in) :: dxt(lx, lx)
3904 real(kind=rp), intent(in) :: dyt(lx, lx)
3905 real(kind=rp), intent(in) :: dzt(lx, lx)
3906 real(kind=rp) :: ur(lx, lx, lx)
3907 real(kind=rp) :: us(lx, lx, lx)
3908 real(kind=rp) :: ut(lx, lx, lx)
3909 real(kind=rp) :: vr(lx, lx, lx)
3910 real(kind=rp) :: vs(lx, lx, lx)
3911 real(kind=rp) :: vt(lx, lx, lx)
3912 real(kind=rp) :: wr(lx, lx, lx)
3913 real(kind=rp) :: ws(lx, lx, lx)
3914 real(kind=rp) :: wt(lx, lx, lx)
3915 real(kind=rp) :: wur(lx, lx, lx)
3916 real(kind=rp) :: wus(lx, lx, lx)
3917 real(kind=rp) :: wut(lx, lx, lx)
3918 real(kind=rp) :: wvr(lx, lx, lx)
3919 real(kind=rp) :: wvs(lx, lx, lx)
3920 real(kind=rp) :: wvt(lx, lx, lx)
3921 real(kind=rp) :: wwr(lx, lx, lx)
3922 real(kind=rp) :: wws(lx, lx, lx)
3923 real(kind=rp) :: wwt(lx, lx, lx)
3924 integer :: e, i, j, k
3925
3926 !$omp do
3927 do e = 1, n
3928 do j = 1, lx * lx
3929 do i = 1, lx
3930 wur(i,j,1) = dx(i,1) * u(1,j,1,e) &
3931 + dx(i,2) * u(2,j,1,e)
3932
3933 wvr(i,j,1) = dx(i,1) * v(1,j,1,e) &
3934 + dx(i,2) * v(2,j,1,e)
3935
3936 wwr(i,j,1) = dx(i,1) * w(1,j,1,e) &
3937 + dx(i,2) * w(2,j,1,e)
3938 end do
3939 end do
3940
3941 do k = 1, lx
3942 do j = 1, lx
3943 do i = 1, lx
3944 wus(i,j,k) = dy(j,1) * u(i,1,k,e) &
3945 + dy(j,2) * u(i,2,k,e)
3946
3947 wvs(i,j,k) = dy(j,1) * v(i,1,k,e) &
3948 + dy(j,2) * v(i,2,k,e)
3949
3950 wws(i,j,k) = dy(j,1) * w(i,1,k,e) &
3951 + dy(j,2) * w(i,2,k,e)
3952
3953 end do
3954 end do
3955 end do
3956
3957 do k = 1, lx
3958 do i = 1, lx*lx
3959 wut(i,1,k) = dz(k,1) * u(i,1,1,e) &
3960 + dz(k,2) * u(i,1,2,e)
3961
3962 wvt(i,1,k) = dz(k,1) * v(i,1,1,e) &
3963 + dz(k,2) * v(i,1,2,e)
3964
3965 wwt(i,1,k) = dz(k,1) * w(i,1,1,e) &
3966 + dz(k,2) * w(i,1,2,e)
3967 end do
3968 end do
3969
3970 do i = 1, lx*lx*lx
3971 ur(i,1,1) = h1(i,1,1,e) &
3972 * ( g11(i,1,1,e) * wur(i,1,1) &
3973 + g12(i,1,1,e) * wus(i,1,1) &
3974 + g13(i,1,1,e) * wut(i,1,1) )
3975 us(i,1,1) = h1(i,1,1,e) &
3976 * ( g12(i,1,1,e) * wur(i,1,1) &
3977 + g22(i,1,1,e) * wus(i,1,1) &
3978 + g23(i,1,1,e) * wut(i,1,1) )
3979 ut(i,1,1) = h1(i,1,1,e) &
3980 * ( g13(i,1,1,e) * wur(i,1,1) &
3981 + g23(i,1,1,e) * wus(i,1,1) &
3982 + g33(i,1,1,e) * wut(i,1,1) )
3983
3984 vr(i,1,1) = h1(i,1,1,e) &
3985 * ( g11(i,1,1,e) * wvr(i,1,1) &
3986 + g12(i,1,1,e) * wvs(i,1,1) &
3987 + g13(i,1,1,e) * wvt(i,1,1) )
3988 vs(i,1,1) = h1(i,1,1,e) &
3989 * ( g12(i,1,1,e) * wvr(i,1,1) &
3990 + g22(i,1,1,e) * wvs(i,1,1) &
3991 + g23(i,1,1,e) * wvt(i,1,1) )
3992 vt(i,1,1) = h1(i,1,1,e) &
3993 * ( g13(i,1,1,e) * wvr(i,1,1) &
3994 + g23(i,1,1,e) * wvs(i,1,1) &
3995 + g33(i,1,1,e) * wvt(i,1,1) )
3996
3997 wr(i,1,1) = h1(i,1,1,e) &
3998 * ( g11(i,1,1,e) * wwr(i,1,1) &
3999 + g12(i,1,1,e) * wws(i,1,1) &
4000 + g13(i,1,1,e) * wwt(i,1,1) )
4001 ws(i,1,1) = h1(i,1,1,e) &
4002 * ( g12(i,1,1,e) * wwr(i,1,1) &
4003 + g22(i,1,1,e) * wws(i,1,1) &
4004 + g23(i,1,1,e) * wwt(i,1,1) )
4005 wt(i,1,1) = h1(i,1,1,e) &
4006 * ( g13(i,1,1,e) * wwr(i,1,1) &
4007 + g23(i,1,1,e) * wws(i,1,1) &
4008 + g33(i,1,1,e) * wwt(i,1,1) )
4009 end do
4010
4011 do j = 1, lx*lx
4012 do i = 1, lx
4013 au(i,j,1,e) = dxt(i,1) * ur(1,j,1) &
4014 + dxt(i,2) * ur(2,j,1)
4015
4016 av(i,j,1,e) = dxt(i,1) * vr(1,j,1) &
4017 + dxt(i,2) * vr(2,j,1)
4018
4019 aw(i,j,1,e) = dxt(i,1) * wr(1,j,1) &
4020 + dxt(i,2) * wr(2,j,1)
4021 end do
4022 end do
4023
4024 do k = 1, lx
4025 do j = 1, lx
4026 do i = 1, lx
4027 au(i,j,k,e) = au(i,j,k,e) &
4028 + dyt(j,1) * us(i,1,k) &
4029 + dyt(j,2) * us(i,2,k)
4030
4031 av(i,j,k,e) = av(i,j,k,e) &
4032 + dyt(j,1) * vs(i,1,k) &
4033 + dyt(j,2) * vs(i,2,k)
4034
4035 aw(i,j,k,e) = aw(i,j,k,e) &
4036 + dyt(j,1) * ws(i,1,k) &
4037 + dyt(j,2) * ws(i,2,k)
4038 end do
4039 end do
4040 end do
4041
4042 do k = 1, lx
4043 do i = 1, lx*lx
4044 au(i,1,k,e) = au(i,1,k,e) &
4045 + dzt(k,1) * ut(i,1,1) &
4046 + dzt(k,2) * ut(i,1,2)
4047
4048 av(i,1,k,e) = av(i,1,k,e) &
4049 + dzt(k,1) * vt(i,1,1) &
4050 + dzt(k,2) * vt(i,1,2)
4051
4052 aw(i,1,k,e) = aw(i,1,k,e) &
4053 + dzt(k,1) * wt(i,1,1) &
4054 + dzt(k,2) * wt(i,1,2)
4055 end do
4056 end do
4057
4058 end do
4059 !$omp end do
4060 end subroutine ax_helm_vector_lx2
4061
4062end submodule ax_helm_vector_cpu
ax_helm_cpu
Definition ax_helm_cpu.f90:33