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