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