Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
pc_jacobi_sx.f90
Go to the documentation of this file.
1! Copyright (c) 2021, 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!
34module sx_jacobi
35 use math
36 use precon
37 use coefs
38 use dofmap
39 use num_types
40 use gather_scatter
41 implicit none
42 private
43
45 type, public, extends(pc_t) :: sx_jacobi_t
46 real(kind=rp), allocatable :: d(:, :, :, :)
47 type(gs_t), pointer :: gs_h
48 type(dofmap_t), pointer :: dof
49 type(coef_t), pointer :: coef
50 contains
51 procedure, pass(this) :: init => sx_jacobi_init
52 procedure, pass(this) :: free => sx_jacobi_free
53 procedure, pass(this) :: solve => sx_jacobi_solve
54 procedure, pass(this) :: update => sx_jacobi_update
55 end type sx_jacobi_t
56
57contains
58
59 subroutine sx_jacobi_init(this, coef, dof, gs_h)
60 class(sx_jacobi_t), intent(inout) :: this
61 type(coef_t), intent(in), target :: coef
62 type(dofmap_t), intent(in), target :: dof
63 type(gs_t), intent(in), target :: gs_h
64
65 call this%free()
66 this%gs_h => gs_h
67 this%dof => dof
68 this%coef => coef
69 allocate(this%d(dof%Xh%lx, dof%Xh%ly, dof%Xh%lz, dof%msh%nelv))
70 call sx_jacobi_update(this)
71
72 end subroutine sx_jacobi_init
73
74 subroutine sx_jacobi_free(this)
75 class(sx_jacobi_t), intent(inout) :: this
76 if (allocated(this%d)) then
77 deallocate(this%d)
78 end if
79 nullify(this%dof)
80 nullify(this%gs_h)
81 nullify(this%coef)
82 end subroutine sx_jacobi_free
83
86 subroutine sx_jacobi_solve(this, z, r, n)
87 integer, intent(in) :: n
88 class(sx_jacobi_t), intent(inout) :: this
89 real(kind=rp), dimension(n), intent(inout) :: z
90 real(kind=rp), dimension(n), intent(inout) :: r
91 call col3(z, r, this%d, n)
92 end subroutine sx_jacobi_solve
93
94 subroutine sx_jacobi_update(this)
95 class(sx_jacobi_t), intent(inout) :: this
96 integer :: lz, ly, lx
97 associate(dof => this%dof, coef => this%coef, &
98 gs_h => this%gs_h, nelv => this%dof%msh%nelv)
99
100 lx = dof%Xh%lx
101 ly = dof%Xh%ly
102 lz = dof%Xh%lz
103
104 this%d = 0d0
105
106 select case(lx)
107 case (14)
108 call sx_update_lx13(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
109 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
110 case (13)
111 call sx_update_lx13(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
112 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
113 case (12)
114 call sx_update_lx12(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
115 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
116 case (11)
117 call sx_update_lx11(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
118 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
119 case (10)
120 call sx_update_lx10(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
121 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
122 case (9)
123 call sx_update_lx9(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
124 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
125 case (8)
126 call sx_update_lx8(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
127 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
128 case (7)
129 call sx_update_lx7(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
130 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
131 case (6)
132 call sx_update_lx6(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
133 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
134 case (5)
135 call sx_update_lx5(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
136 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
137 case (4)
138 call sx_update_lx4(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
139 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
140 case (2)
141 call sx_update_lx2(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
142 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, nelv)
143 case default
144 call sx_update_lx(this%d, coef%Xh%dxt, coef%Xh%dyt, coef%Xh%dzt, &
145 coef%G11, coef%G22, coef%G33, coef%G12, coef%G13, coef%G23, &
146 nelv, lx)
147 end select
148
149 call col2(this%d, coef%h1, coef%dof%size())
150 if (coef%ifh2) call addcol3(this%d, coef%h2, coef%B, coef%dof%size())
151 call gs_h%op(this%d, dof%size(), gs_op_add)
152 if (.not. coef%ifh2) call col2(this%d, coef%mult, coef%dof%size())
153 call invcol1(this%d, dof%size())
154 end associate
155 end subroutine sx_jacobi_update
156
157 subroutine sx_update_lx(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n, lx)
158 integer, intent(in) :: n, lx
159 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
160 real(kind=rp), intent(in) :: dxt(lx, lx)
161 real(kind=rp), intent(in) :: dyt(lx, lx)
162 real(kind=rp), intent(in) :: dzt(lx, lx)
163 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
164 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
165 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
166 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
167 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
168 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
169 integer :: i, j, k, l, e
170
171
172 do l = 1, lx
173 do k = 1, lx
174 do j = 1, lx
175 do i = 1, lx
176 do e = 1, n
177 d(i, j, k, e) = d(i, j, k, e) + &
178 g11(l, j, k, e) * dxt(i, l)**2
179
180 d(i, j, k, e) = d(i, j, k, e) + &
181 g22(i, l, k, e) * dyt(j, l)**2
182
183 d(i, j, k, e) = d(i, j, k, e) + &
184 g33(i, j, l, e) * dzt(k, l)**2
185
186 end do
187 end do
188 end do
189 end do
190 end do
191
192 do j = 1, lx, lx-1
193 do k = 1, lx, lx-1
194 do e = 1, n
195 d(1, j, k, e) = d(1, j, k, e) &
196 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
197 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
198 d(lx, j, k, e) = d(lx, j, k, e) &
199 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
200 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
201 end do
202 end do
203 end do
204
205 do i = 1, lx, lx-1
206 do k = 1, lx, lx-1
207 do e = 1, n
208 d(i, 1, k, e) = d(i, 1, k, e) &
209 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
210 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
211 d(i, lx, k, e) = d(i, lx, k, e) &
212 + g12(i, lx, k, e) * dyt(lx, lx)*dxt(i, i) &
213 + g23(i, lx, k, e) * dyt(lx, lx)*dzt(k, k)
214 end do
215 end do
216 end do
217
218 do i = 1, lx, lx-1
219 do j = 1, lx, lx-1
220 do e = 1, n
221 d(i, j, 1, e) = d(i, j, 1, e) &
222 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
223 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
224 d(i, j, lx, e) = d(i, j, lx, e) &
225 + g13(i, j, lx, e) * dzt(lx, lx)*dxt(i, i) &
226 + g23(i, j, lx, e) * dzt(lx, lx)*dyt(j, j)
227 end do
228 end do
229 end do
230
231 end subroutine sx_update_lx
232
233 subroutine sx_update_lx14(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
234 integer, parameter :: lx = 14
235 integer, parameter :: ly = 14
236 integer, parameter :: lz = 14
237 integer, intent(in) :: n
238 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
239 real(kind=rp), intent(in) :: dxt(lx, lx)
240 real(kind=rp), intent(in) :: dyt(lx, lx)
241 real(kind=rp), intent(in) :: dzt(lx, lx)
242 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
243 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
244 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
245 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
246 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
247 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
248 integer :: i, j, k, l, e
249
250
251 do l = 1, lx
252 do k = 1, lz
253 do j = 1, ly
254 do i = 1, lx
255 do e = 1, n
256 d(i, j, k, e) = d(i, j, k, e) + &
257 g11(l, j, k, e) * dxt(i, l)**2
258
259 d(i, j, k, e) = d(i, j, k, e) + &
260 g22(i, l, k, e) * dyt(j, l)**2
261
262 d(i, j, k, e) = d(i, j, k, e) + &
263 g33(i, j, l, e) * dzt(k, l)**2
264
265 end do
266 end do
267 end do
268 end do
269 end do
270
271 do j = 1, ly, ly-1
272 do k = 1, lz, lz-1
273 do e = 1, n
274 d(1, j, k, e) = d(1, j, k, e) &
275 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
276 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
277 d(lx, j, k, e) = d(lx, j, k, e) &
278 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
279 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
280 end do
281 end do
282 end do
283
284 do i = 1, lx, lx-1
285 do k = 1, lz, lz-1
286 do e = 1, n
287 d(i, 1, k, e) = d(i, 1, k, e) &
288 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
289 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
290 d(i, ly, k, e) = d(i, ly, k, e) &
291 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
292 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
293 end do
294 end do
295 end do
296
297 do i = 1, lx, lx-1
298 do j = 1, ly, ly-1
299 do e = 1, n
300 d(i, j, 1, e) = d(i, j, 1, e) &
301 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
302 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
303 d(i, j, lz, e) = d(i, j, lz, e) &
304 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
305 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
306 end do
307 end do
308 end do
309
310 end subroutine sx_update_lx14
311
312 subroutine sx_update_lx13(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
313 integer, parameter :: lx = 13
314 integer, parameter :: ly = 13
315 integer, parameter :: lz = 13
316 integer, intent(in) :: n
317 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
318 real(kind=rp), intent(in) :: dxt(lx, lx)
319 real(kind=rp), intent(in) :: dyt(lx, lx)
320 real(kind=rp), intent(in) :: dzt(lx, lx)
321 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
322 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
323 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
324 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
325 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
326 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
327 integer :: i, j, k, l, e
328
329
330 do l = 1, lx
331 do k = 1, lz
332 do j = 1, ly
333 do i = 1, lx
334 do e = 1, n
335 d(i, j, k, e) = d(i, j, k, e) + &
336 g11(l, j, k, e) * dxt(i, l)**2
337
338 d(i, j, k, e) = d(i, j, k, e) + &
339 g22(i, l, k, e) * dyt(j, l)**2
340
341 d(i, j, k, e) = d(i, j, k, e) + &
342 g33(i, j, l, e) * dzt(k, l)**2
343
344 end do
345 end do
346 end do
347 end do
348 end do
349
350 do j = 1, ly, ly-1
351 do k = 1, lz, lz-1
352 do e = 1, n
353 d(1, j, k, e) = d(1, j, k, e) &
354 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
355 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
356 d(lx, j, k, e) = d(lx, j, k, e) &
357 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
358 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
359 end do
360 end do
361 end do
362
363 do i = 1, lx, lx-1
364 do k = 1, lz, lz-1
365 do e = 1, n
366 d(i, 1, k, e) = d(i, 1, k, e) &
367 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
368 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
369 d(i, ly, k, e) = d(i, ly, k, e) &
370 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
371 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
372 end do
373 end do
374 end do
375
376 do i = 1, lx, lx-1
377 do j = 1, ly, ly-1
378 do e = 1, n
379 d(i, j, 1, e) = d(i, j, 1, e) &
380 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
381 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
382 d(i, j, lz, e) = d(i, j, lz, e) &
383 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
384 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
385 end do
386 end do
387 end do
388
389 end subroutine sx_update_lx13
390
391 subroutine sx_update_lx12(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
392 integer, parameter :: lx = 12
393 integer, parameter :: ly = 12
394 integer, parameter :: lz = 12
395 integer, intent(in) :: n
396 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
397 real(kind=rp), intent(in) :: dxt(lx, lx)
398 real(kind=rp), intent(in) :: dyt(lx, lx)
399 real(kind=rp), intent(in) :: dzt(lx, lx)
400 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
401 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
402 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
403 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
404 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
405 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
406 integer :: i, j, k, l, e
407
408
409 do l = 1, lx
410 do k = 1, lz
411 do j = 1, ly
412 do i = 1, lx
413 do e = 1, n
414 d(i, j, k, e) = d(i, j, k, e) + &
415 g11(l, j, k, e) * dxt(i, l)**2
416
417 d(i, j, k, e) = d(i, j, k, e) + &
418 g22(i, l, k, e) * dyt(j, l)**2
419
420 d(i, j, k, e) = d(i, j, k, e) + &
421 g33(i, j, l, e) * dzt(k, l)**2
422
423 end do
424 end do
425 end do
426 end do
427 end do
428
429 do j = 1, ly, ly-1
430 do k = 1, lz, lz-1
431 do e = 1, n
432 d(1, j, k, e) = d(1, j, k, e) &
433 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
434 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
435 d(lx, j, k, e) = d(lx, j, k, e) &
436 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
437 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
438 end do
439 end do
440 end do
441
442 do i = 1, lx, lx-1
443 do k = 1, lz, lz-1
444 do e = 1, n
445 d(i, 1, k, e) = d(i, 1, k, e) &
446 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
447 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
448 d(i, ly, k, e) = d(i, ly, k, e) &
449 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
450 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
451 end do
452 end do
453 end do
454
455 do i = 1, lx, lx-1
456 do j = 1, ly, ly-1
457 do e = 1, n
458 d(i, j, 1, e) = d(i, j, 1, e) &
459 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
460 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
461 d(i, j, lz, e) = d(i, j, lz, e) &
462 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
463 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
464 end do
465 end do
466 end do
467
468 end subroutine sx_update_lx12
469
470 subroutine sx_update_lx11(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
471 integer, parameter :: lx = 11
472 integer, parameter :: ly = 11
473 integer, parameter :: lz = 11
474 integer, intent(in) :: n
475 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
476 real(kind=rp), intent(in) :: dxt(lx, lx)
477 real(kind=rp), intent(in) :: dyt(lx, lx)
478 real(kind=rp), intent(in) :: dzt(lx, lx)
479 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
480 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
481 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
482 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
483 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
484 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
485 integer :: i, j, k, l, e
486
487
488 do l = 1, lx
489 do k = 1, lz
490 do j = 1, ly
491 do i = 1, lx
492 do e = 1, n
493 d(i, j, k, e) = d(i, j, k, e) + &
494 g11(l, j, k, e) * dxt(i, l)**2
495
496 d(i, j, k, e) = d(i, j, k, e) + &
497 g22(i, l, k, e) * dyt(j, l)**2
498
499 d(i, j, k, e) = d(i, j, k, e) + &
500 g33(i, j, l, e) * dzt(k, l)**2
501
502 end do
503 end do
504 end do
505 end do
506 end do
507
508 do j = 1, ly, ly-1
509 do k = 1, lz, lz-1
510 do e = 1, n
511 d(1, j, k, e) = d(1, j, k, e) &
512 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
513 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
514 d(lx, j, k, e) = d(lx, j, k, e) &
515 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
516 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
517 end do
518 end do
519 end do
520
521 do i = 1, lx, lx-1
522 do k = 1, lz, lz-1
523 do e = 1, n
524 d(i, 1, k, e) = d(i, 1, k, e) &
525 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
526 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
527 d(i, ly, k, e) = d(i, ly, k, e) &
528 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
529 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
530 end do
531 end do
532 end do
533
534 do i = 1, lx, lx-1
535 do j = 1, ly, ly-1
536 do e = 1, n
537 d(i, j, 1, e) = d(i, j, 1, e) &
538 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
539 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
540 d(i, j, lz, e) = d(i, j, lz, e) &
541 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
542 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
543 end do
544 end do
545 end do
546
547 end subroutine sx_update_lx11
548
549 subroutine sx_update_lx10(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
550 integer, parameter :: lx = 10
551 integer, parameter :: ly = 10
552 integer, parameter :: lz = 10
553 integer, intent(in) :: n
554 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
555 real(kind=rp), intent(in) :: dxt(lx, lx)
556 real(kind=rp), intent(in) :: dyt(lx, lx)
557 real(kind=rp), intent(in) :: dzt(lx, lx)
558 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
559 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
560 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
561 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
562 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
563 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
564 integer :: i, j, k, l, e
565
566
567 do l = 1, lx
568 do k = 1, lz
569 do j = 1, ly
570 do i = 1, lx
571 do e = 1, n
572 d(i, j, k, e) = d(i, j, k, e) + &
573 g11(l, j, k, e) * dxt(i, l)**2
574
575 d(i, j, k, e) = d(i, j, k, e) + &
576 g22(i, l, k, e) * dyt(j, l)**2
577
578 d(i, j, k, e) = d(i, j, k, e) + &
579 g33(i, j, l, e) * dzt(k, l)**2
580
581 end do
582 end do
583 end do
584 end do
585 end do
586
587 do j = 1, ly, ly-1
588 do k = 1, lz, lz-1
589 do e = 1, n
590 d(1, j, k, e) = d(1, j, k, e) &
591 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
592 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
593 d(lx, j, k, e) = d(lx, j, k, e) &
594 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
595 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
596 end do
597 end do
598 end do
599
600 do i = 1, lx, lx-1
601 do k = 1, lz, lz-1
602 do e = 1, n
603 d(i, 1, k, e) = d(i, 1, k, e) &
604 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
605 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
606 d(i, ly, k, e) = d(i, ly, k, e) &
607 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
608 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
609 end do
610 end do
611 end do
612
613 do i = 1, lx, lx-1
614 do j = 1, ly, ly-1
615 do e = 1, n
616 d(i, j, 1, e) = d(i, j, 1, e) &
617 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
618 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
619 d(i, j, lz, e) = d(i, j, lz, e) &
620 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
621 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
622 end do
623 end do
624 end do
625
626 end subroutine sx_update_lx10
627
628 subroutine sx_update_lx9(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
629 integer, parameter :: lx = 9
630 integer, parameter :: ly = 9
631 integer, parameter :: lz = 9
632 integer, intent(in) :: n
633 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
634 real(kind=rp), intent(in) :: dxt(lx, lx)
635 real(kind=rp), intent(in) :: dyt(lx, lx)
636 real(kind=rp), intent(in) :: dzt(lx, lx)
637 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
638 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
639 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
640 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
641 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
642 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
643 integer :: i, j, k, l, e
644
645
646 do l = 1, lx
647 do k = 1, lz
648 do j = 1, ly
649 do i = 1, lx
650 do e = 1, n
651 d(i, j, k, e) = d(i, j, k, e) + &
652 g11(l, j, k, e) * dxt(i, l)**2
653
654 d(i, j, k, e) = d(i, j, k, e) + &
655 g22(i, l, k, e) * dyt(j, l)**2
656
657 d(i, j, k, e) = d(i, j, k, e) + &
658 g33(i, j, l, e) * dzt(k, l)**2
659
660 end do
661 end do
662 end do
663 end do
664 end do
665
666 do j = 1, ly, ly-1
667 do k = 1, lz, lz-1
668 do e = 1, n
669 d(1, j, k, e) = d(1, j, k, e) &
670 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
671 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
672 d(lx, j, k, e) = d(lx, j, k, e) &
673 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
674 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
675 end do
676 end do
677 end do
678
679 do i = 1, lx, lx-1
680 do k = 1, lz, lz-1
681 do e = 1, n
682 d(i, 1, k, e) = d(i, 1, k, e) &
683 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
684 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
685 d(i, ly, k, e) = d(i, ly, k, e) &
686 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
687 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
688 end do
689 end do
690 end do
691
692 do i = 1, lx, lx-1
693 do j = 1, ly, ly-1
694 do e = 1, n
695 d(i, j, 1, e) = d(i, j, 1, e) &
696 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
697 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
698 d(i, j, lz, e) = d(i, j, lz, e) &
699 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
700 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
701 end do
702 end do
703 end do
704
705 end subroutine sx_update_lx9
706
707 subroutine sx_update_lx8(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
708 integer, parameter :: lx = 8
709 integer, parameter :: ly = 8
710 integer, parameter :: lz = 8
711 integer, intent(in) :: n
712 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
713 real(kind=rp), intent(in) :: dxt(lx, lx)
714 real(kind=rp), intent(in) :: dyt(lx, lx)
715 real(kind=rp), intent(in) :: dzt(lx, lx)
716 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
717 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
718 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
719 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
720 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
721 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
722 integer :: i, j, k, l, e
723
724
725 do l = 1, lx
726 do k = 1, lz
727 do j = 1, ly
728 do i = 1, lx
729 do e = 1, n
730 d(i, j, k, e) = d(i, j, k, e) + &
731 g11(l, j, k, e) * dxt(i, l)**2
732
733 d(i, j, k, e) = d(i, j, k, e) + &
734 g22(i, l, k, e) * dyt(j, l)**2
735
736 d(i, j, k, e) = d(i, j, k, e) + &
737 g33(i, j, l, e) * dzt(k, l)**2
738
739 end do
740 end do
741 end do
742 end do
743 end do
744
745 do j = 1, ly, ly-1
746 do k = 1, lz, lz-1
747 do e = 1, n
748 d(1, j, k, e) = d(1, j, k, e) &
749 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
750 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
751 d(lx, j, k, e) = d(lx, j, k, e) &
752 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
753 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
754 end do
755 end do
756 end do
757
758 do i = 1, lx, lx-1
759 do k = 1, lz, lz-1
760 do e = 1, n
761 d(i, 1, k, e) = d(i, 1, k, e) &
762 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
763 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
764 d(i, ly, k, e) = d(i, ly, k, e) &
765 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
766 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
767 end do
768 end do
769 end do
770
771 do i = 1, lx, lx-1
772 do j = 1, ly, ly-1
773 do e = 1, n
774 d(i, j, 1, e) = d(i, j, 1, e) &
775 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
776 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
777 d(i, j, lz, e) = d(i, j, lz, e) &
778 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
779 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
780 end do
781 end do
782 end do
783
784 end subroutine sx_update_lx8
785
786 subroutine sx_update_lx7(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
787 integer, parameter :: lx = 7
788 integer, parameter :: ly = 7
789 integer, parameter :: lz = 7
790 integer, intent(in) :: n
791 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
792 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
793 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
794 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
795 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
796 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
797 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
798 real(kind=rp), intent(in) :: dxt(lx, lx)
799 real(kind=rp), intent(in) :: dyt(lx, lx)
800 real(kind=rp), intent(in) :: dzt(lx, lx)
801 integer :: i, j, k, l, e
802
803 do l = 1, lx
804 do k = 1, lz
805 do j = 1, ly
806 do i = 1, lx
807 do e = 1, n
808 d(i, j, k, e) = d(i, j, k, e) + &
809 g11(l, j, k, e) * dxt(i, l)**2
810
811 d(i, j, k, e) = d(i, j, k, e) + &
812 g22(i, l, k, e) * dyt(j, l)**2
813
814 d(i, j, k, e) = d(i, j, k, e) + &
815 g33(i, j, l, e) * dzt(k, l)**2
816
817 end do
818 end do
819 end do
820 end do
821 end do
822
823 do j = 1, ly, ly-1
824 do k = 1, lz, lz-1
825 do e = 1, n
826 d(1, j, k, e) = d(1, j, k, e) &
827 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
828 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
829 d(lx, j, k, e) = d(lx, j, k, e) &
830 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
831 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
832 end do
833 end do
834 end do
835
836 do i = 1, lx, lx-1
837 do k = 1, lz, lz-1
838 do e = 1, n
839 d(i, 1, k, e) = d(i, 1, k, e) &
840 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
841 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
842 d(i, ly, k, e) = d(i, ly, k, e) &
843 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
844 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
845 end do
846 end do
847 end do
848
849 do i = 1, lx, lx-1
850 do j = 1, ly, ly-1
851 do e = 1, n
852 d(i, j, 1, e) = d(i, j, 1, e) &
853 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
854 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
855 d(i, j, lz, e) = d(i, j, lz, e) &
856 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
857 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
858 end do
859 end do
860 end do
861
862
863 end subroutine sx_update_lx7
864
865 subroutine sx_update_lx6(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
866 integer, parameter :: lx = 6
867 integer, parameter :: ly = 6
868 integer, parameter :: lz = 6
869 integer, intent(in) :: n
870 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
871 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
872 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
873 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
874 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
875 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
876 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
877 real(kind=rp), intent(in) :: dxt(lx, lx)
878 real(kind=rp), intent(in) :: dyt(lx, lx)
879 real(kind=rp), intent(in) :: dzt(lx, lx)
880 integer :: i, j, k, l, e
881
882 do l = 1, lx
883 do k = 1, lz
884 do j = 1, ly
885 do i = 1, lx
886 do e = 1, n
887 d(i, j, k, e) = d(i, j, k, e) + &
888 g11(l, j, k, e) * dxt(i, l)**2
889
890 d(i, j, k, e) = d(i, j, k, e) + &
891 g22(i, l, k, e) * dyt(j, l)**2
892
893 d(i, j, k, e) = d(i, j, k, e) + &
894 g33(i, j, l, e) * dzt(k, l)**2
895
896 end do
897 end do
898 end do
899 end do
900 end do
901
902 do j = 1, ly, ly-1
903 do k = 1, lz, lz-1
904 do e = 1, n
905 d(1, j, k, e) = d(1, j, k, e) &
906 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
907 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
908 d(lx, j, k, e) = d(lx, j, k, e) &
909 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
910 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
911 end do
912 end do
913 end do
914
915 do i = 1, lx, lx-1
916 do k = 1, lz, lz-1
917 do e = 1, n
918 d(i, 1, k, e) = d(i, 1, k, e) &
919 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
920 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
921 d(i, ly, k, e) = d(i, ly, k, e) &
922 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
923 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
924 end do
925 end do
926 end do
927
928 do i = 1, lx, lx-1
929 do j = 1, ly, ly-1
930 do e = 1, n
931 d(i, j, 1, e) = d(i, j, 1, e) &
932 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
933 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
934 d(i, j, lz, e) = d(i, j, lz, e) &
935 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
936 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
937 end do
938 end do
939 end do
940
941
942 end subroutine sx_update_lx6
943
944 subroutine sx_update_lx5(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
945 integer, parameter :: lx = 5
946 integer, parameter :: ly = 5
947 integer, parameter :: lz = 5
948 integer, intent(in) :: n
949 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
950 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
951 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
952 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
953 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
954 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
955 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
956 real(kind=rp), intent(in) :: dxt(lx, lx)
957 real(kind=rp), intent(in) :: dyt(lx, lx)
958 real(kind=rp), intent(in) :: dzt(lx, lx)
959 integer :: i, j, k, l, e
960
961 do l = 1, lx
962 do k = 1, lz
963 do j = 1, ly
964 do i = 1, lx
965 do e = 1, n
966 d(i, j, k, e) = d(i, j, k, e) + &
967 g11(l, j, k, e) * dxt(i, l)**2
968
969 d(i, j, k, e) = d(i, j, k, e) + &
970 g22(i, l, k, e) * dyt(j, l)**2
971
972 d(i, j, k, e) = d(i, j, k, e) + &
973 g33(i, j, l, e) * dzt(k, l)**2
974
975 end do
976 end do
977 end do
978 end do
979 end do
980
981 do j = 1, ly, ly-1
982 do k = 1, lz, lz-1
983 do e = 1, n
984 d(1, j, k, e) = d(1, j, k, e) &
985 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
986 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
987 d(lx, j, k, e) = d(lx, j, k, e) &
988 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
989 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
990 end do
991 end do
992 end do
993
994 do i = 1, lx, lx-1
995 do k = 1, lz, lz-1
996 do e = 1, n
997 d(i, 1, k, e) = d(i, 1, k, e) &
998 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
999 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1000 d(i, ly, k, e) = d(i, ly, k, e) &
1001 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1002 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1003 end do
1004 end do
1005 end do
1006
1007 do i = 1, lx, lx-1
1008 do j = 1, ly, ly-1
1009 do e = 1, n
1010 d(i, j, 1, e) = d(i, j, 1, e) &
1011 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1012 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1013 d(i, j, lz, e) = d(i, j, lz, e) &
1014 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1015 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1016 end do
1017 end do
1018 end do
1019
1020
1021 end subroutine sx_update_lx5
1022
1023 subroutine sx_update_lx4(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
1024 integer, parameter :: lx = 4
1025 integer, parameter :: ly = 4
1026 integer, parameter :: lz = 4
1027 integer, intent(in) :: n
1028 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1029 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1030 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1031 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1032 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1033 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1034 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1035 real(kind=rp), intent(in) :: dxt(lx, lx)
1036 real(kind=rp), intent(in) :: dyt(lx, lx)
1037 real(kind=rp), intent(in) :: dzt(lx, lx)
1038 integer :: i, j, k, l, e
1039
1040 do l = 1, lx
1041 do k = 1, lz
1042 do j = 1, ly
1043 do i = 1, lx
1044 do e = 1, n
1045 d(i, j, k, e) = d(i, j, k, e) + &
1046 g11(l, j, k, e) * dxt(i, l)**2
1047
1048 d(i, j, k, e) = d(i, j, k, e) + &
1049 g22(i, l, k, e) * dyt(j, l)**2
1050
1051 d(i, j, k, e) = d(i, j, k, e) + &
1052 g33(i, j, l, e) * dzt(k, l)**2
1053
1054 end do
1055 end do
1056 end do
1057 end do
1058 end do
1059
1060 do j = 1, ly, ly-1
1061 do k = 1, lz, lz-1
1062 do e = 1, n
1063 d(1, j, k, e) = d(1, j, k, e) &
1064 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1065 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1066 d(lx, j, k, e) = d(lx, j, k, e) &
1067 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1068 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1069 end do
1070 end do
1071 end do
1072
1073 do i = 1, lx, lx-1
1074 do k = 1, lz, lz-1
1075 do e = 1, n
1076 d(i, 1, k, e) = d(i, 1, k, e) &
1077 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1078 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1079 d(i, ly, k, e) = d(i, ly, k, e) &
1080 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1081 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1082 end do
1083 end do
1084 end do
1085
1086 do i = 1, lx, lx-1
1087 do j = 1, ly, ly-1
1088 do e = 1, n
1089 d(i, j, 1, e) = d(i, j, 1, e) &
1090 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1091 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1092 d(i, j, lz, e) = d(i, j, lz, e) &
1093 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1094 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1095 end do
1096 end do
1097 end do
1098
1099
1100 end subroutine sx_update_lx4
1101
1102 subroutine sx_update_lx3(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
1103 integer, parameter :: lx = 3
1104 integer, parameter :: ly = 3
1105 integer, parameter :: lz = 3
1106 integer, intent(in) :: n
1107 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1108 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1109 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1110 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1111 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1112 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1113 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1114 real(kind=rp), intent(in) :: dxt(lx, lx)
1115 real(kind=rp), intent(in) :: dyt(lx, lx)
1116 real(kind=rp), intent(in) :: dzt(lx, lx)
1117 integer :: i, j, k, l, e
1118
1119
1120 do l = 1, lx
1121 do k = 1, lz
1122 do j = 1, ly
1123 do i = 1, lx
1124 do e = 1, n
1125 d(i, j, k, e) = d(i, j, k, e) + &
1126 g11(l, j, k, e) * dxt(i, l)**2
1127
1128 d(i, j, k, e) = d(i, j, k, e) + &
1129 g22(i, l, k, e) * dyt(j, l)**2
1130
1131 d(i, j, k, e) = d(i, j, k, e) + &
1132 g33(i, j, l, e) * dzt(k, l)**2
1133
1134 end do
1135 end do
1136 end do
1137 end do
1138 end do
1139
1140 do j = 1, ly, ly-1
1141 do k = 1, lz, lz-1
1142 do e = 1, n
1143 d(1, j, k, e) = d(1, j, k, e) &
1144 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1145 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1146 d(lx, j, k, e) = d(lx, j, k, e) &
1147 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1148 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1149 end do
1150 end do
1151 end do
1152
1153 do i = 1, lx, lx-1
1154 do k = 1, lz, lz-1
1155 do e = 1, n
1156 d(i, 1, k, e) = d(i, 1, k, e) &
1157 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1158 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1159 d(i, ly, k, e) = d(i, ly, k, e) &
1160 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1161 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1162 end do
1163 end do
1164 end do
1165
1166 do i = 1, lx, lx-1
1167 do j = 1, ly, ly-1
1168 do e = 1, n
1169 d(i, j, 1, e) = d(i, j, 1, e) &
1170 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1171 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1172 d(i, j, lz, e) = d(i, j, lz, e) &
1173 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1174 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1175 end do
1176 end do
1177 end do
1178
1179 end subroutine sx_update_lx3
1180
1181 subroutine sx_update_lx2(d, dxt, dyt, dzt, G11, G22, G33, G12, G13, G23, n)
1182 integer, parameter :: lx = 2
1183 integer, parameter :: ly = 2
1184 integer, parameter :: lz = 2
1185 integer, intent(in) :: n
1186 real(kind=rp), intent(inout) :: d(lx, lx, lx, n)
1187 real(kind=rp), intent(in) :: g11(lx, lx, lx, n)
1188 real(kind=rp), intent(in) :: g22(lx, lx, lx, n)
1189 real(kind=rp), intent(in) :: g33(lx, lx, lx, n)
1190 real(kind=rp), intent(in) :: g12(lx, lx, lx, n)
1191 real(kind=rp), intent(in) :: g13(lx, lx, lx, n)
1192 real(kind=rp), intent(in) :: g23(lx, lx, lx, n)
1193 real(kind=rp), intent(in) :: dxt(lx, lx)
1194 real(kind=rp), intent(in) :: dyt(lx, lx)
1195 real(kind=rp), intent(in) :: dzt(lx, lx)
1196 integer :: i, j, k, l, e
1197
1198
1199 do l = 1, lx
1200 do k = 1, lz
1201 do j = 1, ly
1202 do i = 1, lx
1203 do e = 1, n
1204 d(i, j, k, e) = d(i, j, k, e) + &
1205 g11(l, j, k, e) * dxt(i, l)**2
1206
1207 d(i, j, k, e) = d(i, j, k, e) + &
1208 g22(i, l, k, e) * dyt(j, l)**2
1209
1210 d(i, j, k, e) = d(i, j, k, e) + &
1211 g33(i, j, l, e) * dzt(k, l)**2
1212
1213 end do
1214 end do
1215 end do
1216 end do
1217 end do
1218
1219 do j = 1, ly, ly-1
1220 do k = 1, lz, lz-1
1221 do e = 1, n
1222 d(1, j, k, e) = d(1, j, k, e) &
1223 + g12(1, j, k, e) * dxt(1, 1)*dyt(j, j) &
1224 + g13(1, j, k, e) * dxt(1, 1)*dzt(k, k)
1225 d(lx, j, k, e) = d(lx, j, k, e) &
1226 + g12(lx, j, k, e) * dxt(lx, lx)*dyt(j, j) &
1227 + g13(lx, j, k, e) * dxt(lx, lx)*dzt(k, k)
1228 end do
1229 end do
1230 end do
1231
1232 do i = 1, lx, lx-1
1233 do k = 1, lz, lz-1
1234 do e = 1, n
1235 d(i, 1, k, e) = d(i, 1, k, e) &
1236 + g12(i, 1, k, e) * dyt(1, 1)*dxt(i, i) &
1237 + g23(i, 1, k, e) * dyt(1, 1)*dzt(k, k)
1238 d(i, ly, k, e) = d(i, ly, k, e) &
1239 + g12(i, ly, k, e) * dyt(ly, ly)*dxt(i, i) &
1240 + g23(i, ly, k, e) * dyt(ly, ly)*dzt(k, k)
1241 end do
1242 end do
1243 end do
1244
1245 do i = 1, lx, lx-1
1246 do j = 1, ly, ly-1
1247 do e = 1, n
1248 d(i, j, 1, e) = d(i, j, 1, e) &
1249 + g13(i, j, 1, e) * dzt(1, 1)*dxt(i, i) &
1250 + g23(i, j, 1, e) * dzt(1, 1)*dyt(j, j)
1251 d(i, j, lz, e) = d(i, j, lz, e) &
1252 + g13(i, j, lz, e) * dzt(lz, lz)*dxt(i, i) &
1253 + g23(i, j, lz, e) * dzt(lz, lz)*dyt(j, j)
1254 end do
1255 end do
1256 end do
1257
1258 end subroutine sx_update_lx2
1259
1260end module sx_jacobi
1261
coefs
Coefficients.
Definition coef.f90:34
dofmap
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
gather_scatter
Gather-scatter.
Definition gather_scatter.f90:34
math
Definition math.f90:60
math::addcol3
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:800
math::invcol1
subroutine, public invcol1(a, n)
Invert a vector .
Definition math.f90:474
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:728
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:741
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
precon
Krylov preconditioner.
Definition precon.f90:34
sx_jacobi
Jacobi preconditioner SX-Aurora backend.
Definition pc_jacobi_sx.f90:34
sx_jacobi::sx_update_lx3
subroutine sx_update_lx3(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:1103
sx_jacobi::sx_jacobi_free
subroutine sx_jacobi_free(this)
Definition pc_jacobi_sx.f90:75
sx_jacobi::sx_jacobi_solve
subroutine sx_jacobi_solve(this, z, r, n)
The jacobi preconditioner where .
Definition pc_jacobi_sx.f90:87
sx_jacobi::sx_update_lx5
subroutine sx_update_lx5(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:945
sx_jacobi::sx_update_lx12
subroutine sx_update_lx12(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:392
sx_jacobi::sx_update_lx11
subroutine sx_update_lx11(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:471
sx_jacobi::sx_update_lx2
subroutine sx_update_lx2(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:1182
sx_jacobi::sx_update_lx7
subroutine sx_update_lx7(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:787
sx_jacobi::sx_jacobi_init
subroutine sx_jacobi_init(this, coef, dof, gs_h)
Definition pc_jacobi_sx.f90:60
sx_jacobi::sx_update_lx13
subroutine sx_update_lx13(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:313
sx_jacobi::sx_update_lx4
subroutine sx_update_lx4(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:1024
sx_jacobi::sx_update_lx6
subroutine sx_update_lx6(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:866
sx_jacobi::sx_update_lx
subroutine sx_update_lx(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n, lx)
Definition pc_jacobi_sx.f90:158
sx_jacobi::sx_jacobi_update
subroutine sx_jacobi_update(this)
Definition pc_jacobi_sx.f90:95
sx_jacobi::sx_update_lx9
subroutine sx_update_lx9(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:629
sx_jacobi::sx_update_lx14
subroutine sx_update_lx14(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:234
sx_jacobi::sx_update_lx8
subroutine sx_update_lx8(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:708
sx_jacobi::sx_update_lx10
subroutine sx_update_lx10(d, dxt, dyt, dzt, g11, g22, g33, g12, g13, g23, n)
Definition pc_jacobi_sx.f90:550
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
dofmap::dofmap_t
Definition dofmap.f90:53
gather_scatter::gs_t
Definition gather_scatter.f90:58
precon::pc_t
Defines a canonical Krylov preconditioner.
Definition precon.f90:40
sx_jacobi::sx_jacobi_t
Defines a jacobi preconditioner for SX-Aurora.
Definition pc_jacobi_sx.f90:45