Neko  0.9.99
A portable framework for high-order spectral element flow simulations
cpu_opgrad.f90
Go to the documentation of this file.
1 ! Copyright (c) 2021-2024, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
34 submodule(opr_cpu) cpu_opgrad
35  implicit none
36 
37 contains
38 
39  module subroutine opr_cpu_opgrad(ux, uy, uz, u, coef, e_start, e_end)
40  type(coef_t), intent(in) :: coef
41  integer, intent(in) :: e_start, e_end
42  real(kind=rp), dimension(coef%Xh%lxyz, e_end - e_start + 1), &
43  intent(inout) :: ux
44  real(kind=rp), dimension(coef%Xh%lxyz, e_end - e_start + 1), &
45  intent(inout) :: uy
46  real(kind=rp), dimension(coef%Xh%lxyz, e_end - e_start + 1), &
47  intent(inout) :: uz
48  real(kind=rp), dimension(coef%Xh%lxyz, e_end - e_start + 1), &
49  intent(in) :: u
50  integer :: e_len
51 
52  e_len = e_end-e_start+1
53 
54  if (e_len .eq. 1) then
55  call opr_cpu_opgrad_single(ux, uy, uz, u, coef, e_start)
56  else
57  call opr_cpu_opgrad_many(ux, uy, uz, u, coef, e_start, e_len)
58  end if
59  end subroutine opr_cpu_opgrad
60 
61  subroutine opr_cpu_opgrad_many(ux, uy, uz, u, coef, e_start, e_len)
62  type(coef_t), intent(in) :: coef
63  integer, intent(in) :: e_start, e_len
64  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(inout) :: ux
65  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(inout) :: uy
66  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(inout) :: uz
67  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(in) :: u
68 
69  associate(xh => coef%Xh, msh => coef%msh, &
70  drdx => coef%drdx, drdy => coef%drdy, drdz => coef%drdz, &
71  dsdx => coef%dsdx, dsdy => coef%dsdy, dsdz => coef%dsdz, &
72  dtdx => coef%dtdx, dtdy => coef%dtdy, dtdz => coef%dtdz)
73 
74  select case (xh%lx)
75  case (18)
76  call cpu_opgrad_lx18(ux, uy, uz, u, &
77  xh%dx, xh%dy, xh%dz, &
78  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
79  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
80  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
81  xh%w3, e_len)
82  case (17)
83  call cpu_opgrad_lx17(ux, uy, uz, u, &
84  xh%dx, xh%dy, xh%dz, &
85  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
86  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
87  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
88  xh%w3, e_len)
89  case (16)
90  call cpu_opgrad_lx16(ux, uy, uz, u, &
91  xh%dx, xh%dy, xh%dz, &
92  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
93  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
94  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
95  xh%w3, e_len)
96  case (15)
97  call cpu_opgrad_lx15(ux, uy, uz, u, &
98  xh%dx, xh%dy, xh%dz, &
99  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
100  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
101  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
102  xh%w3, e_len)
103  case (14)
104  call cpu_opgrad_lx14(ux, uy, uz, u, &
105  xh%dx, xh%dy, xh%dz, &
106  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
107  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
108  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
109  xh%w3, e_len)
110  case (13)
111  call cpu_opgrad_lx13(ux, uy, uz, u, &
112  xh%dx, xh%dy, xh%dz, &
113  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
114  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
115  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
116  xh%w3, e_len)
117  case (12)
118  call cpu_opgrad_lx12(ux, uy, uz, u, &
119  xh%dx, xh%dy, xh%dz, &
120  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
121  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
122  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
123  xh%w3, e_len)
124  case (11)
125  call cpu_opgrad_lx11(ux, uy, uz, u, &
126  xh%dx, xh%dy, xh%dz, &
127  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
128  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
129  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
130  xh%w3, e_len)
131  case (10)
132  call cpu_opgrad_lx10(ux, uy, uz, u, &
133  xh%dx, xh%dy, xh%dz, &
134  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
135  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
136  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
137  xh%w3, e_len)
138 
139  case (9)
140  call cpu_opgrad_lx9(ux, uy, uz, u, &
141  xh%dx, xh%dy, xh%dz, &
142  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
143  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
144  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
145  xh%w3, e_len)
146  case (8)
147  call cpu_opgrad_lx8(ux, uy, uz, u, &
148  xh%dx, xh%dy, xh%dz, &
149  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
150  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
151  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
152  xh%w3, e_len)
153  case (7)
154  call cpu_opgrad_lx7(ux, uy, uz, u, &
155  xh%dx, xh%dy, xh%dz, &
156  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
157  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
158  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
159  xh%w3, e_len)
160  case (6)
161  call cpu_opgrad_lx6(ux, uy, uz, u, &
162  xh%dx, xh%dy, xh%dz, &
163  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
164  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
165  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
166  xh%w3, e_len)
167  case (5)
168  call cpu_opgrad_lx5(ux, uy, uz, u, &
169  xh%dx, xh%dy, xh%dz, &
170  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
171  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
172  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
173  xh%w3, e_len)
174  case (4)
175  call cpu_opgrad_lx4(ux, uy, uz, u, &
176  xh%dx, xh%dy, xh%dz, &
177  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
178  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
179  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
180  xh%w3, e_len)
181  case (3)
182  call cpu_opgrad_lx3(ux, uy, uz, u, &
183  xh%dx, xh%dy, xh%dz, &
184  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
185  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
186  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
187  xh%w3, e_len)
188  case (2)
189  call cpu_opgrad_lx2(ux, uy, uz, u, &
190  xh%dx, xh%dy, xh%dz, &
191  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
192  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
193  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
194  xh%w3, e_len)
195  case default
196  call cpu_opgrad_lx(ux, uy, uz, u, &
197  xh%dx, xh%dy, xh%dz, &
198  drdx(1,1,1,e_start), dsdx(1,1,1,e_start), dtdx(1,1,1,e_start), &
199  drdy(1,1,1,e_start), dsdy(1,1,1,e_start), dtdy(1,1,1,e_start), &
200  drdz(1,1,1,e_start), dsdz(1,1,1,e_start), dtdz(1,1,1,e_start), &
201  xh%w3, e_len, xh%lx)
202  end select
203  end associate
204 
205  end subroutine opr_cpu_opgrad_many
206 
207  subroutine cpu_opgrad_lx(ux, uy, uz, u, dx, dy, dz, &
208  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n, lx)
209  integer, intent(in) :: n, lx
210  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
211  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
212  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
213  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
214  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
215  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
216  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
217  real(kind=rp) :: ur(lx, lx, lx)
218  real(kind=rp) :: us(lx, lx, lx)
219  real(kind=rp) :: ut(lx, lx, lx)
220  real(kind=rp) :: tmp
221  integer :: e, i, j, k, l
222 
223  do e = 1, n
224  do j = 1, lx * lx
225  do i = 1, lx
226  tmp = 0.0_rp
227  do k = 1, lx
228  tmp = tmp + dx(i,k) * u(k,j,1,e)
229  end do
230  ur(i,j,1) = tmp
231  end do
232  end do
233 
234  do k = 1, lx
235  do j = 1, lx
236  do i = 1, lx
237  tmp = 0.0_rp
238  do l = 1, lx
239  tmp = tmp + dy(j,l) * u(i,l,k,e)
240  end do
241  us(i,j,k) = tmp
242  end do
243  end do
244  end do
245 
246  do k = 1, lx
247  do i = 1, lx*lx
248  tmp = 0.0_rp
249  do l = 1, lx
250  tmp = tmp + dz(k,l) * u(i,1,l,e)
251  end do
252  ut(i,1,k) = tmp
253  end do
254  end do
255 
256  do i = 1, lx * lx * lx
257  ux(i,1,1,e) = w3(i,1,1) &
258  * ( drdx(i,1,1,e) * ur(i,1,1) &
259  + dsdx(i,1,1,e) * us(i,1,1) &
260  + dtdx(i,1,1,e) * ut(i,1,1) )
261  uy(i,1,1,e) = w3(i,1,1) &
262  * ( dsdy(i,1,1,e) * us(i,1,1) &
263  + drdy(i,1,1,e) * ur(i,1,1) &
264  + dtdy(i,1,1,e) * ut(i,1,1) )
265  uz(i,1,1,e) = w3(i,1,1) &
266  * ( dtdz(i,1,1,e) * ut(i,1,1) &
267  + drdz(i,1,1,e) * ur(i,1,1) &
268  + dsdz(i,1,1,e) * us(i,1,1) )
269  end do
270  end do
271  end subroutine cpu_opgrad_lx
272 
273  subroutine cpu_opgrad_lx18(ux, uy, uz, u, dx, dy, dz, &
274  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
275  integer, parameter :: lx = 18
276  integer, intent(in) :: n
277  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
278  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
279  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
280  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
281  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
282  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
283  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
284  real(kind=rp) :: ur(lx, lx, lx)
285  real(kind=rp) :: us(lx, lx, lx)
286  real(kind=rp) :: ut(lx, lx, lx)
287  integer :: e, i, j, k
288 
289  do e = 1, n
290  do j = 1, lx * lx
291  do i = 1, lx
292  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
293  + dx(i,2) * u(2,j,1,e) &
294  + dx(i,3) * u(3,j,1,e) &
295  + dx(i,4) * u(4,j,1,e) &
296  + dx(i,5) * u(5,j,1,e) &
297  + dx(i,6) * u(6,j,1,e) &
298  + dx(i,7) * u(7,j,1,e) &
299  + dx(i,8) * u(8,j,1,e) &
300  + dx(i,9) * u(9,j,1,e) &
301  + dx(i,10) * u(10,j,1,e) &
302  + dx(i,11) * u(11,j,1,e) &
303  + dx(i,12) * u(12,j,1,e) &
304  + dx(i,13) * u(13,j,1,e) &
305  + dx(i,14) * u(14,j,1,e) &
306  + dx(i,15) * u(15,j,1,e) &
307  + dx(i,16) * u(16,j,1,e) &
308  + dx(i,17) * u(17,j,1,e) &
309  + dx(i,18) * u(18,j,1,e)
310  end do
311  end do
312 
313  do k = 1, lx
314  do j = 1, lx
315  do i = 1, lx
316  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
317  + dy(j,2) * u(i,2,k,e) &
318  + dy(j,3) * u(i,3,k,e) &
319  + dy(j,4) * u(i,4,k,e) &
320  + dy(j,5) * u(i,5,k,e) &
321  + dy(j,6) * u(i,6,k,e) &
322  + dy(j,7) * u(i,7,k,e) &
323  + dy(j,8) * u(i,8,k,e) &
324  + dy(j,9) * u(i,9,k,e) &
325  + dy(j,10) * u(i,10,k,e) &
326  + dy(j,11) * u(i,11,k,e) &
327  + dy(j,12) * u(i,12,k,e) &
328  + dy(j,13) * u(i,13,k,e) &
329  + dy(j,14) * u(i,14,k,e) &
330  + dy(j,15) * u(i,15,k,e) &
331  + dy(j,16) * u(i,16,k,e) &
332  + dy(j,17) * u(i,17,k,e) &
333  + dy(j,18) * u(i,18,k,e)
334  end do
335  end do
336  end do
337 
338  do k = 1, lx
339  do i = 1, lx*lx
340  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
341  + dz(k,2) * u(i,1,2,e) &
342  + dz(k,3) * u(i,1,3,e) &
343  + dz(k,4) * u(i,1,4,e) &
344  + dz(k,5) * u(i,1,5,e) &
345  + dz(k,6) * u(i,1,6,e) &
346  + dz(k,7) * u(i,1,7,e) &
347  + dz(k,8) * u(i,1,8,e) &
348  + dz(k,9) * u(i,1,9,e) &
349  + dz(k,10) * u(i,1,10,e) &
350  + dz(k,11) * u(i,1,11,e) &
351  + dz(k,12) * u(i,1,12,e) &
352  + dz(k,13) * u(i,1,13,e) &
353  + dz(k,14) * u(i,1,14,e) &
354  + dz(k,15) * u(i,1,15,e) &
355  + dz(k,16) * u(i,1,16,e) &
356  + dz(k,17) * u(i,1,17,e) &
357  + dz(k,18) * u(i,1,18,e)
358  end do
359  end do
360 
361  do i = 1, lx * lx * lx
362  ux(i,1,1,e) = w3(i,1,1) &
363  * ( drdx(i,1,1,e) * ur(i,1,1) &
364  + dsdx(i,1,1,e) * us(i,1,1) &
365  + dtdx(i,1,1,e) * ut(i,1,1) )
366  uy(i,1,1,e) = w3(i,1,1) &
367  * ( dsdy(i,1,1,e) * us(i,1,1) &
368  + drdy(i,1,1,e) * ur(i,1,1) &
369  + dtdy(i,1,1,e) * ut(i,1,1) )
370  uz(i,1,1,e) = w3(i,1,1) &
371  * ( dtdz(i,1,1,e) * ut(i,1,1) &
372  + drdz(i,1,1,e) * ur(i,1,1) &
373  + dsdz(i,1,1,e) * us(i,1,1) )
374  end do
375  end do
376  end subroutine cpu_opgrad_lx18
377 
378  subroutine cpu_opgrad_lx17(ux, uy, uz, u, dx, dy, dz, &
379  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
380  integer, parameter :: lx = 17
381  integer, intent(in) :: n
382  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
383  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
384  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
385  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
386  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
387  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
388  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
389  real(kind=rp) :: ur(lx, lx, lx)
390  real(kind=rp) :: us(lx, lx, lx)
391  real(kind=rp) :: ut(lx, lx, lx)
392  integer :: e, i, j, k
393 
394  do e = 1, n
395  do j = 1, lx * lx
396  do i = 1, lx
397  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
398  + dx(i,2) * u(2,j,1,e) &
399  + dx(i,3) * u(3,j,1,e) &
400  + dx(i,4) * u(4,j,1,e) &
401  + dx(i,5) * u(5,j,1,e) &
402  + dx(i,6) * u(6,j,1,e) &
403  + dx(i,7) * u(7,j,1,e) &
404  + dx(i,8) * u(8,j,1,e) &
405  + dx(i,9) * u(9,j,1,e) &
406  + dx(i,10) * u(10,j,1,e) &
407  + dx(i,11) * u(11,j,1,e) &
408  + dx(i,12) * u(12,j,1,e) &
409  + dx(i,13) * u(13,j,1,e) &
410  + dx(i,14) * u(14,j,1,e) &
411  + dx(i,15) * u(15,j,1,e) &
412  + dx(i,16) * u(16,j,1,e) &
413  + dx(i,17) * u(17,j,1,e)
414  end do
415  end do
416 
417  do k = 1, lx
418  do j = 1, lx
419  do i = 1, lx
420  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
421  + dy(j,2) * u(i,2,k,e) &
422  + dy(j,3) * u(i,3,k,e) &
423  + dy(j,4) * u(i,4,k,e) &
424  + dy(j,5) * u(i,5,k,e) &
425  + dy(j,6) * u(i,6,k,e) &
426  + dy(j,7) * u(i,7,k,e) &
427  + dy(j,8) * u(i,8,k,e) &
428  + dy(j,9) * u(i,9,k,e) &
429  + dy(j,10) * u(i,10,k,e) &
430  + dy(j,11) * u(i,11,k,e) &
431  + dy(j,12) * u(i,12,k,e) &
432  + dy(j,13) * u(i,13,k,e) &
433  + dy(j,14) * u(i,14,k,e) &
434  + dy(j,15) * u(i,15,k,e) &
435  + dy(j,16) * u(i,16,k,e) &
436  + dy(j,17) * u(i,17,k,e)
437  end do
438  end do
439  end do
440 
441  do k = 1, lx
442  do i = 1, lx*lx
443  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
444  + dz(k,2) * u(i,1,2,e) &
445  + dz(k,3) * u(i,1,3,e) &
446  + dz(k,4) * u(i,1,4,e) &
447  + dz(k,5) * u(i,1,5,e) &
448  + dz(k,6) * u(i,1,6,e) &
449  + dz(k,7) * u(i,1,7,e) &
450  + dz(k,8) * u(i,1,8,e) &
451  + dz(k,9) * u(i,1,9,e) &
452  + dz(k,10) * u(i,1,10,e) &
453  + dz(k,11) * u(i,1,11,e) &
454  + dz(k,12) * u(i,1,12,e) &
455  + dz(k,13) * u(i,1,13,e) &
456  + dz(k,14) * u(i,1,14,e) &
457  + dz(k,15) * u(i,1,15,e) &
458  + dz(k,16) * u(i,1,16,e) &
459  + dz(k,17) * u(i,1,17,e)
460  end do
461  end do
462 
463  do i = 1, lx * lx * lx
464  ux(i,1,1,e) = w3(i,1,1) &
465  * ( drdx(i,1,1,e) * ur(i,1,1) &
466  + dsdx(i,1,1,e) * us(i,1,1) &
467  + dtdx(i,1,1,e) * ut(i,1,1) )
468  uy(i,1,1,e) = w3(i,1,1) &
469  * ( dsdy(i,1,1,e) * us(i,1,1) &
470  + drdy(i,1,1,e) * ur(i,1,1) &
471  + dtdy(i,1,1,e) * ut(i,1,1) )
472  uz(i,1,1,e) = w3(i,1,1) &
473  * ( dtdz(i,1,1,e) * ut(i,1,1) &
474  + drdz(i,1,1,e) * ur(i,1,1) &
475  + dsdz(i,1,1,e) * us(i,1,1) )
476  end do
477  end do
478  end subroutine cpu_opgrad_lx17
479 
480  subroutine cpu_opgrad_lx16(ux, uy, uz, u, dx, dy, dz, &
481  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
482  integer, parameter :: lx = 16
483  integer, intent(in) :: n
484  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
485  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
486  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
487  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
488  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
489  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
490  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
491  real(kind=rp) :: ur(lx, lx, lx)
492  real(kind=rp) :: us(lx, lx, lx)
493  real(kind=rp) :: ut(lx, lx, lx)
494  integer :: e, i, j, k
495 
496  do e = 1, n
497  do j = 1, lx * lx
498  do i = 1, lx
499  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
500  + dx(i,2) * u(2,j,1,e) &
501  + dx(i,3) * u(3,j,1,e) &
502  + dx(i,4) * u(4,j,1,e) &
503  + dx(i,5) * u(5,j,1,e) &
504  + dx(i,6) * u(6,j,1,e) &
505  + dx(i,7) * u(7,j,1,e) &
506  + dx(i,8) * u(8,j,1,e) &
507  + dx(i,9) * u(9,j,1,e) &
508  + dx(i,10) * u(10,j,1,e) &
509  + dx(i,11) * u(11,j,1,e) &
510  + dx(i,12) * u(12,j,1,e) &
511  + dx(i,13) * u(13,j,1,e) &
512  + dx(i,14) * u(14,j,1,e) &
513  + dx(i,15) * u(15,j,1,e) &
514  + dx(i,16) * u(16,j,1,e)
515  end do
516  end do
517 
518  do k = 1, lx
519  do j = 1, lx
520  do i = 1, lx
521  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
522  + dy(j,2) * u(i,2,k,e) &
523  + dy(j,3) * u(i,3,k,e) &
524  + dy(j,4) * u(i,4,k,e) &
525  + dy(j,5) * u(i,5,k,e) &
526  + dy(j,6) * u(i,6,k,e) &
527  + dy(j,7) * u(i,7,k,e) &
528  + dy(j,8) * u(i,8,k,e) &
529  + dy(j,9) * u(i,9,k,e) &
530  + dy(j,10) * u(i,10,k,e) &
531  + dy(j,11) * u(i,11,k,e) &
532  + dy(j,12) * u(i,12,k,e) &
533  + dy(j,13) * u(i,13,k,e) &
534  + dy(j,14) * u(i,14,k,e) &
535  + dy(j,15) * u(i,15,k,e) &
536  + dy(j,16) * u(i,16,k,e)
537  end do
538  end do
539  end do
540 
541  do k = 1, lx
542  do i = 1, lx*lx
543  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
544  + dz(k,2) * u(i,1,2,e) &
545  + dz(k,3) * u(i,1,3,e) &
546  + dz(k,4) * u(i,1,4,e) &
547  + dz(k,5) * u(i,1,5,e) &
548  + dz(k,6) * u(i,1,6,e) &
549  + dz(k,7) * u(i,1,7,e) &
550  + dz(k,8) * u(i,1,8,e) &
551  + dz(k,9) * u(i,1,9,e) &
552  + dz(k,10) * u(i,1,10,e) &
553  + dz(k,11) * u(i,1,11,e) &
554  + dz(k,12) * u(i,1,12,e) &
555  + dz(k,13) * u(i,1,13,e) &
556  + dz(k,14) * u(i,1,14,e) &
557  + dz(k,15) * u(i,1,15,e) &
558  + dz(k,16) * u(i,1,16,e)
559  end do
560  end do
561 
562  do i = 1, lx * lx * lx
563  ux(i,1,1,e) = w3(i,1,1) &
564  * ( drdx(i,1,1,e) * ur(i,1,1) &
565  + dsdx(i,1,1,e) * us(i,1,1) &
566  + dtdx(i,1,1,e) * ut(i,1,1) )
567  uy(i,1,1,e) = w3(i,1,1) &
568  * ( dsdy(i,1,1,e) * us(i,1,1) &
569  + drdy(i,1,1,e) * ur(i,1,1) &
570  + dtdy(i,1,1,e) * ut(i,1,1) )
571  uz(i,1,1,e) = w3(i,1,1) &
572  * ( dtdz(i,1,1,e) * ut(i,1,1) &
573  + drdz(i,1,1,e) * ur(i,1,1) &
574  + dsdz(i,1,1,e) * us(i,1,1) )
575  end do
576  end do
577  end subroutine cpu_opgrad_lx16
578 
579  subroutine cpu_opgrad_lx15(ux, uy, uz, u, dx, dy, dz, &
580  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
581  integer, parameter :: lx = 15
582  integer, intent(in) :: n
583  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
584  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
585  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
586  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
587  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
588  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
589  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
590  real(kind=rp) :: ur(lx, lx, lx)
591  real(kind=rp) :: us(lx, lx, lx)
592  real(kind=rp) :: ut(lx, lx, lx)
593  integer :: e, i, j, k
594 
595  do e = 1, n
596  do j = 1, lx * lx
597  do i = 1, lx
598  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
599  + dx(i,2) * u(2,j,1,e) &
600  + dx(i,3) * u(3,j,1,e) &
601  + dx(i,4) * u(4,j,1,e) &
602  + dx(i,5) * u(5,j,1,e) &
603  + dx(i,6) * u(6,j,1,e) &
604  + dx(i,7) * u(7,j,1,e) &
605  + dx(i,8) * u(8,j,1,e) &
606  + dx(i,9) * u(9,j,1,e) &
607  + dx(i,10) * u(10,j,1,e) &
608  + dx(i,11) * u(11,j,1,e) &
609  + dx(i,12) * u(12,j,1,e) &
610  + dx(i,13) * u(13,j,1,e) &
611  + dx(i,14) * u(14,j,1,e) &
612  + dx(i,15) * u(15,j,1,e)
613  end do
614  end do
615 
616  do k = 1, lx
617  do j = 1, lx
618  do i = 1, lx
619  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
620  + dy(j,2) * u(i,2,k,e) &
621  + dy(j,3) * u(i,3,k,e) &
622  + dy(j,4) * u(i,4,k,e) &
623  + dy(j,5) * u(i,5,k,e) &
624  + dy(j,6) * u(i,6,k,e) &
625  + dy(j,7) * u(i,7,k,e) &
626  + dy(j,8) * u(i,8,k,e) &
627  + dy(j,9) * u(i,9,k,e) &
628  + dy(j,10) * u(i,10,k,e) &
629  + dy(j,11) * u(i,11,k,e) &
630  + dy(j,12) * u(i,12,k,e) &
631  + dy(j,13) * u(i,13,k,e) &
632  + dy(j,14) * u(i,14,k,e) &
633  + dy(j,15) * u(i,15,k,e)
634  end do
635  end do
636  end do
637 
638  do k = 1, lx
639  do i = 1, lx*lx
640  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
641  + dz(k,2) * u(i,1,2,e) &
642  + dz(k,3) * u(i,1,3,e) &
643  + dz(k,4) * u(i,1,4,e) &
644  + dz(k,5) * u(i,1,5,e) &
645  + dz(k,6) * u(i,1,6,e) &
646  + dz(k,7) * u(i,1,7,e) &
647  + dz(k,8) * u(i,1,8,e) &
648  + dz(k,9) * u(i,1,9,e) &
649  + dz(k,10) * u(i,1,10,e) &
650  + dz(k,11) * u(i,1,11,e) &
651  + dz(k,12) * u(i,1,12,e) &
652  + dz(k,13) * u(i,1,13,e) &
653  + dz(k,14) * u(i,1,14,e) &
654  + dz(k,15) * u(i,1,15,e)
655  end do
656  end do
657 
658  do i = 1, lx * lx * lx
659  ux(i,1,1,e) = w3(i,1,1) &
660  * ( drdx(i,1,1,e) * ur(i,1,1) &
661  + dsdx(i,1,1,e) * us(i,1,1) &
662  + dtdx(i,1,1,e) * ut(i,1,1) )
663  uy(i,1,1,e) = w3(i,1,1) &
664  * ( dsdy(i,1,1,e) * us(i,1,1) &
665  + drdy(i,1,1,e) * ur(i,1,1) &
666  + dtdy(i,1,1,e) * ut(i,1,1) )
667  uz(i,1,1,e) = w3(i,1,1) &
668  * ( dtdz(i,1,1,e) * ut(i,1,1) &
669  + drdz(i,1,1,e) * ur(i,1,1) &
670  + dsdz(i,1,1,e) * us(i,1,1) )
671  end do
672  end do
673  end subroutine cpu_opgrad_lx15
674 
675  subroutine cpu_opgrad_lx14(ux, uy, uz, u, dx, dy, dz, &
676  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
677  integer, parameter :: lx = 14
678  integer, intent(in) :: n
679  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
680  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
681  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
682  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
683  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
684  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
685  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
686  real(kind=rp) :: ur(lx, lx, lx)
687  real(kind=rp) :: us(lx, lx, lx)
688  real(kind=rp) :: ut(lx, lx, lx)
689  integer :: e, i, j, k
690 
691  do e = 1, n
692  do j = 1, lx * lx
693  do i = 1, lx
694  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
695  + dx(i,2) * u(2,j,1,e) &
696  + dx(i,3) * u(3,j,1,e) &
697  + dx(i,4) * u(4,j,1,e) &
698  + dx(i,5) * u(5,j,1,e) &
699  + dx(i,6) * u(6,j,1,e) &
700  + dx(i,7) * u(7,j,1,e) &
701  + dx(i,8) * u(8,j,1,e) &
702  + dx(i,9) * u(9,j,1,e) &
703  + dx(i,10) * u(10,j,1,e) &
704  + dx(i,11) * u(11,j,1,e) &
705  + dx(i,12) * u(12,j,1,e) &
706  + dx(i,13) * u(13,j,1,e) &
707  + dx(i,14) * u(14,j,1,e)
708  end do
709  end do
710 
711  do k = 1, lx
712  do j = 1, lx
713  do i = 1, lx
714  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
715  + dy(j,2) * u(i,2,k,e) &
716  + dy(j,3) * u(i,3,k,e) &
717  + dy(j,4) * u(i,4,k,e) &
718  + dy(j,5) * u(i,5,k,e) &
719  + dy(j,6) * u(i,6,k,e) &
720  + dy(j,7) * u(i,7,k,e) &
721  + dy(j,8) * u(i,8,k,e) &
722  + dy(j,9) * u(i,9,k,e) &
723  + dy(j,10) * u(i,10,k,e) &
724  + dy(j,11) * u(i,11,k,e) &
725  + dy(j,12) * u(i,12,k,e) &
726  + dy(j,13) * u(i,13,k,e) &
727  + dy(j,14) * u(i,14,k,e)
728  end do
729  end do
730  end do
731 
732  do k = 1, lx
733  do i = 1, lx*lx
734  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
735  + dz(k,2) * u(i,1,2,e) &
736  + dz(k,3) * u(i,1,3,e) &
737  + dz(k,4) * u(i,1,4,e) &
738  + dz(k,5) * u(i,1,5,e) &
739  + dz(k,6) * u(i,1,6,e) &
740  + dz(k,7) * u(i,1,7,e) &
741  + dz(k,8) * u(i,1,8,e) &
742  + dz(k,9) * u(i,1,9,e) &
743  + dz(k,10) * u(i,1,10,e) &
744  + dz(k,11) * u(i,1,11,e) &
745  + dz(k,12) * u(i,1,12,e) &
746  + dz(k,13) * u(i,1,13,e) &
747  + dz(k,14) * u(i,1,14,e)
748  end do
749  end do
750 
751  do i = 1, lx * lx * lx
752  ux(i,1,1,e) = w3(i,1,1) &
753  * ( drdx(i,1,1,e) * ur(i,1,1) &
754  + dsdx(i,1,1,e) * us(i,1,1) &
755  + dtdx(i,1,1,e) * ut(i,1,1) )
756  uy(i,1,1,e) = w3(i,1,1) &
757  * ( dsdy(i,1,1,e) * us(i,1,1) &
758  + drdy(i,1,1,e) * ur(i,1,1) &
759  + dtdy(i,1,1,e) * ut(i,1,1) )
760  uz(i,1,1,e) = w3(i,1,1) &
761  * ( dtdz(i,1,1,e) * ut(i,1,1) &
762  + drdz(i,1,1,e) * ur(i,1,1) &
763  + dsdz(i,1,1,e) * us(i,1,1) )
764  end do
765  end do
766  end subroutine cpu_opgrad_lx14
767 
768  subroutine cpu_opgrad_lx13(ux, uy, uz, u, dx, dy, dz, &
769  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
770  integer, parameter :: lx = 13
771  integer, intent(in) :: n
772  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
773  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
774  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
775  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
776  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
777  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
778  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
779  real(kind=rp) :: ur(lx, lx, lx)
780  real(kind=rp) :: us(lx, lx, lx)
781  real(kind=rp) :: ut(lx, lx, lx)
782  integer :: e, i, j, k
783 
784  do e = 1, n
785  do j = 1, lx * lx
786  do i = 1, lx
787  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
788  + dx(i,2) * u(2,j,1,e) &
789  + dx(i,3) * u(3,j,1,e) &
790  + dx(i,4) * u(4,j,1,e) &
791  + dx(i,5) * u(5,j,1,e) &
792  + dx(i,6) * u(6,j,1,e) &
793  + dx(i,7) * u(7,j,1,e) &
794  + dx(i,8) * u(8,j,1,e) &
795  + dx(i,9) * u(9,j,1,e) &
796  + dx(i,10) * u(10,j,1,e) &
797  + dx(i,11) * u(11,j,1,e) &
798  + dx(i,12) * u(12,j,1,e) &
799  + dx(i,13) * u(13,j,1,e)
800  end do
801  end do
802 
803  do k = 1, lx
804  do j = 1, lx
805  do i = 1, lx
806  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
807  + dy(j,2) * u(i,2,k,e) &
808  + dy(j,3) * u(i,3,k,e) &
809  + dy(j,4) * u(i,4,k,e) &
810  + dy(j,5) * u(i,5,k,e) &
811  + dy(j,6) * u(i,6,k,e) &
812  + dy(j,7) * u(i,7,k,e) &
813  + dy(j,8) * u(i,8,k,e) &
814  + dy(j,9) * u(i,9,k,e) &
815  + dy(j,10) * u(i,10,k,e) &
816  + dy(j,11) * u(i,11,k,e) &
817  + dy(j,12) * u(i,12,k,e) &
818  + dy(j,13) * u(i,13,k,e)
819  end do
820  end do
821  end do
822 
823  do k = 1, lx
824  do i = 1, lx*lx
825  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
826  + dz(k,2) * u(i,1,2,e) &
827  + dz(k,3) * u(i,1,3,e) &
828  + dz(k,4) * u(i,1,4,e) &
829  + dz(k,5) * u(i,1,5,e) &
830  + dz(k,6) * u(i,1,6,e) &
831  + dz(k,7) * u(i,1,7,e) &
832  + dz(k,8) * u(i,1,8,e) &
833  + dz(k,9) * u(i,1,9,e) &
834  + dz(k,10) * u(i,1,10,e) &
835  + dz(k,11) * u(i,1,11,e) &
836  + dz(k,12) * u(i,1,12,e) &
837  + dz(k,13) * u(i,1,13,e)
838  end do
839  end do
840 
841  do i = 1, lx * lx * lx
842  ux(i,1,1,e) = w3(i,1,1) &
843  * ( drdx(i,1,1,e) * ur(i,1,1) &
844  + dsdx(i,1,1,e) * us(i,1,1) &
845  + dtdx(i,1,1,e) * ut(i,1,1) )
846  uy(i,1,1,e) = w3(i,1,1) &
847  * ( dsdy(i,1,1,e) * us(i,1,1) &
848  + drdy(i,1,1,e) * ur(i,1,1) &
849  + dtdy(i,1,1,e) * ut(i,1,1) )
850  uz(i,1,1,e) = w3(i,1,1) &
851  * ( dtdz(i,1,1,e) * ut(i,1,1) &
852  + drdz(i,1,1,e) * ur(i,1,1) &
853  + dsdz(i,1,1,e) * us(i,1,1) )
854  end do
855  end do
856  end subroutine cpu_opgrad_lx13
857 
858  subroutine cpu_opgrad_lx12(ux, uy, uz, u, dx, dy, dz, &
859  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
860  integer, parameter :: lx = 12
861  integer, intent(in) :: n
862  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
863  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
864  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
865  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
866  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
867  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
868  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
869  real(kind=rp) :: ur(lx, lx, lx)
870  real(kind=rp) :: us(lx, lx, lx)
871  real(kind=rp) :: ut(lx, lx, lx)
872  integer :: e, i, j, k
873 
874  do e = 1, n
875  do j = 1, lx * lx
876  do i = 1, lx
877  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
878  + dx(i,2) * u(2,j,1,e) &
879  + dx(i,3) * u(3,j,1,e) &
880  + dx(i,4) * u(4,j,1,e) &
881  + dx(i,5) * u(5,j,1,e) &
882  + dx(i,6) * u(6,j,1,e) &
883  + dx(i,7) * u(7,j,1,e) &
884  + dx(i,8) * u(8,j,1,e) &
885  + dx(i,9) * u(9,j,1,e) &
886  + dx(i,10) * u(10,j,1,e) &
887  + dx(i,11) * u(11,j,1,e) &
888  + dx(i,12) * u(12,j,1,e)
889  end do
890  end do
891 
892  do k = 1, lx
893  do j = 1, lx
894  do i = 1, lx
895  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
896  + dy(j,2) * u(i,2,k,e) &
897  + dy(j,3) * u(i,3,k,e) &
898  + dy(j,4) * u(i,4,k,e) &
899  + dy(j,5) * u(i,5,k,e) &
900  + dy(j,6) * u(i,6,k,e) &
901  + dy(j,7) * u(i,7,k,e) &
902  + dy(j,8) * u(i,8,k,e) &
903  + dy(j,9) * u(i,9,k,e) &
904  + dy(j,10) * u(i,10,k,e) &
905  + dy(j,11) * u(i,11,k,e) &
906  + dy(j,12) * u(i,12,k,e)
907  end do
908  end do
909  end do
910 
911  do k = 1, lx
912  do i = 1, lx*lx
913  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
914  + dz(k,2) * u(i,1,2,e) &
915  + dz(k,3) * u(i,1,3,e) &
916  + dz(k,4) * u(i,1,4,e) &
917  + dz(k,5) * u(i,1,5,e) &
918  + dz(k,6) * u(i,1,6,e) &
919  + dz(k,7) * u(i,1,7,e) &
920  + dz(k,8) * u(i,1,8,e) &
921  + dz(k,9) * u(i,1,9,e) &
922  + dz(k,10) * u(i,1,10,e) &
923  + dz(k,11) * u(i,1,11,e) &
924  + dz(k,12) * u(i,1,12,e)
925  end do
926  end do
927 
928  do i = 1, lx * lx * lx
929  ux(i,1,1,e) = w3(i,1,1) &
930  * ( drdx(i,1,1,e) * ur(i,1,1) &
931  + dsdx(i,1,1,e) * us(i,1,1) &
932  + dtdx(i,1,1,e) * ut(i,1,1) )
933  uy(i,1,1,e) = w3(i,1,1) &
934  * ( dsdy(i,1,1,e) * us(i,1,1) &
935  + drdy(i,1,1,e) * ur(i,1,1) &
936  + dtdy(i,1,1,e) * ut(i,1,1) )
937  uz(i,1,1,e) = w3(i,1,1) &
938  * ( dtdz(i,1,1,e) * ut(i,1,1) &
939  + drdz(i,1,1,e) * ur(i,1,1) &
940  + dsdz(i,1,1,e) * us(i,1,1) )
941  end do
942  end do
943  end subroutine cpu_opgrad_lx12
944 
945  subroutine cpu_opgrad_lx11(ux, uy, uz, u, dx, dy, dz, &
946  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
947  integer, parameter :: lx = 11
948  integer, intent(in) :: n
949  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
950  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
951  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
952  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
953  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
954  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
955  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
956  real(kind=rp) :: ur(lx, lx, lx)
957  real(kind=rp) :: us(lx, lx, lx)
958  real(kind=rp) :: ut(lx, lx, lx)
959  integer :: e, i, j, k
960 
961  do e = 1, n
962  do j = 1, lx * lx
963  do i = 1, lx
964  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
965  + dx(i,2) * u(2,j,1,e) &
966  + dx(i,3) * u(3,j,1,e) &
967  + dx(i,4) * u(4,j,1,e) &
968  + dx(i,5) * u(5,j,1,e) &
969  + dx(i,6) * u(6,j,1,e) &
970  + dx(i,7) * u(7,j,1,e) &
971  + dx(i,8) * u(8,j,1,e) &
972  + dx(i,9) * u(9,j,1,e) &
973  + dx(i,10) * u(10,j,1,e) &
974  + dx(i,11) * u(11,j,1,e)
975  end do
976  end do
977 
978  do k = 1, lx
979  do j = 1, lx
980  do i = 1, lx
981  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
982  + dy(j,2) * u(i,2,k,e) &
983  + dy(j,3) * u(i,3,k,e) &
984  + dy(j,4) * u(i,4,k,e) &
985  + dy(j,5) * u(i,5,k,e) &
986  + dy(j,6) * u(i,6,k,e) &
987  + dy(j,7) * u(i,7,k,e) &
988  + dy(j,8) * u(i,8,k,e) &
989  + dy(j,9) * u(i,9,k,e) &
990  + dy(j,10) * u(i,10,k,e) &
991  + dy(j,11) * u(i,11,k,e)
992  end do
993  end do
994  end do
995 
996  do k = 1, lx
997  do i = 1, lx*lx
998  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
999  + dz(k,2) * u(i,1,2,e) &
1000  + dz(k,3) * u(i,1,3,e) &
1001  + dz(k,4) * u(i,1,4,e) &
1002  + dz(k,5) * u(i,1,5,e) &
1003  + dz(k,6) * u(i,1,6,e) &
1004  + dz(k,7) * u(i,1,7,e) &
1005  + dz(k,8) * u(i,1,8,e) &
1006  + dz(k,9) * u(i,1,9,e) &
1007  + dz(k,10) * u(i,1,10,e) &
1008  + dz(k,11) * u(i,1,11,e)
1009  end do
1010  end do
1011 
1012  do i = 1, lx * lx * lx
1013  ux(i,1,1,e) = w3(i,1,1) &
1014  * ( drdx(i,1,1,e) * ur(i,1,1) &
1015  + dsdx(i,1,1,e) * us(i,1,1) &
1016  + dtdx(i,1,1,e) * ut(i,1,1) )
1017  uy(i,1,1,e) = w3(i,1,1) &
1018  * ( dsdy(i,1,1,e) * us(i,1,1) &
1019  + drdy(i,1,1,e) * ur(i,1,1) &
1020  + dtdy(i,1,1,e) * ut(i,1,1) )
1021  uz(i,1,1,e) = w3(i,1,1) &
1022  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1023  + drdz(i,1,1,e) * ur(i,1,1) &
1024  + dsdz(i,1,1,e) * us(i,1,1) )
1025  end do
1026  end do
1027  end subroutine cpu_opgrad_lx11
1028 
1029  subroutine cpu_opgrad_lx10(ux, uy, uz, u, dx, dy, dz, &
1030  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1031  integer, parameter :: lx = 10
1032  integer, intent(in) :: n
1033  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1034  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1035  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1036  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1037  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1038  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1039  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1040  real(kind=rp) :: ur(lx, lx, lx)
1041  real(kind=rp) :: us(lx, lx, lx)
1042  real(kind=rp) :: ut(lx, lx, lx)
1043  integer :: e, i, j, k
1044 
1045  do e = 1, n
1046  do j = 1, lx * lx
1047  do i = 1, lx
1048  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1049  + dx(i,2) * u(2,j,1,e) &
1050  + dx(i,3) * u(3,j,1,e) &
1051  + dx(i,4) * u(4,j,1,e) &
1052  + dx(i,5) * u(5,j,1,e) &
1053  + dx(i,6) * u(6,j,1,e) &
1054  + dx(i,7) * u(7,j,1,e) &
1055  + dx(i,8) * u(8,j,1,e) &
1056  + dx(i,9) * u(9,j,1,e) &
1057  + dx(i,10) * u(10,j,1,e)
1058  end do
1059  end do
1060 
1061  do k = 1, lx
1062  do j = 1, lx
1063  do i = 1, lx
1064  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1065  + dy(j,2) * u(i,2,k,e) &
1066  + dy(j,3) * u(i,3,k,e) &
1067  + dy(j,4) * u(i,4,k,e) &
1068  + dy(j,5) * u(i,5,k,e) &
1069  + dy(j,6) * u(i,6,k,e) &
1070  + dy(j,7) * u(i,7,k,e) &
1071  + dy(j,8) * u(i,8,k,e) &
1072  + dy(j,9) * u(i,9,k,e) &
1073  + dy(j,10) * u(i,10,k,e)
1074  end do
1075  end do
1076  end do
1077 
1078  do k = 1, lx
1079  do i = 1, lx*lx
1080  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1081  + dz(k,2) * u(i,1,2,e) &
1082  + dz(k,3) * u(i,1,3,e) &
1083  + dz(k,4) * u(i,1,4,e) &
1084  + dz(k,5) * u(i,1,5,e) &
1085  + dz(k,6) * u(i,1,6,e) &
1086  + dz(k,7) * u(i,1,7,e) &
1087  + dz(k,8) * u(i,1,8,e) &
1088  + dz(k,9) * u(i,1,9,e) &
1089  + dz(k,10) * u(i,1,10,e)
1090  end do
1091  end do
1092 
1093  do i = 1, lx * lx * lx
1094  ux(i,1,1,e) = w3(i,1,1) &
1095  * ( drdx(i,1,1,e) * ur(i,1,1) &
1096  + dsdx(i,1,1,e) * us(i,1,1) &
1097  + dtdx(i,1,1,e) * ut(i,1,1) )
1098  uy(i,1,1,e) = w3(i,1,1) &
1099  * ( dsdy(i,1,1,e) * us(i,1,1) &
1100  + drdy(i,1,1,e) * ur(i,1,1) &
1101  + dtdy(i,1,1,e) * ut(i,1,1) )
1102  uz(i,1,1,e) = w3(i,1,1) &
1103  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1104  + drdz(i,1,1,e) * ur(i,1,1) &
1105  + dsdz(i,1,1,e) * us(i,1,1) )
1106  end do
1107  end do
1108  end subroutine cpu_opgrad_lx10
1109 
1110  subroutine cpu_opgrad_lx9(ux, uy, uz, u, dx, dy, dz, &
1111  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1112  integer, parameter :: lx = 9
1113  integer, intent(in) :: n
1114  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1115  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1116  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1117  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1118  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1119  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1120  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1121  real(kind=rp) :: ur(lx, lx, lx)
1122  real(kind=rp) :: us(lx, lx, lx)
1123  real(kind=rp) :: ut(lx, lx, lx)
1124  integer :: e, i, j, k
1125 
1126  do e = 1, n
1127  do j = 1, lx * lx
1128  do i = 1, lx
1129  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1130  + dx(i,2) * u(2,j,1,e) &
1131  + dx(i,3) * u(3,j,1,e) &
1132  + dx(i,4) * u(4,j,1,e) &
1133  + dx(i,5) * u(5,j,1,e) &
1134  + dx(i,6) * u(6,j,1,e) &
1135  + dx(i,7) * u(7,j,1,e) &
1136  + dx(i,8) * u(8,j,1,e) &
1137  + dx(i,9) * u(9,j,1,e)
1138  end do
1139  end do
1140 
1141  do k = 1, lx
1142  do j = 1, lx
1143  do i = 1, lx
1144  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1145  + dy(j,2) * u(i,2,k,e) &
1146  + dy(j,3) * u(i,3,k,e) &
1147  + dy(j,4) * u(i,4,k,e) &
1148  + dy(j,5) * u(i,5,k,e) &
1149  + dy(j,6) * u(i,6,k,e) &
1150  + dy(j,7) * u(i,7,k,e) &
1151  + dy(j,8) * u(i,8,k,e) &
1152  + dy(j,9) * u(i,9,k,e)
1153  end do
1154  end do
1155  end do
1156 
1157  do k = 1, lx
1158  do i = 1, lx*lx
1159  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1160  + dz(k,2) * u(i,1,2,e) &
1161  + dz(k,3) * u(i,1,3,e) &
1162  + dz(k,4) * u(i,1,4,e) &
1163  + dz(k,5) * u(i,1,5,e) &
1164  + dz(k,6) * u(i,1,6,e) &
1165  + dz(k,7) * u(i,1,7,e) &
1166  + dz(k,8) * u(i,1,8,e) &
1167  + dz(k,9) * u(i,1,9,e)
1168  end do
1169  end do
1170 
1171  do i = 1, lx * lx * lx
1172  ux(i,1,1,e) = w3(i,1,1) &
1173  * ( drdx(i,1,1,e) * ur(i,1,1) &
1174  + dsdx(i,1,1,e) * us(i,1,1) &
1175  + dtdx(i,1,1,e) * ut(i,1,1) )
1176  uy(i,1,1,e) = w3(i,1,1) &
1177  * ( dsdy(i,1,1,e) * us(i,1,1) &
1178  + drdy(i,1,1,e) * ur(i,1,1) &
1179  + dtdy(i,1,1,e) * ut(i,1,1) )
1180  uz(i,1,1,e) = w3(i,1,1) &
1181  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1182  + drdz(i,1,1,e) * ur(i,1,1) &
1183  + dsdz(i,1,1,e) * us(i,1,1) )
1184  end do
1185  end do
1186  end subroutine cpu_opgrad_lx9
1187 
1188  subroutine cpu_opgrad_lx8(ux, uy, uz, u, dx, dy, dz, &
1189  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1190  integer, parameter :: lx = 8
1191  integer, intent(in) :: n
1192  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1193  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1194  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1195  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1196  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1197  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1198  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1199  real(kind=rp) :: ur(lx, lx, lx)
1200  real(kind=rp) :: us(lx, lx, lx)
1201  real(kind=rp) :: ut(lx, lx, lx)
1202  integer :: e, i, j, k
1203 
1204  do e = 1, n
1205  do j = 1, lx * lx
1206  do i = 1, lx
1207  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1208  + dx(i,2) * u(2,j,1,e) &
1209  + dx(i,3) * u(3,j,1,e) &
1210  + dx(i,4) * u(4,j,1,e) &
1211  + dx(i,5) * u(5,j,1,e) &
1212  + dx(i,6) * u(6,j,1,e) &
1213  + dx(i,7) * u(7,j,1,e) &
1214  + dx(i,8) * u(8,j,1,e)
1215  end do
1216  end do
1217 
1218  do k = 1, lx
1219  do j = 1, lx
1220  do i = 1, lx
1221  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1222  + dy(j,2) * u(i,2,k,e) &
1223  + dy(j,3) * u(i,3,k,e) &
1224  + dy(j,4) * u(i,4,k,e) &
1225  + dy(j,5) * u(i,5,k,e) &
1226  + dy(j,6) * u(i,6,k,e) &
1227  + dy(j,7) * u(i,7,k,e) &
1228  + dy(j,8) * u(i,8,k,e)
1229  end do
1230  end do
1231  end do
1232 
1233  do k = 1, lx
1234  do i = 1, lx*lx
1235  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1236  + dz(k,2) * u(i,1,2,e) &
1237  + dz(k,3) * u(i,1,3,e) &
1238  + dz(k,4) * u(i,1,4,e) &
1239  + dz(k,5) * u(i,1,5,e) &
1240  + dz(k,6) * u(i,1,6,e) &
1241  + dz(k,7) * u(i,1,7,e) &
1242  + dz(k,8) * u(i,1,8,e)
1243  end do
1244  end do
1245 
1246  do i = 1, lx * lx * lx
1247  ux(i,1,1,e) = w3(i,1,1) &
1248  * ( drdx(i,1,1,e) * ur(i,1,1) &
1249  + dsdx(i,1,1,e) * us(i,1,1) &
1250  + dtdx(i,1,1,e) * ut(i,1,1) )
1251  uy(i,1,1,e) = w3(i,1,1) &
1252  * ( dsdy(i,1,1,e) * us(i,1,1) &
1253  + drdy(i,1,1,e) * ur(i,1,1) &
1254  + dtdy(i,1,1,e) * ut(i,1,1) )
1255  uz(i,1,1,e) = w3(i,1,1) &
1256  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1257  + drdz(i,1,1,e) * ur(i,1,1) &
1258  + dsdz(i,1,1,e) * us(i,1,1) )
1259  end do
1260  end do
1261  end subroutine cpu_opgrad_lx8
1262 
1263  subroutine cpu_opgrad_lx7(ux, uy, uz, u, dx, dy, dz, &
1264  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1265  integer, parameter :: lx = 7
1266  integer, intent(in) :: n
1267  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1268  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1269  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1270  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1271  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1272  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1273  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1274  real(kind=rp) :: ur(lx, lx, lx)
1275  real(kind=rp) :: us(lx, lx, lx)
1276  real(kind=rp) :: ut(lx, lx, lx)
1277  integer :: e, i, j, k
1278 
1279  do e = 1, n
1280  do j = 1, lx * lx
1281  do i = 1, lx
1282  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1283  + dx(i,2) * u(2,j,1,e) &
1284  + dx(i,3) * u(3,j,1,e) &
1285  + dx(i,4) * u(4,j,1,e) &
1286  + dx(i,5) * u(5,j,1,e) &
1287  + dx(i,6) * u(6,j,1,e) &
1288  + dx(i,7) * u(7,j,1,e)
1289  end do
1290  end do
1291 
1292  do k = 1, lx
1293  do j = 1, lx
1294  do i = 1, lx
1295  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1296  + dy(j,2) * u(i,2,k,e) &
1297  + dy(j,3) * u(i,3,k,e) &
1298  + dy(j,4) * u(i,4,k,e) &
1299  + dy(j,5) * u(i,5,k,e) &
1300  + dy(j,6) * u(i,6,k,e) &
1301  + dy(j,7) * u(i,7,k,e)
1302  end do
1303  end do
1304  end do
1305 
1306  do k = 1, lx
1307  do i = 1, lx*lx
1308  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1309  + dz(k,2) * u(i,1,2,e) &
1310  + dz(k,3) * u(i,1,3,e) &
1311  + dz(k,4) * u(i,1,4,e) &
1312  + dz(k,5) * u(i,1,5,e) &
1313  + dz(k,6) * u(i,1,6,e) &
1314  + dz(k,7) * u(i,1,7,e)
1315  end do
1316  end do
1317 
1318  do i = 1, lx * lx * lx
1319  ux(i,1,1,e) = w3(i,1,1) &
1320  * ( drdx(i,1,1,e) * ur(i,1,1) &
1321  + dsdx(i,1,1,e) * us(i,1,1) &
1322  + dtdx(i,1,1,e) * ut(i,1,1) )
1323  uy(i,1,1,e) = w3(i,1,1) &
1324  * ( dsdy(i,1,1,e) * us(i,1,1) &
1325  + drdy(i,1,1,e) * ur(i,1,1) &
1326  + dtdy(i,1,1,e) * ut(i,1,1) )
1327  uz(i,1,1,e) = w3(i,1,1) &
1328  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1329  + drdz(i,1,1,e) * ur(i,1,1) &
1330  + dsdz(i,1,1,e) * us(i,1,1) )
1331  end do
1332  end do
1333  end subroutine cpu_opgrad_lx7
1334 
1335  subroutine cpu_opgrad_lx6(ux, uy, uz, u, dx, dy, dz, &
1336  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1337  integer, parameter :: lx = 6
1338  integer, intent(in) :: n
1339  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1340  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1341  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1342  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1343  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1344  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1345  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1346  real(kind=rp) :: ur(lx, lx, lx)
1347  real(kind=rp) :: us(lx, lx, lx)
1348  real(kind=rp) :: ut(lx, lx, lx)
1349  integer :: e, i, j, k
1350 
1351  do e = 1, n
1352  do j = 1, lx * lx
1353  do i = 1, lx
1354  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1355  + dx(i,2) * u(2,j,1,e) &
1356  + dx(i,3) * u(3,j,1,e) &
1357  + dx(i,4) * u(4,j,1,e) &
1358  + dx(i,5) * u(5,j,1,e) &
1359  + dx(i,6) * u(6,j,1,e)
1360  end do
1361  end do
1362 
1363  do k = 1, lx
1364  do j = 1, lx
1365  do i = 1, lx
1366  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1367  + dy(j,2) * u(i,2,k,e) &
1368  + dy(j,3) * u(i,3,k,e) &
1369  + dy(j,4) * u(i,4,k,e) &
1370  + dy(j,5) * u(i,5,k,e) &
1371  + dy(j,6) * u(i,6,k,e)
1372  end do
1373  end do
1374  end do
1375 
1376  do k = 1, lx
1377  do i = 1, lx*lx
1378  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1379  + dz(k,2) * u(i,1,2,e) &
1380  + dz(k,3) * u(i,1,3,e) &
1381  + dz(k,4) * u(i,1,4,e) &
1382  + dz(k,5) * u(i,1,5,e) &
1383  + dz(k,6) * u(i,1,6,e)
1384  end do
1385  end do
1386 
1387  do i = 1, lx * lx * lx
1388  ux(i,1,1,e) = w3(i,1,1) &
1389  * ( drdx(i,1,1,e) * ur(i,1,1) &
1390  + dsdx(i,1,1,e) * us(i,1,1) &
1391  + dtdx(i,1,1,e) * ut(i,1,1) )
1392  uy(i,1,1,e) = w3(i,1,1) &
1393  * ( dsdy(i,1,1,e) * us(i,1,1) &
1394  + drdy(i,1,1,e) * ur(i,1,1) &
1395  + dtdy(i,1,1,e) * ut(i,1,1) )
1396  uz(i,1,1,e) = w3(i,1,1) &
1397  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1398  + drdz(i,1,1,e) * ur(i,1,1) &
1399  + dsdz(i,1,1,e) * us(i,1,1) )
1400  end do
1401  end do
1402  end subroutine cpu_opgrad_lx6
1403 
1404  subroutine cpu_opgrad_lx5(ux, uy, uz, u, dx, dy, dz, &
1405  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1406  integer, parameter :: lx = 5
1407  integer, intent(in) :: n
1408  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1409  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1410  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1411  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1412  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1413  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1414  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1415  real(kind=rp) :: ur(lx, lx, lx)
1416  real(kind=rp) :: us(lx, lx, lx)
1417  real(kind=rp) :: ut(lx, lx, lx)
1418  integer :: e, i, j, k
1419 
1420  do e = 1, n
1421  do j = 1, lx * lx
1422  do i = 1, lx
1423  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1424  + dx(i,2) * u(2,j,1,e) &
1425  + dx(i,3) * u(3,j,1,e) &
1426  + dx(i,4) * u(4,j,1,e) &
1427  + dx(i,5) * u(5,j,1,e)
1428  end do
1429  end do
1430 
1431  do k = 1, lx
1432  do j = 1, lx
1433  do i = 1, lx
1434  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1435  + dy(j,2) * u(i,2,k,e) &
1436  + dy(j,3) * u(i,3,k,e) &
1437  + dy(j,4) * u(i,4,k,e) &
1438  + dy(j,5) * u(i,5,k,e)
1439  end do
1440  end do
1441  end do
1442 
1443  do k = 1, lx
1444  do i = 1, lx*lx
1445  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1446  + dz(k,2) * u(i,1,2,e) &
1447  + dz(k,3) * u(i,1,3,e) &
1448  + dz(k,4) * u(i,1,4,e) &
1449  + dz(k,5) * u(i,1,5,e)
1450  end do
1451  end do
1452 
1453  do i = 1, lx * lx * lx
1454  ux(i,1,1,e) = w3(i,1,1) &
1455  * ( drdx(i,1,1,e) * ur(i,1,1) &
1456  + dsdx(i,1,1,e) * us(i,1,1) &
1457  + dtdx(i,1,1,e) * ut(i,1,1) )
1458  uy(i,1,1,e) = w3(i,1,1) &
1459  * ( dsdy(i,1,1,e) * us(i,1,1) &
1460  + drdy(i,1,1,e) * ur(i,1,1) &
1461  + dtdy(i,1,1,e) * ut(i,1,1) )
1462  uz(i,1,1,e) = w3(i,1,1) &
1463  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1464  + drdz(i,1,1,e) * ur(i,1,1) &
1465  + dsdz(i,1,1,e) * us(i,1,1) )
1466  end do
1467  end do
1468  end subroutine cpu_opgrad_lx5
1469 
1470  subroutine cpu_opgrad_lx4(ux, uy, uz, u, dx, dy, dz, &
1471  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1472  integer, parameter :: lx = 4
1473  integer, intent(in) :: n
1474  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1475  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1476  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1477  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1478  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1479  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1480  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1481  real(kind=rp) :: ur(lx, lx, lx)
1482  real(kind=rp) :: us(lx, lx, lx)
1483  real(kind=rp) :: ut(lx, lx, lx)
1484  integer :: e, i, j, k
1485 
1486  do e = 1, n
1487  do j = 1, lx * lx
1488  do i = 1, lx
1489  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1490  + dx(i,2) * u(2,j,1,e) &
1491  + dx(i,3) * u(3,j,1,e) &
1492  + dx(i,4) * u(4,j,1,e)
1493  end do
1494  end do
1495 
1496  do k = 1, lx
1497  do j = 1, lx
1498  do i = 1, lx
1499  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1500  + dy(j,2) * u(i,2,k,e) &
1501  + dy(j,3) * u(i,3,k,e) &
1502  + dy(j,4) * u(i,4,k,e)
1503  end do
1504  end do
1505  end do
1506 
1507  do k = 1, lx
1508  do i = 1, lx*lx
1509  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1510  + dz(k,2) * u(i,1,2,e) &
1511  + dz(k,3) * u(i,1,3,e) &
1512  + dz(k,4) * u(i,1,4,e)
1513  end do
1514  end do
1515 
1516  do i = 1, lx * lx * lx
1517  ux(i,1,1,e) = w3(i,1,1) &
1518  * ( drdx(i,1,1,e) * ur(i,1,1) &
1519  + dsdx(i,1,1,e) * us(i,1,1) &
1520  + dtdx(i,1,1,e) * ut(i,1,1) )
1521  uy(i,1,1,e) = w3(i,1,1) &
1522  * ( dsdy(i,1,1,e) * us(i,1,1) &
1523  + drdy(i,1,1,e) * ur(i,1,1) &
1524  + dtdy(i,1,1,e) * ut(i,1,1) )
1525  uz(i,1,1,e) = w3(i,1,1) &
1526  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1527  + drdz(i,1,1,e) * ur(i,1,1) &
1528  + dsdz(i,1,1,e) * us(i,1,1) )
1529  end do
1530  end do
1531  end subroutine cpu_opgrad_lx4
1532 
1533  subroutine cpu_opgrad_lx3(ux, uy, uz, u, dx, dy, dz, &
1534  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1535  integer, parameter :: lx = 3
1536  integer, intent(in) :: n
1537  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1538  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1539  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1540  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1541  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1542  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1543  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1544  real(kind=rp) :: ur(lx, lx, lx)
1545  real(kind=rp) :: us(lx, lx, lx)
1546  real(kind=rp) :: ut(lx, lx, lx)
1547  integer :: e, i, j, k
1548 
1549  do e = 1, n
1550  do j = 1, lx * lx
1551  do i = 1, lx
1552  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1553  + dx(i,2) * u(2,j,1,e) &
1554  + dx(i,3) * u(3,j,1,e)
1555  end do
1556  end do
1557 
1558  do k = 1, lx
1559  do j = 1, lx
1560  do i = 1, lx
1561  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1562  + dy(j,2) * u(i,2,k,e) &
1563  + dy(j,3) * u(i,3,k,e)
1564  end do
1565  end do
1566  end do
1567 
1568  do k = 1, lx
1569  do i = 1, lx*lx
1570  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1571  + dz(k,2) * u(i,1,2,e) &
1572  + dz(k,3) * u(i,1,3,e)
1573  end do
1574  end do
1575 
1576  do i = 1, lx * lx * lx
1577  ux(i,1,1,e) = w3(i,1,1) &
1578  * ( drdx(i,1,1,e) * ur(i,1,1) &
1579  + dsdx(i,1,1,e) * us(i,1,1) &
1580  + dtdx(i,1,1,e) * ut(i,1,1) )
1581  uy(i,1,1,e) = w3(i,1,1) &
1582  * ( dsdy(i,1,1,e) * us(i,1,1) &
1583  + drdy(i,1,1,e) * ur(i,1,1) &
1584  + dtdy(i,1,1,e) * ut(i,1,1) )
1585  uz(i,1,1,e) = w3(i,1,1) &
1586  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1587  + drdz(i,1,1,e) * ur(i,1,1) &
1588  + dsdz(i,1,1,e) * us(i,1,1) )
1589  end do
1590  end do
1591  end subroutine cpu_opgrad_lx3
1592 
1593  subroutine cpu_opgrad_lx2(ux, uy, uz, u, dx, dy, dz, &
1594  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, n)
1595  integer, parameter :: lx = 2
1596  integer, intent(in) :: n
1597  real(kind=rp), dimension(lx, lx, lx, n), intent(inout) :: ux, uy, uz
1598  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: u
1599  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1600  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdx, dsdx, dtdx
1601  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdy, dsdy, dtdy
1602  real(kind=rp), dimension(lx, lx, lx, n), intent(in) :: drdz, dsdz, dtdz
1603  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1604  real(kind=rp) :: ur(lx, lx, lx)
1605  real(kind=rp) :: us(lx, lx, lx)
1606  real(kind=rp) :: ut(lx, lx, lx)
1607  integer :: e, i, j, k
1608 
1609  do e = 1, n
1610  do j = 1, lx * lx
1611  do i = 1, lx
1612  ur(i,j,1) = dx(i,1) * u(1,j,1,e) &
1613  + dx(i,2) * u(2,j,1,e)
1614  end do
1615  end do
1616 
1617  do k = 1, lx
1618  do j = 1, lx
1619  do i = 1, lx
1620  us(i,j,k) = dy(j,1) * u(i,1,k,e) &
1621  + dy(j,2) * u(i,2,k,e)
1622  end do
1623  end do
1624  end do
1625 
1626  do k = 1, lx
1627  do i = 1, lx*lx
1628  ut(i,1,k) = dz(k,1) * u(i,1,1,e) &
1629  + dz(k,2) * u(i,1,2,e)
1630  end do
1631  end do
1632 
1633  do i = 1, lx * lx * lx
1634  ux(i,1,1,e) = w3(i,1,1) &
1635  * ( drdx(i,1,1,e) * ur(i,1,1) &
1636  + dsdx(i,1,1,e) * us(i,1,1) &
1637  + dtdx(i,1,1,e) * ut(i,1,1) )
1638  uy(i,1,1,e) = w3(i,1,1) &
1639  * ( dsdy(i,1,1,e) * us(i,1,1) &
1640  + drdy(i,1,1,e) * ur(i,1,1) &
1641  + dtdy(i,1,1,e) * ut(i,1,1) )
1642  uz(i,1,1,e) = w3(i,1,1) &
1643  * ( dtdz(i,1,1,e) * ut(i,1,1) &
1644  + drdz(i,1,1,e) * ur(i,1,1) &
1645  + dsdz(i,1,1,e) * us(i,1,1) )
1646  end do
1647  end do
1648  end subroutine cpu_opgrad_lx2
1649 
1650  subroutine opr_cpu_opgrad_single(ux, uy, uz, u, coef, e)
1651  integer, parameter :: e_len = 1
1652  type(coef_t), intent(in) :: coef
1653  integer, intent(in) :: e
1654  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(inout) :: ux
1655  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(inout) :: uy
1656  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(inout) :: uz
1657  real(kind=rp), dimension(coef%Xh%lxyz, e_len), intent(in) :: u
1658 
1659  associate(xh => coef%Xh, msh => coef%msh, &
1660  drdx => coef%drdx, drdy => coef%drdy, drdz => coef%drdz, &
1661  dsdx => coef%dsdx, dsdy => coef%dsdy, dsdz => coef%dsdz, &
1662  dtdx => coef%dtdx, dtdy => coef%dtdy, dtdz => coef%dtdz)
1663 
1664  select case (xh%lx)
1665  case (18)
1666  call cpu_opgrad_lx18_single(ux, uy, uz, u, &
1667  xh%dx, xh%dy, xh%dz, &
1668  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1669  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1670  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1671  xh%w3)
1672  case (17)
1673  call cpu_opgrad_lx17_single(ux, uy, uz, u, &
1674  xh%dx, xh%dy, xh%dz, &
1675  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1676  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1677  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1678  xh%w3)
1679  case (16)
1680  call cpu_opgrad_lx16_single(ux, uy, uz, u, &
1681  xh%dx, xh%dy, xh%dz, &
1682  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1683  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1684  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1685  xh%w3)
1686  case (15)
1687  call cpu_opgrad_lx15_single(ux, uy, uz, u, &
1688  xh%dx, xh%dy, xh%dz, &
1689  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1690  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1691  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1692  xh%w3)
1693  case (14)
1694  call cpu_opgrad_lx14_single(ux, uy, uz, u, &
1695  xh%dx, xh%dy, xh%dz, &
1696  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1697  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1698  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1699  xh%w3)
1700  case (13)
1701  call cpu_opgrad_lx13_single(ux, uy, uz, u, &
1702  xh%dx, xh%dy, xh%dz, &
1703  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1704  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1705  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1706  xh%w3)
1707  case (12)
1708  call cpu_opgrad_lx12_single(ux, uy, uz, u, &
1709  xh%dx, xh%dy, xh%dz, &
1710  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1711  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1712  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1713  xh%w3)
1714  case (11)
1715  call cpu_opgrad_lx11_single(ux, uy, uz, u, &
1716  xh%dx, xh%dy, xh%dz, &
1717  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1718  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1719  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1720  xh%w3)
1721  case (10)
1722  call cpu_opgrad_lx10_single(ux, uy, uz, u, &
1723  xh%dx, xh%dy, xh%dz, &
1724  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1725  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1726  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1727  xh%w3)
1728 
1729  case (9)
1730  call cpu_opgrad_lx9_single(ux, uy, uz, u, &
1731  xh%dx, xh%dy, xh%dz, &
1732  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1733  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1734  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1735  xh%w3)
1736  case (8)
1737  call cpu_opgrad_lx8_single(ux, uy, uz, u, &
1738  xh%dx, xh%dy, xh%dz, &
1739  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1740  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1741  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1742  xh%w3)
1743  case (7)
1744  call cpu_opgrad_lx7_single(ux, uy, uz, u, &
1745  xh%dx, xh%dy, xh%dz, &
1746  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1747  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1748  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1749  xh%w3)
1750  case (6)
1751  call cpu_opgrad_lx6_single(ux, uy, uz, u, &
1752  xh%dx, xh%dy, xh%dz, &
1753  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1754  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1755  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1756  xh%w3)
1757  case (5)
1758  call cpu_opgrad_lx5_single(ux, uy, uz, u, &
1759  xh%dx, xh%dy, xh%dz, &
1760  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1761  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1762  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1763  xh%w3)
1764  case (4)
1765  call cpu_opgrad_lx4_single(ux, uy, uz, u, &
1766  xh%dx, xh%dy, xh%dz, &
1767  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1768  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1769  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1770  xh%w3)
1771  case (3)
1772  call cpu_opgrad_lx3_single(ux, uy, uz, u, &
1773  xh%dx, xh%dy, xh%dz, &
1774  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1775  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1776  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1777  xh%w3)
1778  case (2)
1779  call cpu_opgrad_lx2_single(ux, uy, uz, u, &
1780  xh%dx, xh%dy, xh%dz, &
1781  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1782  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1783  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1784  xh%w3)
1785  case default
1786  call cpu_opgrad_lx_single(ux, uy, uz, u, &
1787  xh%dx, xh%dy, xh%dz, &
1788  drdx(1,1,1,e), dsdx(1,1,1,e), dtdx(1,1,1,e), &
1789  drdy(1,1,1,e), dsdy(1,1,1,e), dtdy(1,1,1,e), &
1790  drdz(1,1,1,e), dsdz(1,1,1,e), dtdz(1,1,1,e), &
1791  xh%w3, xh%lx)
1792  end select
1793  end associate
1794 
1795  end subroutine opr_cpu_opgrad_single
1796 
1797  subroutine cpu_opgrad_lx_single(ux, uy, uz, u, dx, dy, dz, &
1798  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3, lx)
1799  integer, intent(in) :: lx
1800  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
1801  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
1802  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1803  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
1804  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
1805  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
1806  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1807  real(kind=rp) :: ur(lx, lx, lx)
1808  real(kind=rp) :: us(lx, lx, lx)
1809  real(kind=rp) :: ut(lx, lx, lx)
1810  real(kind=rp) :: tmp
1811  integer :: i, j, k, l
1812 
1813  do j = 1, lx * lx
1814  do i = 1, lx
1815  tmp = 0.0_rp
1816  do k = 1, lx
1817  tmp = tmp + dx(i,k) * u(k,j,1)
1818  end do
1819  ur(i,j,1) = tmp
1820  end do
1821  end do
1822 
1823  do k = 1, lx
1824  do j = 1, lx
1825  do i = 1, lx
1826  tmp = 0.0_rp
1827  do l = 1, lx
1828  tmp = tmp + dy(j,l) * u(i,l,k)
1829  end do
1830  us(i,j,k) = tmp
1831  end do
1832  end do
1833  end do
1834 
1835  do k = 1, lx
1836  do i = 1, lx*lx
1837  tmp = 0.0_rp
1838  do l = 1, lx
1839  tmp = tmp + dz(k,l) * u(i,1,l)
1840  end do
1841  ut(i,1,k) = tmp
1842  end do
1843  end do
1844 
1845  do i = 1, lx * lx * lx
1846  ux(i,1,1) = w3(i,1,1) &
1847  * ( drdx(i,1,1) * ur(i,1,1) &
1848  + dsdx(i,1,1) * us(i,1,1) &
1849  + dtdx(i,1,1) * ut(i,1,1) )
1850  uy(i,1,1) = w3(i,1,1) &
1851  * ( dsdy(i,1,1) * us(i,1,1) &
1852  + drdy(i,1,1) * ur(i,1,1) &
1853  + dtdy(i,1,1) * ut(i,1,1) )
1854  uz(i,1,1) = w3(i,1,1) &
1855  * ( dtdz(i,1,1) * ut(i,1,1) &
1856  + drdz(i,1,1) * ur(i,1,1) &
1857  + dsdz(i,1,1) * us(i,1,1) )
1858  end do
1859 
1860  end subroutine cpu_opgrad_lx_single
1861 
1862  subroutine cpu_opgrad_lx18_single(ux, uy, uz, u, dx, dy, dz, &
1863  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
1864  integer, parameter :: lx = 18
1865  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
1866  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
1867  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1868  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
1869  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
1870  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
1871  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1872  real(kind=rp) :: ur(lx, lx, lx)
1873  real(kind=rp) :: us(lx, lx, lx)
1874  real(kind=rp) :: ut(lx, lx, lx)
1875  integer :: i, j, k
1876 
1877  do j = 1, lx * lx
1878  do i = 1, lx
1879  ur(i,j,1) = dx(i,1) * u(1,j,1) &
1880  + dx(i,2) * u(2,j,1) &
1881  + dx(i,3) * u(3,j,1) &
1882  + dx(i,4) * u(4,j,1) &
1883  + dx(i,5) * u(5,j,1) &
1884  + dx(i,6) * u(6,j,1) &
1885  + dx(i,7) * u(7,j,1) &
1886  + dx(i,8) * u(8,j,1) &
1887  + dx(i,9) * u(9,j,1) &
1888  + dx(i,10) * u(10,j,1) &
1889  + dx(i,11) * u(11,j,1) &
1890  + dx(i,12) * u(12,j,1) &
1891  + dx(i,13) * u(13,j,1) &
1892  + dx(i,14) * u(14,j,1) &
1893  + dx(i,15) * u(15,j,1) &
1894  + dx(i,16) * u(16,j,1) &
1895  + dx(i,17) * u(17,j,1) &
1896  + dx(i,18) * u(18,j,1)
1897  end do
1898  end do
1899 
1900  do k = 1, lx
1901  do j = 1, lx
1902  do i = 1, lx
1903  us(i,j,k) = dy(j,1) * u(i,1,k) &
1904  + dy(j,2) * u(i,2,k) &
1905  + dy(j,3) * u(i,3,k) &
1906  + dy(j,4) * u(i,4,k) &
1907  + dy(j,5) * u(i,5,k) &
1908  + dy(j,6) * u(i,6,k) &
1909  + dy(j,7) * u(i,7,k) &
1910  + dy(j,8) * u(i,8,k) &
1911  + dy(j,9) * u(i,9,k) &
1912  + dy(j,10) * u(i,10,k) &
1913  + dy(j,11) * u(i,11,k) &
1914  + dy(j,12) * u(i,12,k) &
1915  + dy(j,13) * u(i,13,k) &
1916  + dy(j,14) * u(i,14,k) &
1917  + dy(j,15) * u(i,15,k) &
1918  + dy(j,16) * u(i,16,k) &
1919  + dy(j,17) * u(i,17,k) &
1920  + dy(j,18) * u(i,18,k)
1921  end do
1922  end do
1923  end do
1924 
1925  do k = 1, lx
1926  do i = 1, lx*lx
1927  ut(i,1,k) = dz(k,1) * u(i,1,1) &
1928  + dz(k,2) * u(i,1,2) &
1929  + dz(k,3) * u(i,1,3) &
1930  + dz(k,4) * u(i,1,4) &
1931  + dz(k,5) * u(i,1,5) &
1932  + dz(k,6) * u(i,1,6) &
1933  + dz(k,7) * u(i,1,7) &
1934  + dz(k,8) * u(i,1,8) &
1935  + dz(k,9) * u(i,1,9) &
1936  + dz(k,10) * u(i,1,10) &
1937  + dz(k,11) * u(i,1,11) &
1938  + dz(k,12) * u(i,1,12) &
1939  + dz(k,13) * u(i,1,13) &
1940  + dz(k,14) * u(i,1,14) &
1941  + dz(k,15) * u(i,1,15) &
1942  + dz(k,16) * u(i,1,16) &
1943  + dz(k,17) * u(i,1,17) &
1944  + dz(k,18) * u(i,1,18)
1945  end do
1946  end do
1947 
1948  do i = 1, lx * lx * lx
1949  ux(i,1,1) = w3(i,1,1) &
1950  * ( drdx(i,1,1) * ur(i,1,1) &
1951  + dsdx(i,1,1) * us(i,1,1) &
1952  + dtdx(i,1,1) * ut(i,1,1) )
1953  uy(i,1,1) = w3(i,1,1) &
1954  * ( dsdy(i,1,1) * us(i,1,1) &
1955  + drdy(i,1,1) * ur(i,1,1) &
1956  + dtdy(i,1,1) * ut(i,1,1) )
1957  uz(i,1,1) = w3(i,1,1) &
1958  * ( dtdz(i,1,1) * ut(i,1,1) &
1959  + drdz(i,1,1) * ur(i,1,1) &
1960  + dsdz(i,1,1) * us(i,1,1) )
1961  end do
1962 
1963  end subroutine cpu_opgrad_lx18_single
1964 
1965  subroutine cpu_opgrad_lx17_single(ux, uy, uz, u, dx, dy, dz, &
1966  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
1967  integer, parameter :: lx = 17
1968  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
1969  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
1970  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
1971  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
1972  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
1973  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
1974  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
1975  real(kind=rp) :: ur(lx, lx, lx)
1976  real(kind=rp) :: us(lx, lx, lx)
1977  real(kind=rp) :: ut(lx, lx, lx)
1978  integer :: i, j, k
1979 
1980  do j = 1, lx * lx
1981  do i = 1, lx
1982  ur(i,j,1) = dx(i,1) * u(1,j,1) &
1983  + dx(i,2) * u(2,j,1) &
1984  + dx(i,3) * u(3,j,1) &
1985  + dx(i,4) * u(4,j,1) &
1986  + dx(i,5) * u(5,j,1) &
1987  + dx(i,6) * u(6,j,1) &
1988  + dx(i,7) * u(7,j,1) &
1989  + dx(i,8) * u(8,j,1) &
1990  + dx(i,9) * u(9,j,1) &
1991  + dx(i,10) * u(10,j,1) &
1992  + dx(i,11) * u(11,j,1) &
1993  + dx(i,12) * u(12,j,1) &
1994  + dx(i,13) * u(13,j,1) &
1995  + dx(i,14) * u(14,j,1) &
1996  + dx(i,15) * u(15,j,1) &
1997  + dx(i,16) * u(16,j,1) &
1998  + dx(i,17) * u(17,j,1)
1999  end do
2000  end do
2001 
2002  do k = 1, lx
2003  do j = 1, lx
2004  do i = 1, lx
2005  us(i,j,k) = dy(j,1) * u(i,1,k) &
2006  + dy(j,2) * u(i,2,k) &
2007  + dy(j,3) * u(i,3,k) &
2008  + dy(j,4) * u(i,4,k) &
2009  + dy(j,5) * u(i,5,k) &
2010  + dy(j,6) * u(i,6,k) &
2011  + dy(j,7) * u(i,7,k) &
2012  + dy(j,8) * u(i,8,k) &
2013  + dy(j,9) * u(i,9,k) &
2014  + dy(j,10) * u(i,10,k) &
2015  + dy(j,11) * u(i,11,k) &
2016  + dy(j,12) * u(i,12,k) &
2017  + dy(j,13) * u(i,13,k) &
2018  + dy(j,14) * u(i,14,k) &
2019  + dy(j,15) * u(i,15,k) &
2020  + dy(j,16) * u(i,16,k) &
2021  + dy(j,17) * u(i,17,k)
2022  end do
2023  end do
2024  end do
2025 
2026  do k = 1, lx
2027  do i = 1, lx*lx
2028  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2029  + dz(k,2) * u(i,1,2) &
2030  + dz(k,3) * u(i,1,3) &
2031  + dz(k,4) * u(i,1,4) &
2032  + dz(k,5) * u(i,1,5) &
2033  + dz(k,6) * u(i,1,6) &
2034  + dz(k,7) * u(i,1,7) &
2035  + dz(k,8) * u(i,1,8) &
2036  + dz(k,9) * u(i,1,9) &
2037  + dz(k,10) * u(i,1,10) &
2038  + dz(k,11) * u(i,1,11) &
2039  + dz(k,12) * u(i,1,12) &
2040  + dz(k,13) * u(i,1,13) &
2041  + dz(k,14) * u(i,1,14) &
2042  + dz(k,15) * u(i,1,15) &
2043  + dz(k,16) * u(i,1,16) &
2044  + dz(k,17) * u(i,1,17)
2045  end do
2046  end do
2047 
2048  do i = 1, lx * lx * lx
2049  ux(i,1,1) = w3(i,1,1) &
2050  * ( drdx(i,1,1) * ur(i,1,1) &
2051  + dsdx(i,1,1) * us(i,1,1) &
2052  + dtdx(i,1,1) * ut(i,1,1) )
2053  uy(i,1,1) = w3(i,1,1) &
2054  * ( dsdy(i,1,1) * us(i,1,1) &
2055  + drdy(i,1,1) * ur(i,1,1) &
2056  + dtdy(i,1,1) * ut(i,1,1) )
2057  uz(i,1,1) = w3(i,1,1) &
2058  * ( dtdz(i,1,1) * ut(i,1,1) &
2059  + drdz(i,1,1) * ur(i,1,1) &
2060  + dsdz(i,1,1) * us(i,1,1) )
2061  end do
2062  end subroutine cpu_opgrad_lx17_single
2063 
2064  subroutine cpu_opgrad_lx16_single(ux, uy, uz, u, dx, dy, dz, &
2065  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2066  integer, parameter :: lx = 16
2067  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2068  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2069  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2070  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2071  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2072  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2073  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2074  real(kind=rp) :: ur(lx, lx, lx)
2075  real(kind=rp) :: us(lx, lx, lx)
2076  real(kind=rp) :: ut(lx, lx, lx)
2077  integer :: i, j, k
2078 
2079  do j = 1, lx * lx
2080  do i = 1, lx
2081  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2082  + dx(i,2) * u(2,j,1) &
2083  + dx(i,3) * u(3,j,1) &
2084  + dx(i,4) * u(4,j,1) &
2085  + dx(i,5) * u(5,j,1) &
2086  + dx(i,6) * u(6,j,1) &
2087  + dx(i,7) * u(7,j,1) &
2088  + dx(i,8) * u(8,j,1) &
2089  + dx(i,9) * u(9,j,1) &
2090  + dx(i,10) * u(10,j,1) &
2091  + dx(i,11) * u(11,j,1) &
2092  + dx(i,12) * u(12,j,1) &
2093  + dx(i,13) * u(13,j,1) &
2094  + dx(i,14) * u(14,j,1) &
2095  + dx(i,15) * u(15,j,1) &
2096  + dx(i,16) * u(16,j,1)
2097  end do
2098  end do
2099 
2100  do k = 1, lx
2101  do j = 1, lx
2102  do i = 1, lx
2103  us(i,j,k) = dy(j,1) * u(i,1,k) &
2104  + dy(j,2) * u(i,2,k) &
2105  + dy(j,3) * u(i,3,k) &
2106  + dy(j,4) * u(i,4,k) &
2107  + dy(j,5) * u(i,5,k) &
2108  + dy(j,6) * u(i,6,k) &
2109  + dy(j,7) * u(i,7,k) &
2110  + dy(j,8) * u(i,8,k) &
2111  + dy(j,9) * u(i,9,k) &
2112  + dy(j,10) * u(i,10,k) &
2113  + dy(j,11) * u(i,11,k) &
2114  + dy(j,12) * u(i,12,k) &
2115  + dy(j,13) * u(i,13,k) &
2116  + dy(j,14) * u(i,14,k) &
2117  + dy(j,15) * u(i,15,k) &
2118  + dy(j,16) * u(i,16,k)
2119  end do
2120  end do
2121  end do
2122 
2123  do k = 1, lx
2124  do i = 1, lx*lx
2125  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2126  + dz(k,2) * u(i,1,2) &
2127  + dz(k,3) * u(i,1,3) &
2128  + dz(k,4) * u(i,1,4) &
2129  + dz(k,5) * u(i,1,5) &
2130  + dz(k,6) * u(i,1,6) &
2131  + dz(k,7) * u(i,1,7) &
2132  + dz(k,8) * u(i,1,8) &
2133  + dz(k,9) * u(i,1,9) &
2134  + dz(k,10) * u(i,1,10) &
2135  + dz(k,11) * u(i,1,11) &
2136  + dz(k,12) * u(i,1,12) &
2137  + dz(k,13) * u(i,1,13) &
2138  + dz(k,14) * u(i,1,14) &
2139  + dz(k,15) * u(i,1,15) &
2140  + dz(k,16) * u(i,1,16)
2141  end do
2142  end do
2143 
2144  do i = 1, lx * lx * lx
2145  ux(i,1,1) = w3(i,1,1) &
2146  * ( drdx(i,1,1) * ur(i,1,1) &
2147  + dsdx(i,1,1) * us(i,1,1) &
2148  + dtdx(i,1,1) * ut(i,1,1) )
2149  uy(i,1,1) = w3(i,1,1) &
2150  * ( dsdy(i,1,1) * us(i,1,1) &
2151  + drdy(i,1,1) * ur(i,1,1) &
2152  + dtdy(i,1,1) * ut(i,1,1) )
2153  uz(i,1,1) = w3(i,1,1) &
2154  * ( dtdz(i,1,1) * ut(i,1,1) &
2155  + drdz(i,1,1) * ur(i,1,1) &
2156  + dsdz(i,1,1) * us(i,1,1) )
2157  end do
2158  end subroutine cpu_opgrad_lx16_single
2159 
2160  subroutine cpu_opgrad_lx15_single(ux, uy, uz, u, dx, dy, dz, &
2161  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2162  integer, parameter :: lx = 15
2163  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2164  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2165  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2166  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2167  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2168  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2169  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2170  real(kind=rp) :: ur(lx, lx, lx)
2171  real(kind=rp) :: us(lx, lx, lx)
2172  real(kind=rp) :: ut(lx, lx, lx)
2173  integer :: i, j, k
2174 
2175  do j = 1, lx * lx
2176  do i = 1, lx
2177  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2178  + dx(i,2) * u(2,j,1) &
2179  + dx(i,3) * u(3,j,1) &
2180  + dx(i,4) * u(4,j,1) &
2181  + dx(i,5) * u(5,j,1) &
2182  + dx(i,6) * u(6,j,1) &
2183  + dx(i,7) * u(7,j,1) &
2184  + dx(i,8) * u(8,j,1) &
2185  + dx(i,9) * u(9,j,1) &
2186  + dx(i,10) * u(10,j,1) &
2187  + dx(i,11) * u(11,j,1) &
2188  + dx(i,12) * u(12,j,1) &
2189  + dx(i,13) * u(13,j,1) &
2190  + dx(i,14) * u(14,j,1) &
2191  + dx(i,15) * u(15,j,1)
2192  end do
2193  end do
2194 
2195  do k = 1, lx
2196  do j = 1, lx
2197  do i = 1, lx
2198  us(i,j,k) = dy(j,1) * u(i,1,k) &
2199  + dy(j,2) * u(i,2,k) &
2200  + dy(j,3) * u(i,3,k) &
2201  + dy(j,4) * u(i,4,k) &
2202  + dy(j,5) * u(i,5,k) &
2203  + dy(j,6) * u(i,6,k) &
2204  + dy(j,7) * u(i,7,k) &
2205  + dy(j,8) * u(i,8,k) &
2206  + dy(j,9) * u(i,9,k) &
2207  + dy(j,10) * u(i,10,k) &
2208  + dy(j,11) * u(i,11,k) &
2209  + dy(j,12) * u(i,12,k) &
2210  + dy(j,13) * u(i,13,k) &
2211  + dy(j,14) * u(i,14,k) &
2212  + dy(j,15) * u(i,15,k)
2213  end do
2214  end do
2215  end do
2216 
2217  do k = 1, lx
2218  do i = 1, lx*lx
2219  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2220  + dz(k,2) * u(i,1,2) &
2221  + dz(k,3) * u(i,1,3) &
2222  + dz(k,4) * u(i,1,4) &
2223  + dz(k,5) * u(i,1,5) &
2224  + dz(k,6) * u(i,1,6) &
2225  + dz(k,7) * u(i,1,7) &
2226  + dz(k,8) * u(i,1,8) &
2227  + dz(k,9) * u(i,1,9) &
2228  + dz(k,10) * u(i,1,10) &
2229  + dz(k,11) * u(i,1,11) &
2230  + dz(k,12) * u(i,1,12) &
2231  + dz(k,13) * u(i,1,13) &
2232  + dz(k,14) * u(i,1,14) &
2233  + dz(k,15) * u(i,1,15)
2234  end do
2235  end do
2236 
2237  do i = 1, lx * lx * lx
2238  ux(i,1,1) = w3(i,1,1) &
2239  * ( drdx(i,1,1) * ur(i,1,1) &
2240  + dsdx(i,1,1) * us(i,1,1) &
2241  + dtdx(i,1,1) * ut(i,1,1) )
2242  uy(i,1,1) = w3(i,1,1) &
2243  * ( dsdy(i,1,1) * us(i,1,1) &
2244  + drdy(i,1,1) * ur(i,1,1) &
2245  + dtdy(i,1,1) * ut(i,1,1) )
2246  uz(i,1,1) = w3(i,1,1) &
2247  * ( dtdz(i,1,1) * ut(i,1,1) &
2248  + drdz(i,1,1) * ur(i,1,1) &
2249  + dsdz(i,1,1) * us(i,1,1) )
2250  end do
2251  end subroutine cpu_opgrad_lx15_single
2252 
2253  subroutine cpu_opgrad_lx14_single(ux, uy, uz, u, dx, dy, dz, &
2254  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2255  integer, parameter :: lx = 14
2256  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2257  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2258  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2259  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2260  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2261  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2262  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2263  real(kind=rp) :: ur(lx, lx, lx)
2264  real(kind=rp) :: us(lx, lx, lx)
2265  real(kind=rp) :: ut(lx, lx, lx)
2266  integer :: i, j, k
2267 
2268  do j = 1, lx * lx
2269  do i = 1, lx
2270  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2271  + dx(i,2) * u(2,j,1) &
2272  + dx(i,3) * u(3,j,1) &
2273  + dx(i,4) * u(4,j,1) &
2274  + dx(i,5) * u(5,j,1) &
2275  + dx(i,6) * u(6,j,1) &
2276  + dx(i,7) * u(7,j,1) &
2277  + dx(i,8) * u(8,j,1) &
2278  + dx(i,9) * u(9,j,1) &
2279  + dx(i,10) * u(10,j,1) &
2280  + dx(i,11) * u(11,j,1) &
2281  + dx(i,12) * u(12,j,1) &
2282  + dx(i,13) * u(13,j,1) &
2283  + dx(i,14) * u(14,j,1)
2284  end do
2285  end do
2286 
2287  do k = 1, lx
2288  do j = 1, lx
2289  do i = 1, lx
2290  us(i,j,k) = dy(j,1) * u(i,1,k) &
2291  + dy(j,2) * u(i,2,k) &
2292  + dy(j,3) * u(i,3,k) &
2293  + dy(j,4) * u(i,4,k) &
2294  + dy(j,5) * u(i,5,k) &
2295  + dy(j,6) * u(i,6,k) &
2296  + dy(j,7) * u(i,7,k) &
2297  + dy(j,8) * u(i,8,k) &
2298  + dy(j,9) * u(i,9,k) &
2299  + dy(j,10) * u(i,10,k) &
2300  + dy(j,11) * u(i,11,k) &
2301  + dy(j,12) * u(i,12,k) &
2302  + dy(j,13) * u(i,13,k) &
2303  + dy(j,14) * u(i,14,k)
2304  end do
2305  end do
2306  end do
2307 
2308  do k = 1, lx
2309  do i = 1, lx*lx
2310  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2311  + dz(k,2) * u(i,1,2) &
2312  + dz(k,3) * u(i,1,3) &
2313  + dz(k,4) * u(i,1,4) &
2314  + dz(k,5) * u(i,1,5) &
2315  + dz(k,6) * u(i,1,6) &
2316  + dz(k,7) * u(i,1,7) &
2317  + dz(k,8) * u(i,1,8) &
2318  + dz(k,9) * u(i,1,9) &
2319  + dz(k,10) * u(i,1,10) &
2320  + dz(k,11) * u(i,1,11) &
2321  + dz(k,12) * u(i,1,12) &
2322  + dz(k,13) * u(i,1,13) &
2323  + dz(k,14) * u(i,1,14)
2324  end do
2325  end do
2326 
2327  do i = 1, lx * lx * lx
2328  ux(i,1,1) = w3(i,1,1) &
2329  * ( drdx(i,1,1) * ur(i,1,1) &
2330  + dsdx(i,1,1) * us(i,1,1) &
2331  + dtdx(i,1,1) * ut(i,1,1) )
2332  uy(i,1,1) = w3(i,1,1) &
2333  * ( dsdy(i,1,1) * us(i,1,1) &
2334  + drdy(i,1,1) * ur(i,1,1) &
2335  + dtdy(i,1,1) * ut(i,1,1) )
2336  uz(i,1,1) = w3(i,1,1) &
2337  * ( dtdz(i,1,1) * ut(i,1,1) &
2338  + drdz(i,1,1) * ur(i,1,1) &
2339  + dsdz(i,1,1) * us(i,1,1) )
2340  end do
2341  end subroutine cpu_opgrad_lx14_single
2342 
2343  subroutine cpu_opgrad_lx13_single(ux, uy, uz, u, dx, dy, dz, &
2344  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2345  integer, parameter :: lx = 13
2346  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2347  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2348  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2349  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2350  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2351  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2352  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2353  real(kind=rp) :: ur(lx, lx, lx)
2354  real(kind=rp) :: us(lx, lx, lx)
2355  real(kind=rp) :: ut(lx, lx, lx)
2356  integer :: i, j, k
2357 
2358  do j = 1, lx * lx
2359  do i = 1, lx
2360  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2361  + dx(i,2) * u(2,j,1) &
2362  + dx(i,3) * u(3,j,1) &
2363  + dx(i,4) * u(4,j,1) &
2364  + dx(i,5) * u(5,j,1) &
2365  + dx(i,6) * u(6,j,1) &
2366  + dx(i,7) * u(7,j,1) &
2367  + dx(i,8) * u(8,j,1) &
2368  + dx(i,9) * u(9,j,1) &
2369  + dx(i,10) * u(10,j,1) &
2370  + dx(i,11) * u(11,j,1) &
2371  + dx(i,12) * u(12,j,1) &
2372  + dx(i,13) * u(13,j,1)
2373  end do
2374  end do
2375 
2376  do k = 1, lx
2377  do j = 1, lx
2378  do i = 1, lx
2379  us(i,j,k) = dy(j,1) * u(i,1,k) &
2380  + dy(j,2) * u(i,2,k) &
2381  + dy(j,3) * u(i,3,k) &
2382  + dy(j,4) * u(i,4,k) &
2383  + dy(j,5) * u(i,5,k) &
2384  + dy(j,6) * u(i,6,k) &
2385  + dy(j,7) * u(i,7,k) &
2386  + dy(j,8) * u(i,8,k) &
2387  + dy(j,9) * u(i,9,k) &
2388  + dy(j,10) * u(i,10,k) &
2389  + dy(j,11) * u(i,11,k) &
2390  + dy(j,12) * u(i,12,k) &
2391  + dy(j,13) * u(i,13,k)
2392  end do
2393  end do
2394  end do
2395 
2396  do k = 1, lx
2397  do i = 1, lx*lx
2398  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2399  + dz(k,2) * u(i,1,2) &
2400  + dz(k,3) * u(i,1,3) &
2401  + dz(k,4) * u(i,1,4) &
2402  + dz(k,5) * u(i,1,5) &
2403  + dz(k,6) * u(i,1,6) &
2404  + dz(k,7) * u(i,1,7) &
2405  + dz(k,8) * u(i,1,8) &
2406  + dz(k,9) * u(i,1,9) &
2407  + dz(k,10) * u(i,1,10) &
2408  + dz(k,11) * u(i,1,11) &
2409  + dz(k,12) * u(i,1,12) &
2410  + dz(k,13) * u(i,1,13)
2411  end do
2412  end do
2413 
2414  do i = 1, lx * lx * lx
2415  ux(i,1,1) = w3(i,1,1) &
2416  * ( drdx(i,1,1) * ur(i,1,1) &
2417  + dsdx(i,1,1) * us(i,1,1) &
2418  + dtdx(i,1,1) * ut(i,1,1) )
2419  uy(i,1,1) = w3(i,1,1) &
2420  * ( dsdy(i,1,1) * us(i,1,1) &
2421  + drdy(i,1,1) * ur(i,1,1) &
2422  + dtdy(i,1,1) * ut(i,1,1) )
2423  uz(i,1,1) = w3(i,1,1) &
2424  * ( dtdz(i,1,1) * ut(i,1,1) &
2425  + drdz(i,1,1) * ur(i,1,1) &
2426  + dsdz(i,1,1) * us(i,1,1) )
2427  end do
2428  end subroutine cpu_opgrad_lx13_single
2429 
2430  subroutine cpu_opgrad_lx12_single(ux, uy, uz, u, dx, dy, dz, &
2431  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2432  integer, parameter :: lx = 12
2433  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2434  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2435  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2436  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2437  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2438  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2439  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2440  real(kind=rp) :: ur(lx, lx, lx)
2441  real(kind=rp) :: us(lx, lx, lx)
2442  real(kind=rp) :: ut(lx, lx, lx)
2443  integer :: i, j, k
2444 
2445  do j = 1, lx * lx
2446  do i = 1, lx
2447  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2448  + dx(i,2) * u(2,j,1) &
2449  + dx(i,3) * u(3,j,1) &
2450  + dx(i,4) * u(4,j,1) &
2451  + dx(i,5) * u(5,j,1) &
2452  + dx(i,6) * u(6,j,1) &
2453  + dx(i,7) * u(7,j,1) &
2454  + dx(i,8) * u(8,j,1) &
2455  + dx(i,9) * u(9,j,1) &
2456  + dx(i,10) * u(10,j,1) &
2457  + dx(i,11) * u(11,j,1) &
2458  + dx(i,12) * u(12,j,1)
2459  end do
2460  end do
2461 
2462  do k = 1, lx
2463  do j = 1, lx
2464  do i = 1, lx
2465  us(i,j,k) = dy(j,1) * u(i,1,k) &
2466  + dy(j,2) * u(i,2,k) &
2467  + dy(j,3) * u(i,3,k) &
2468  + dy(j,4) * u(i,4,k) &
2469  + dy(j,5) * u(i,5,k) &
2470  + dy(j,6) * u(i,6,k) &
2471  + dy(j,7) * u(i,7,k) &
2472  + dy(j,8) * u(i,8,k) &
2473  + dy(j,9) * u(i,9,k) &
2474  + dy(j,10) * u(i,10,k) &
2475  + dy(j,11) * u(i,11,k) &
2476  + dy(j,12) * u(i,12,k)
2477  end do
2478  end do
2479  end do
2480 
2481  do k = 1, lx
2482  do i = 1, lx*lx
2483  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2484  + dz(k,2) * u(i,1,2) &
2485  + dz(k,3) * u(i,1,3) &
2486  + dz(k,4) * u(i,1,4) &
2487  + dz(k,5) * u(i,1,5) &
2488  + dz(k,6) * u(i,1,6) &
2489  + dz(k,7) * u(i,1,7) &
2490  + dz(k,8) * u(i,1,8) &
2491  + dz(k,9) * u(i,1,9) &
2492  + dz(k,10) * u(i,1,10) &
2493  + dz(k,11) * u(i,1,11) &
2494  + dz(k,12) * u(i,1,12)
2495  end do
2496  end do
2497 
2498  do i = 1, lx * lx * lx
2499  ux(i,1,1) = w3(i,1,1) &
2500  * ( drdx(i,1,1) * ur(i,1,1) &
2501  + dsdx(i,1,1) * us(i,1,1) &
2502  + dtdx(i,1,1) * ut(i,1,1) )
2503  uy(i,1,1) = w3(i,1,1) &
2504  * ( dsdy(i,1,1) * us(i,1,1) &
2505  + drdy(i,1,1) * ur(i,1,1) &
2506  + dtdy(i,1,1) * ut(i,1,1) )
2507  uz(i,1,1) = w3(i,1,1) &
2508  * ( dtdz(i,1,1) * ut(i,1,1) &
2509  + drdz(i,1,1) * ur(i,1,1) &
2510  + dsdz(i,1,1) * us(i,1,1) )
2511  end do
2512  end subroutine cpu_opgrad_lx12_single
2513 
2514  subroutine cpu_opgrad_lx11_single(ux, uy, uz, u, dx, dy, dz, &
2515  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2516  integer, parameter :: lx = 11
2517  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2518  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2519  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2520  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2521  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2522  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2523  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2524  real(kind=rp) :: ur(lx, lx, lx)
2525  real(kind=rp) :: us(lx, lx, lx)
2526  real(kind=rp) :: ut(lx, lx, lx)
2527  integer :: i, j, k
2528 
2529  do j = 1, lx * lx
2530  do i = 1, lx
2531  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2532  + dx(i,2) * u(2,j,1) &
2533  + dx(i,3) * u(3,j,1) &
2534  + dx(i,4) * u(4,j,1) &
2535  + dx(i,5) * u(5,j,1) &
2536  + dx(i,6) * u(6,j,1) &
2537  + dx(i,7) * u(7,j,1) &
2538  + dx(i,8) * u(8,j,1) &
2539  + dx(i,9) * u(9,j,1) &
2540  + dx(i,10) * u(10,j,1) &
2541  + dx(i,11) * u(11,j,1)
2542  end do
2543  end do
2544 
2545  do k = 1, lx
2546  do j = 1, lx
2547  do i = 1, lx
2548  us(i,j,k) = dy(j,1) * u(i,1,k) &
2549  + dy(j,2) * u(i,2,k) &
2550  + dy(j,3) * u(i,3,k) &
2551  + dy(j,4) * u(i,4,k) &
2552  + dy(j,5) * u(i,5,k) &
2553  + dy(j,6) * u(i,6,k) &
2554  + dy(j,7) * u(i,7,k) &
2555  + dy(j,8) * u(i,8,k) &
2556  + dy(j,9) * u(i,9,k) &
2557  + dy(j,10) * u(i,10,k) &
2558  + dy(j,11) * u(i,11,k)
2559  end do
2560  end do
2561  end do
2562 
2563  do k = 1, lx
2564  do i = 1, lx*lx
2565  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2566  + dz(k,2) * u(i,1,2) &
2567  + dz(k,3) * u(i,1,3) &
2568  + dz(k,4) * u(i,1,4) &
2569  + dz(k,5) * u(i,1,5) &
2570  + dz(k,6) * u(i,1,6) &
2571  + dz(k,7) * u(i,1,7) &
2572  + dz(k,8) * u(i,1,8) &
2573  + dz(k,9) * u(i,1,9) &
2574  + dz(k,10) * u(i,1,10) &
2575  + dz(k,11) * u(i,1,11)
2576  end do
2577  end do
2578 
2579  do i = 1, lx * lx * lx
2580  ux(i,1,1) = w3(i,1,1) &
2581  * ( drdx(i,1,1) * ur(i,1,1) &
2582  + dsdx(i,1,1) * us(i,1,1) &
2583  + dtdx(i,1,1) * ut(i,1,1) )
2584  uy(i,1,1) = w3(i,1,1) &
2585  * ( dsdy(i,1,1) * us(i,1,1) &
2586  + drdy(i,1,1) * ur(i,1,1) &
2587  + dtdy(i,1,1) * ut(i,1,1) )
2588  uz(i,1,1) = w3(i,1,1) &
2589  * ( dtdz(i,1,1) * ut(i,1,1) &
2590  + drdz(i,1,1) * ur(i,1,1) &
2591  + dsdz(i,1,1) * us(i,1,1) )
2592  end do
2593  end subroutine cpu_opgrad_lx11_single
2594 
2595  subroutine cpu_opgrad_lx10_single(ux, uy, uz, u, dx, dy, dz, &
2596  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2597  integer, parameter :: lx = 10
2598  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2599  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2600  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2601  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2602  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2603  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2604  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2605  real(kind=rp) :: ur(lx, lx, lx)
2606  real(kind=rp) :: us(lx, lx, lx)
2607  real(kind=rp) :: ut(lx, lx, lx)
2608  integer :: i, j, k
2609 
2610  do j = 1, lx * lx
2611  do i = 1, lx
2612  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2613  + dx(i,2) * u(2,j,1) &
2614  + dx(i,3) * u(3,j,1) &
2615  + dx(i,4) * u(4,j,1) &
2616  + dx(i,5) * u(5,j,1) &
2617  + dx(i,6) * u(6,j,1) &
2618  + dx(i,7) * u(7,j,1) &
2619  + dx(i,8) * u(8,j,1) &
2620  + dx(i,9) * u(9,j,1) &
2621  + dx(i,10) * u(10,j,1)
2622  end do
2623  end do
2624 
2625  do k = 1, lx
2626  do j = 1, lx
2627  do i = 1, lx
2628  us(i,j,k) = dy(j,1) * u(i,1,k) &
2629  + dy(j,2) * u(i,2,k) &
2630  + dy(j,3) * u(i,3,k) &
2631  + dy(j,4) * u(i,4,k) &
2632  + dy(j,5) * u(i,5,k) &
2633  + dy(j,6) * u(i,6,k) &
2634  + dy(j,7) * u(i,7,k) &
2635  + dy(j,8) * u(i,8,k) &
2636  + dy(j,9) * u(i,9,k) &
2637  + dy(j,10) * u(i,10,k)
2638  end do
2639  end do
2640  end do
2641 
2642  do k = 1, lx
2643  do i = 1, lx*lx
2644  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2645  + dz(k,2) * u(i,1,2) &
2646  + dz(k,3) * u(i,1,3) &
2647  + dz(k,4) * u(i,1,4) &
2648  + dz(k,5) * u(i,1,5) &
2649  + dz(k,6) * u(i,1,6) &
2650  + dz(k,7) * u(i,1,7) &
2651  + dz(k,8) * u(i,1,8) &
2652  + dz(k,9) * u(i,1,9) &
2653  + dz(k,10) * u(i,1,10)
2654  end do
2655  end do
2656 
2657  do i = 1, lx * lx * lx
2658  ux(i,1,1) = w3(i,1,1) &
2659  * ( drdx(i,1,1) * ur(i,1,1) &
2660  + dsdx(i,1,1) * us(i,1,1) &
2661  + dtdx(i,1,1) * ut(i,1,1) )
2662  uy(i,1,1) = w3(i,1,1) &
2663  * ( dsdy(i,1,1) * us(i,1,1) &
2664  + drdy(i,1,1) * ur(i,1,1) &
2665  + dtdy(i,1,1) * ut(i,1,1) )
2666  uz(i,1,1) = w3(i,1,1) &
2667  * ( dtdz(i,1,1) * ut(i,1,1) &
2668  + drdz(i,1,1) * ur(i,1,1) &
2669  + dsdz(i,1,1) * us(i,1,1) )
2670  end do
2671  end subroutine cpu_opgrad_lx10_single
2672 
2673  subroutine cpu_opgrad_lx9_single(ux, uy, uz, u, dx, dy, dz, &
2674  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2675  integer, parameter :: lx = 9
2676  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2677  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2678  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2679  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2680  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2681  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2682  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2683  real(kind=rp) :: ur(lx, lx, lx)
2684  real(kind=rp) :: us(lx, lx, lx)
2685  real(kind=rp) :: ut(lx, lx, lx)
2686  integer :: i, j, k
2687 
2688  do j = 1, lx * lx
2689  do i = 1, lx
2690  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2691  + dx(i,2) * u(2,j,1) &
2692  + dx(i,3) * u(3,j,1) &
2693  + dx(i,4) * u(4,j,1) &
2694  + dx(i,5) * u(5,j,1) &
2695  + dx(i,6) * u(6,j,1) &
2696  + dx(i,7) * u(7,j,1) &
2697  + dx(i,8) * u(8,j,1) &
2698  + dx(i,9) * u(9,j,1)
2699  end do
2700  end do
2701 
2702  do k = 1, lx
2703  do j = 1, lx
2704  do i = 1, lx
2705  us(i,j,k) = dy(j,1) * u(i,1,k) &
2706  + dy(j,2) * u(i,2,k) &
2707  + dy(j,3) * u(i,3,k) &
2708  + dy(j,4) * u(i,4,k) &
2709  + dy(j,5) * u(i,5,k) &
2710  + dy(j,6) * u(i,6,k) &
2711  + dy(j,7) * u(i,7,k) &
2712  + dy(j,8) * u(i,8,k) &
2713  + dy(j,9) * u(i,9,k)
2714  end do
2715  end do
2716  end do
2717 
2718  do k = 1, lx
2719  do i = 1, lx*lx
2720  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2721  + dz(k,2) * u(i,1,2) &
2722  + dz(k,3) * u(i,1,3) &
2723  + dz(k,4) * u(i,1,4) &
2724  + dz(k,5) * u(i,1,5) &
2725  + dz(k,6) * u(i,1,6) &
2726  + dz(k,7) * u(i,1,7) &
2727  + dz(k,8) * u(i,1,8) &
2728  + dz(k,9) * u(i,1,9)
2729  end do
2730  end do
2731 
2732  do i = 1, lx * lx * lx
2733  ux(i,1,1) = w3(i,1,1) &
2734  * ( drdx(i,1,1) * ur(i,1,1) &
2735  + dsdx(i,1,1) * us(i,1,1) &
2736  + dtdx(i,1,1) * ut(i,1,1) )
2737  uy(i,1,1) = w3(i,1,1) &
2738  * ( dsdy(i,1,1) * us(i,1,1) &
2739  + drdy(i,1,1) * ur(i,1,1) &
2740  + dtdy(i,1,1) * ut(i,1,1) )
2741  uz(i,1,1) = w3(i,1,1) &
2742  * ( dtdz(i,1,1) * ut(i,1,1) &
2743  + drdz(i,1,1) * ur(i,1,1) &
2744  + dsdz(i,1,1) * us(i,1,1) )
2745  end do
2746  end subroutine cpu_opgrad_lx9_single
2747 
2748  subroutine cpu_opgrad_lx8_single(ux, uy, uz, u, dx, dy, dz, &
2749  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2750  integer, parameter :: lx = 8
2751  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2752  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2753  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2754  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2755  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2756  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2757  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2758  real(kind=rp) :: ur(lx, lx, lx)
2759  real(kind=rp) :: us(lx, lx, lx)
2760  real(kind=rp) :: ut(lx, lx, lx)
2761  integer :: i, j, k
2762 
2763  do j = 1, lx * lx
2764  do i = 1, lx
2765  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2766  + dx(i,2) * u(2,j,1) &
2767  + dx(i,3) * u(3,j,1) &
2768  + dx(i,4) * u(4,j,1) &
2769  + dx(i,5) * u(5,j,1) &
2770  + dx(i,6) * u(6,j,1) &
2771  + dx(i,7) * u(7,j,1) &
2772  + dx(i,8) * u(8,j,1)
2773  end do
2774  end do
2775 
2776  do k = 1, lx
2777  do j = 1, lx
2778  do i = 1, lx
2779  us(i,j,k) = dy(j,1) * u(i,1,k) &
2780  + dy(j,2) * u(i,2,k) &
2781  + dy(j,3) * u(i,3,k) &
2782  + dy(j,4) * u(i,4,k) &
2783  + dy(j,5) * u(i,5,k) &
2784  + dy(j,6) * u(i,6,k) &
2785  + dy(j,7) * u(i,7,k) &
2786  + dy(j,8) * u(i,8,k)
2787  end do
2788  end do
2789  end do
2790 
2791  do k = 1, lx
2792  do i = 1, lx*lx
2793  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2794  + dz(k,2) * u(i,1,2) &
2795  + dz(k,3) * u(i,1,3) &
2796  + dz(k,4) * u(i,1,4) &
2797  + dz(k,5) * u(i,1,5) &
2798  + dz(k,6) * u(i,1,6) &
2799  + dz(k,7) * u(i,1,7) &
2800  + dz(k,8) * u(i,1,8)
2801  end do
2802  end do
2803 
2804  do i = 1, lx * lx * lx
2805  ux(i,1,1) = w3(i,1,1) &
2806  * ( drdx(i,1,1) * ur(i,1,1) &
2807  + dsdx(i,1,1) * us(i,1,1) &
2808  + dtdx(i,1,1) * ut(i,1,1) )
2809  uy(i,1,1) = w3(i,1,1) &
2810  * ( dsdy(i,1,1) * us(i,1,1) &
2811  + drdy(i,1,1) * ur(i,1,1) &
2812  + dtdy(i,1,1) * ut(i,1,1) )
2813  uz(i,1,1) = w3(i,1,1) &
2814  * ( dtdz(i,1,1) * ut(i,1,1) &
2815  + drdz(i,1,1) * ur(i,1,1) &
2816  + dsdz(i,1,1) * us(i,1,1) )
2817  end do
2818  end subroutine cpu_opgrad_lx8_single
2819 
2820  subroutine cpu_opgrad_lx7_single(ux, uy, uz, u, dx, dy, dz, &
2821  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2822  integer, parameter :: lx = 7
2823  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2824  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2825  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2826  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2827  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2828  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2829  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2830  real(kind=rp) :: ur(lx, lx, lx)
2831  real(kind=rp) :: us(lx, lx, lx)
2832  real(kind=rp) :: ut(lx, lx, lx)
2833  integer :: i, j, k
2834 
2835  do j = 1, lx * lx
2836  do i = 1, lx
2837  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2838  + dx(i,2) * u(2,j,1) &
2839  + dx(i,3) * u(3,j,1) &
2840  + dx(i,4) * u(4,j,1) &
2841  + dx(i,5) * u(5,j,1) &
2842  + dx(i,6) * u(6,j,1) &
2843  + dx(i,7) * u(7,j,1)
2844  end do
2845  end do
2846 
2847  do k = 1, lx
2848  do j = 1, lx
2849  do i = 1, lx
2850  us(i,j,k) = dy(j,1) * u(i,1,k) &
2851  + dy(j,2) * u(i,2,k) &
2852  + dy(j,3) * u(i,3,k) &
2853  + dy(j,4) * u(i,4,k) &
2854  + dy(j,5) * u(i,5,k) &
2855  + dy(j,6) * u(i,6,k) &
2856  + dy(j,7) * u(i,7,k)
2857  end do
2858  end do
2859  end do
2860 
2861  do k = 1, lx
2862  do i = 1, lx*lx
2863  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2864  + dz(k,2) * u(i,1,2) &
2865  + dz(k,3) * u(i,1,3) &
2866  + dz(k,4) * u(i,1,4) &
2867  + dz(k,5) * u(i,1,5) &
2868  + dz(k,6) * u(i,1,6) &
2869  + dz(k,7) * u(i,1,7)
2870  end do
2871  end do
2872 
2873  do i = 1, lx * lx * lx
2874  ux(i,1,1) = w3(i,1,1) &
2875  * ( drdx(i,1,1) * ur(i,1,1) &
2876  + dsdx(i,1,1) * us(i,1,1) &
2877  + dtdx(i,1,1) * ut(i,1,1) )
2878  uy(i,1,1) = w3(i,1,1) &
2879  * ( dsdy(i,1,1) * us(i,1,1) &
2880  + drdy(i,1,1) * ur(i,1,1) &
2881  + dtdy(i,1,1) * ut(i,1,1) )
2882  uz(i,1,1) = w3(i,1,1) &
2883  * ( dtdz(i,1,1) * ut(i,1,1) &
2884  + drdz(i,1,1) * ur(i,1,1) &
2885  + dsdz(i,1,1) * us(i,1,1) )
2886  end do
2887  end subroutine cpu_opgrad_lx7_single
2888 
2889  subroutine cpu_opgrad_lx6_single(ux, uy, uz, u, dx, dy, dz, &
2890  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2891  integer, parameter :: lx = 6
2892  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2893  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2894  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2895  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2896  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2897  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2898  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2899  real(kind=rp) :: ur(lx, lx, lx)
2900  real(kind=rp) :: us(lx, lx, lx)
2901  real(kind=rp) :: ut(lx, lx, lx)
2902  integer :: i, j, k
2903 
2904  do j = 1, lx * lx
2905  do i = 1, lx
2906  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2907  + dx(i,2) * u(2,j,1) &
2908  + dx(i,3) * u(3,j,1) &
2909  + dx(i,4) * u(4,j,1) &
2910  + dx(i,5) * u(5,j,1) &
2911  + dx(i,6) * u(6,j,1)
2912  end do
2913  end do
2914 
2915  do k = 1, lx
2916  do j = 1, lx
2917  do i = 1, lx
2918  us(i,j,k) = dy(j,1) * u(i,1,k) &
2919  + dy(j,2) * u(i,2,k) &
2920  + dy(j,3) * u(i,3,k) &
2921  + dy(j,4) * u(i,4,k) &
2922  + dy(j,5) * u(i,5,k) &
2923  + dy(j,6) * u(i,6,k)
2924  end do
2925  end do
2926  end do
2927 
2928  do k = 1, lx
2929  do i = 1, lx*lx
2930  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2931  + dz(k,2) * u(i,1,2) &
2932  + dz(k,3) * u(i,1,3) &
2933  + dz(k,4) * u(i,1,4) &
2934  + dz(k,5) * u(i,1,5) &
2935  + dz(k,6) * u(i,1,6)
2936  end do
2937  end do
2938 
2939  do i = 1, lx * lx * lx
2940  ux(i,1,1) = w3(i,1,1) &
2941  * ( drdx(i,1,1) * ur(i,1,1) &
2942  + dsdx(i,1,1) * us(i,1,1) &
2943  + dtdx(i,1,1) * ut(i,1,1) )
2944  uy(i,1,1) = w3(i,1,1) &
2945  * ( dsdy(i,1,1) * us(i,1,1) &
2946  + drdy(i,1,1) * ur(i,1,1) &
2947  + dtdy(i,1,1) * ut(i,1,1) )
2948  uz(i,1,1) = w3(i,1,1) &
2949  * ( dtdz(i,1,1) * ut(i,1,1) &
2950  + drdz(i,1,1) * ur(i,1,1) &
2951  + dsdz(i,1,1) * us(i,1,1) )
2952  end do
2953  end subroutine cpu_opgrad_lx6_single
2954 
2955  subroutine cpu_opgrad_lx5_single(ux, uy, uz, u, dx, dy, dz, &
2956  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
2957  integer, parameter :: lx = 5
2958  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
2959  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
2960  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
2961  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
2962  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
2963  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
2964  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
2965  real(kind=rp) :: ur(lx, lx, lx)
2966  real(kind=rp) :: us(lx, lx, lx)
2967  real(kind=rp) :: ut(lx, lx, lx)
2968  integer :: i, j, k
2969 
2970  do j = 1, lx * lx
2971  do i = 1, lx
2972  ur(i,j,1) = dx(i,1) * u(1,j,1) &
2973  + dx(i,2) * u(2,j,1) &
2974  + dx(i,3) * u(3,j,1) &
2975  + dx(i,4) * u(4,j,1) &
2976  + dx(i,5) * u(5,j,1)
2977  end do
2978  end do
2979 
2980  do k = 1, lx
2981  do j = 1, lx
2982  do i = 1, lx
2983  us(i,j,k) = dy(j,1) * u(i,1,k) &
2984  + dy(j,2) * u(i,2,k) &
2985  + dy(j,3) * u(i,3,k) &
2986  + dy(j,4) * u(i,4,k) &
2987  + dy(j,5) * u(i,5,k)
2988  end do
2989  end do
2990  end do
2991 
2992  do k = 1, lx
2993  do i = 1, lx*lx
2994  ut(i,1,k) = dz(k,1) * u(i,1,1) &
2995  + dz(k,2) * u(i,1,2) &
2996  + dz(k,3) * u(i,1,3) &
2997  + dz(k,4) * u(i,1,4) &
2998  + dz(k,5) * u(i,1,5)
2999  end do
3000  end do
3001 
3002  do i = 1, lx * lx * lx
3003  ux(i,1,1) = w3(i,1,1) &
3004  * ( drdx(i,1,1) * ur(i,1,1) &
3005  + dsdx(i,1,1) * us(i,1,1) &
3006  + dtdx(i,1,1) * ut(i,1,1) )
3007  uy(i,1,1) = w3(i,1,1) &
3008  * ( dsdy(i,1,1) * us(i,1,1) &
3009  + drdy(i,1,1) * ur(i,1,1) &
3010  + dtdy(i,1,1) * ut(i,1,1) )
3011  uz(i,1,1) = w3(i,1,1) &
3012  * ( dtdz(i,1,1) * ut(i,1,1) &
3013  + drdz(i,1,1) * ur(i,1,1) &
3014  + dsdz(i,1,1) * us(i,1,1) )
3015  end do
3016  end subroutine cpu_opgrad_lx5_single
3017 
3018  subroutine cpu_opgrad_lx4_single(ux, uy, uz, u, dx, dy, dz, &
3019  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
3020  integer, parameter :: lx = 4
3021  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
3022  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
3023  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
3024  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
3025  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
3026  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
3027  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
3028  real(kind=rp) :: ur(lx, lx, lx)
3029  real(kind=rp) :: us(lx, lx, lx)
3030  real(kind=rp) :: ut(lx, lx, lx)
3031  integer :: i, j, k
3032 
3033  do j = 1, lx * lx
3034  do i = 1, lx
3035  ur(i,j,1) = dx(i,1) * u(1,j,1) &
3036  + dx(i,2) * u(2,j,1) &
3037  + dx(i,3) * u(3,j,1) &
3038  + dx(i,4) * u(4,j,1)
3039  end do
3040  end do
3041 
3042  do k = 1, lx
3043  do j = 1, lx
3044  do i = 1, lx
3045  us(i,j,k) = dy(j,1) * u(i,1,k) &
3046  + dy(j,2) * u(i,2,k) &
3047  + dy(j,3) * u(i,3,k) &
3048  + dy(j,4) * u(i,4,k)
3049  end do
3050  end do
3051  end do
3052 
3053  do k = 1, lx
3054  do i = 1, lx*lx
3055  ut(i,1,k) = dz(k,1) * u(i,1,1) &
3056  + dz(k,2) * u(i,1,2) &
3057  + dz(k,3) * u(i,1,3) &
3058  + dz(k,4) * u(i,1,4)
3059  end do
3060  end do
3061 
3062  do i = 1, lx * lx * lx
3063  ux(i,1,1) = w3(i,1,1) &
3064  * ( drdx(i,1,1) * ur(i,1,1) &
3065  + dsdx(i,1,1) * us(i,1,1) &
3066  + dtdx(i,1,1) * ut(i,1,1) )
3067  uy(i,1,1) = w3(i,1,1) &
3068  * ( dsdy(i,1,1) * us(i,1,1) &
3069  + drdy(i,1,1) * ur(i,1,1) &
3070  + dtdy(i,1,1) * ut(i,1,1) )
3071  uz(i,1,1) = w3(i,1,1) &
3072  * ( dtdz(i,1,1) * ut(i,1,1) &
3073  + drdz(i,1,1) * ur(i,1,1) &
3074  + dsdz(i,1,1) * us(i,1,1) )
3075  end do
3076  end subroutine cpu_opgrad_lx4_single
3077 
3078  subroutine cpu_opgrad_lx3_single(ux, uy, uz, u, dx, dy, dz, &
3079  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
3080  integer, parameter :: lx = 3
3081  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
3082  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
3083  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
3084  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
3085  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
3086  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
3087  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
3088  real(kind=rp) :: ur(lx, lx, lx)
3089  real(kind=rp) :: us(lx, lx, lx)
3090  real(kind=rp) :: ut(lx, lx, lx)
3091  integer :: i, j, k
3092 
3093  do j = 1, lx * lx
3094  do i = 1, lx
3095  ur(i,j,1) = dx(i,1) * u(1,j,1) &
3096  + dx(i,2) * u(2,j,1) &
3097  + dx(i,3) * u(3,j,1)
3098  end do
3099  end do
3100 
3101  do k = 1, lx
3102  do j = 1, lx
3103  do i = 1, lx
3104  us(i,j,k) = dy(j,1) * u(i,1,k) &
3105  + dy(j,2) * u(i,2,k) &
3106  + dy(j,3) * u(i,3,k)
3107  end do
3108  end do
3109  end do
3110 
3111  do k = 1, lx
3112  do i = 1, lx*lx
3113  ut(i,1,k) = dz(k,1) * u(i,1,1) &
3114  + dz(k,2) * u(i,1,2) &
3115  + dz(k,3) * u(i,1,3)
3116  end do
3117  end do
3118 
3119  do i = 1, lx * lx * lx
3120  ux(i,1,1) = w3(i,1,1) &
3121  * ( drdx(i,1,1) * ur(i,1,1) &
3122  + dsdx(i,1,1) * us(i,1,1) &
3123  + dtdx(i,1,1) * ut(i,1,1) )
3124  uy(i,1,1) = w3(i,1,1) &
3125  * ( dsdy(i,1,1) * us(i,1,1) &
3126  + drdy(i,1,1) * ur(i,1,1) &
3127  + dtdy(i,1,1) * ut(i,1,1) )
3128  uz(i,1,1) = w3(i,1,1) &
3129  * ( dtdz(i,1,1) * ut(i,1,1) &
3130  + drdz(i,1,1) * ur(i,1,1) &
3131  + dsdz(i,1,1) * us(i,1,1) )
3132  end do
3133  end subroutine cpu_opgrad_lx3_single
3134 
3135  subroutine cpu_opgrad_lx2_single(ux, uy, uz, u, dx, dy, dz, &
3136  drdx, dsdx, dtdx, drdy, dsdy, dtdy, drdz, dsdz, dtdz, w3)
3137  integer, parameter :: lx = 2
3138  real(kind=rp), dimension(lx, lx, lx), intent(inout) :: ux, uy, uz
3139  real(kind=rp), dimension(lx, lx, lx), intent(in) :: u
3140  real(kind=rp), dimension(lx, lx), intent(in) :: dx, dy, dz
3141  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdx, dsdx, dtdx
3142  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdy, dsdy, dtdy
3143  real(kind=rp), dimension(lx, lx, lx), intent(in) :: drdz, dsdz, dtdz
3144  real(kind=rp), dimension(lx, lx, lx), intent(in) :: w3
3145  real(kind=rp) :: ur(lx, lx, lx)
3146  real(kind=rp) :: us(lx, lx, lx)
3147  real(kind=rp) :: ut(lx, lx, lx)
3148  integer :: i, j, k
3149 
3150  do j = 1, lx * lx
3151  do i = 1, lx
3152  ur(i,j,1) = dx(i,1) * u(1,j,1) &
3153  + dx(i,2) * u(2,j,1)
3154  end do
3155  end do
3156 
3157  do k = 1, lx
3158  do j = 1, lx
3159  do i = 1, lx
3160  us(i,j,k) = dy(j,1) * u(i,1,k) &
3161  + dy(j,2) * u(i,2,k)
3162  end do
3163  end do
3164  end do
3165 
3166  do k = 1, lx
3167  do i = 1, lx*lx
3168  ut(i,1,k) = dz(k,1) * u(i,1,1) &
3169  + dz(k,2) * u(i,1,2)
3170  end do
3171  end do
3172 
3173  do i = 1, lx * lx * lx
3174  ux(i,1,1) = w3(i,1,1) &
3175  * ( drdx(i,1,1) * ur(i,1,1) &
3176  + dsdx(i,1,1) * us(i,1,1) &
3177  + dtdx(i,1,1) * ut(i,1,1) )
3178  uy(i,1,1) = w3(i,1,1) &
3179  * ( dsdy(i,1,1) * us(i,1,1) &
3180  + drdy(i,1,1) * ur(i,1,1) &
3181  + dtdy(i,1,1) * ut(i,1,1) )
3182  uz(i,1,1) = w3(i,1,1) &
3183  * ( dtdz(i,1,1) * ut(i,1,1) &
3184  + drdz(i,1,1) * ur(i,1,1) &
3185  + dsdz(i,1,1) * us(i,1,1) )
3186  end do
3187 
3188  end subroutine cpu_opgrad_lx2_single
3189 
3190 end submodule cpu_opgrad
opr_cpu
Operators CPU backend.
Definition: opr_cpu.f90:34