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