DISPATCH
stagger_2nd.f90
1 !===============================================================================
2 !> $Id: 5d1a00ac7f333dc239ec9c7ca3826d6f0b4a7604 $
3 !===============================================================================
4 MODULE stagger_mod
5  USE io_mod
6  USE trace_mod
7  USE stagger_16
8  USE stagger_20
9  USE stagger_24
10  USE stagger_32
11  USE stagger_36
12  implicit none
13  public
14  integer, private:: verbose=0
15  logical, save:: hardwire=.false.
16  type, public:: stagger_t
17  real(8):: ds(3)
18  integer(kind=8):: flops=0
19  integer(kind=8):: count=0
20  contains
21  procedure, nopass:: ddxdn1
22  procedure, nopass:: ddydn1
23  procedure, nopass:: ddzdn1
24  procedure, nopass:: ddxup1
25  procedure, nopass:: ddyup1
26  procedure, nopass:: ddzup1
27  procedure, nopass:: xdn1
28  procedure, nopass:: ydn1
29  procedure, nopass:: zdn1
30  procedure, nopass:: xup1
31  procedure, nopass:: yup1
32  procedure, nopass:: zup1
33  procedure, nopass:: ddx
34  procedure, nopass:: ddy
35  procedure, nopass:: ddz
36  procedure, nopass:: ddxdn
37  procedure, nopass:: ddydn
38  procedure, nopass:: ddzdn
39  procedure, nopass:: ddxup
40  procedure, nopass:: ddyup
41  procedure, nopass:: ddzup
42  procedure, nopass:: xdn
43  procedure, nopass:: ydn
44  procedure, nopass:: zdn
45  procedure, nopass:: xup
46  procedure, nopass:: yup
47  procedure, nopass:: zup
48  procedure, nopass:: xsm
49  procedure, nopass:: ysm
50  procedure, nopass:: zsm
51  procedure, nopass:: xyzsm
52  procedure, nopass:: stagger_test
53  end type
54  type(stagger_t), public:: stagger
55 CONTAINS
56 
57 !===============================================================================
58 FUNCTION ddxdn1 (ds, a) RESULT (b)
59  real, dimension(:,:,:), intent(in):: a
60  real, dimension(size(a,1),size(a,2),size(a,3)):: b
61  real(8):: ds(3)
62  b= ddxdn(ds, a)
63 END FUNCTION ddxdn1
64 
65 !*******************************************************************************
66 FUNCTION ddydn1 (ds, a) RESULT (b)
67  real, dimension(:,:,:), intent(in):: a
68  real, dimension(size(a,1),size(a,2),size(a,3)):: b
69  real(8):: ds(3)
70  b= ddydn(ds, a)
71 END FUNCTION ddydn1
72 
73 !*******************************************************************************
74 FUNCTION ddzdn1 (ds, a) RESULT (b)
75  real, dimension(:,:,:), intent(in):: a
76  real, dimension(size(a,1),size(a,2),size(a,3)):: b
77  real(8):: ds(3)
78  b= ddzdn(ds, a)
79 END FUNCTION ddzdn1
80 
81 !*******************************************************************************
82 FUNCTION ddxup1 (ds, a) RESULT (b)
83  real, dimension(:,:,:), intent(in):: a
84  real, dimension(size(a,1),size(a,2),size(a,3)):: b
85  real(8):: ds(3)
86  b= ddxup(ds, a)
87 END FUNCTION ddxup1
88 
89 !*******************************************************************************
90 FUNCTION ddyup1 (ds, a) RESULT (b)
91  real, dimension(:,:,:), intent(in):: a
92  real, dimension(size(a,1),size(a,2),size(a,3)):: b
93  real(8):: ds(3)
94  b= ddyup(ds, a)
95 END FUNCTION ddyup1
96 
97 !*******************************************************************************
98 FUNCTION ddzup1 (ds, a) RESULT (b)
99  real, dimension(:,:,:), intent(in):: a
100  real, dimension(size(a,1),size(a,2),size(a,3)):: b
101  real(8):: ds(3)
102  b= ddzup(ds, a)
103 END FUNCTION ddzup1
104 
105 !*******************************************************************************
106 FUNCTION xdn1 (a) RESULT (b)
107  real, dimension(:,:,:), intent(in):: a
108  real, dimension(size(a,1),size(a,2),size(a,3)):: b
109  b = xdn(a)
110 END FUNCTION xdn1
111 
112 !*******************************************************************************
113 FUNCTION ydn1 (a) RESULT (b)
114  real, dimension(:,:,:), intent(in):: a
115  real, dimension(size(a,1),size(a,2),size(a,3)):: b
116  b = ydn(a)
117 END FUNCTION ydn1
118 
119 !*******************************************************************************
120 FUNCTION zdn1 (a) RESULT (b)
121  real, dimension(:,:,:), intent(in):: a
122  real, dimension(size(a,1),size(a,2),size(a,3)):: b
123  b = zdn(a)
124 END FUNCTION zdn1
125 
126 !*******************************************************************************
127 FUNCTION xup1 (a) RESULT (b)
128  real, dimension(:,:,:), intent(in):: a
129  real, dimension(size(a,1),size(a,2),size(a,3)):: b
130  b = xup(a)
131 END FUNCTION xup1
132 
133 !*******************************************************************************
134 FUNCTION yup1 (a) RESULT (b)
135  real, dimension(:,:,:), intent(in):: a
136  real, dimension(size(a,1),size(a,2),size(a,3)):: b
137  b = yup(a)
138 END FUNCTION yup1
139 
140 !*******************************************************************************
141 FUNCTION zup1 (a) RESULT (b)
142  real, dimension(:,:,:), intent(in):: a
143  real, dimension(size(a,1),size(a,2),size(a,3)):: b
144  b = zup(a)
145 END FUNCTION zup1
146 
147 !*******************************************************************************
148 FUNCTION ddx (ds, a) RESULT (b)
149  real, dimension(:,:,:), intent(in):: a
150  real, dimension(size(a,1),size(a,2),size(a,3)):: b
151  real(8):: ds, c
152  integer :: ix, iy, iz, n(3)
153  integer, save:: itimer=0
154 !...............................................................................
155  !call trace%begin ('stagger_mod::ddx', itimer=itimer)
156  n = shape(a)
157 ! if (hardwire) then
158 ! if (all(n==[16,16,16])) then
159 ! call dxdn_16 (ds,a,b)
160 ! else if (all(n==[20,20,20])) then
161 ! call dxdn_20 (ds,a,b)
162 ! else if (all(n==[24,24,24])) then
163 ! call dxdn_24 (ds,a,b)
164 ! else if (all(n==[32,32,32])) then
165 ! call dxdn_32 (ds,a,b)
166 ! else if (all(n==[36,36,36])) then
167 ! call dxdn_36 (ds,a,b)
168 ! end if
169 ! else if (n(1) > 1) then
170  if (n(1) > 1) then
171  b( 1,:,:) = 0.0
172  b(n(1),:,:) = 0.0
173  c = 0.5/ds
174  do iz=1,n(3)
175  do iy=1,n(2)
176  !$vector always assert
177  do ix=2,n(1)-1
178  b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix-1,iy,iz))
179  end do
180  end do
181  end do
182  else
183  b = 0.0
184  endif
185  !call trace%end (itimer)
186  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
187  stagger%count = stagger%count + 1
188 END FUNCTION ddx
189 
190 !*******************************************************************************
191 FUNCTION ddy (ds, a) RESULT (b)
192  real, dimension(:,:,:), intent(in):: a
193  real, dimension(size(a,1),size(a,2),size(a,3)):: b
194  real(8):: ds, c
195  integer :: ix, iy, iz, n(3)
196  integer, save:: itimer=0
197 !...............................................................................
198  !call trace%begin ('stagger_mod::ddy', itimer=itimer)
199  n = shape(a)
200 ! if (hardwire) then
201 ! if (all(n==[16,16,16])) then
202 ! call dxdn_16 (ds,a,b)
203 ! else if (all(n==[20,20,20])) then
204 ! call dxdn_20 (ds,a,b)
205 ! else if (all(n==[24,24,24])) then
206 ! call dxdn_24 (ds,a,b)
207 ! else if (all(n==[32,32,32])) then
208 ! call dxdn_32 (ds,a,b)
209 ! else if (all(n==[36,36,36])) then
210 ! call dxdn_36 (ds,a,b)
211 ! end if
212 ! else if (n(1) > 1) then
213  if (n(2) > 1) then
214  b(:, 1,:) = 0.0
215  b(:,n(2),:) = 0.0
216  c = 0.5/ds
217  do iz=1,n(3)
218  do iy=2,n(2)-1
219  !$vector always assert
220  do ix=1,n(1)
221  b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy-1,iz))
222  end do
223  end do
224  end do
225  else
226  b = 0.0
227  endif
228  !call trace%end (itimer)
229  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
230  stagger%count = stagger%count + 1
231 END FUNCTION ddy
232 
233 !*******************************************************************************
234 FUNCTION ddz (ds, a) RESULT (b)
235  real, dimension(:,:,:), intent(in):: a
236  real, dimension(size(a,1),size(a,2),size(a,3)):: b
237  real(8):: ds, c
238  integer :: ix, iy, iz, n(3)
239  integer, save:: itimer=0
240 !...............................................................................
241  !call trace%begin ('stagger_mod::ddz', itimer=itimer)
242  n = shape(a)
243 ! if (hardwire) then
244 ! if (all(n==[16,16,16])) then
245 ! call dxdn_16 (ds,a,b)
246 ! else if (all(n==[20,20,20])) then
247 ! call dxdn_20 (ds,a,b)
248 ! else if (all(n==[24,24,24])) then
249 ! call dxdn_24 (ds,a,b)
250 ! else if (all(n==[32,32,32])) then
251 ! call dxdn_32 (ds,a,b)
252 ! else if (all(n==[36,36,36])) then
253 ! call dxdn_36 (ds,a,b)
254 ! end if
255 ! else if (n(1) > 1) then
256  if (n(3) > 1) then
257  b(:,:, 1) = 0.0
258  b(:,:,n(3)) = 0.0
259  c = 0.5/ds
260  do iz=2,n(3)-1
261  do iy=1,n(2)
262  !$vector always assert
263  do ix=1,n(1)
264  b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz-1))
265  end do
266  end do
267  end do
268  else
269  b = 0.0
270  endif
271  !call trace%end (itimer)
272  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
273  stagger%count = stagger%count + 1
274 END FUNCTION ddz
275 
276 !*******************************************************************************
277 FUNCTION ddxdn (ds, a) RESULT (b)
278  real, dimension(:,:,:), intent(in):: a
279  real, dimension(size(a,1),size(a,2),size(a,3)):: b
280  real(8):: ds(3), c
281  integer :: ix, iy, iz, n(3)
282  integer, save:: itimer=0
283 !...............................................................................
284 ! call trace%begin ('stagger_mod::ddxdn', itimer=itimer)
285  n = shape(a)
286 ! if (hardwire) then
287 ! if (all(n==[16,16,16])) then
288 ! call dxdn_16 (ds,a,b)
289 ! else if (all(n==[20,20,20])) then
290 ! call dxdn_20 (ds,a,b)
291 ! else if (all(n==[24,24,24])) then
292 ! call dxdn_24 (ds,a,b)
293 ! else if (all(n==[32,32,32])) then
294 ! call dxdn_32 (ds,a,b)
295 ! else if (all(n==[36,36,36])) then
296 ! call dxdn_36 (ds,a,b)
297 ! end if
298 ! else if (n(1) > 1) then
299  if (n(1) > 1) then
300  b(1,:,:) = 0.0
301  c = 1./ds(1)
302  do iz=1,n(3)
303  do iy=1,n(2)
304  !$vector always assert
305  do ix=2,n(1)
306  b(ix,iy,iz) = c*(a(ix ,iy,iz) - a(ix-1,iy,iz))
307  end do
308  end do
309  end do
310  else
311  b = 0.0
312  endif
313 ! call trace%end (itimer)
314  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
315  stagger%count = stagger%count + 1
316 END FUNCTION ddxdn
317 
318 !*******************************************************************************
319 FUNCTION ddydn (ds, a) RESULT (b)
320  real, dimension(:,:,:), intent(in):: a
321  real, dimension(size(a,1),size(a,2),size(a,3)):: b
322  integer :: ix, iy, iz, n(3)
323  real(8):: ds(3), c
324  integer, save:: itimer=0
325 !...............................................................................
326 ! call trace%begin ('stagger_mod::ddydn', itimer=itimer)
327  n = shape(a)
328 ! if (hardwire) then
329 ! if (all(n==[16,16,16])) then
330 ! call dydn_16 (ds,a,b)
331 ! else if (all(n==[20,20,20])) then
332 ! call dydn_20 (ds,a,b)
333 ! else if (all(n==[24,24,24])) then
334 ! call dydn_24 (ds,a,b)
335 ! else if (all(n==[32,32,32])) then
336 ! call dydn_32 (ds,a,b)
337 ! else if (all(n==[36,36,36])) then
338 ! call dydn_36 (ds,a,b)
339 ! end if
340 ! else if (n(2) > 1) then
341  if (n(2) > 1) then
342  b(:,1,:) = 0.0
343  c = 1./ds(2)
344  do iz=1,n(3)
345  do iy=2,n(2)
346  !$vector always assert
347  do ix=1,n(1)
348  b(ix,iy,iz) = c*(a(ix,iy ,iz) - a(ix,iy-1,iz))
349  end do
350  end do
351  end do
352  else
353  b = 0.0
354  endif
355 ! call trace%end (itimer)
356  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
357  stagger%count = stagger%count + 1
358 END FUNCTION ddydn
359 
360 !*******************************************************************************
361 FUNCTION ddzdn (ds, a) RESULT (b)
362  real, dimension(:,:,:), intent(in):: a
363  real, dimension(size(a,1),size(a,2),size(a,3)):: b
364  real(8):: ds(3), c
365  integer :: ix, iy, iz, n(3)
366  integer, save:: itimer=0
367 !...............................................................................
368 ! call trace%begin ('stagger_mod::ddzdn', itimer=itimer)
369  n = shape(a)
370 ! if (hardwire) then
371 ! if (all(n==[16,16,16])) then
372 ! call dzdn_16 (ds,a,b)
373 ! else if (all(n==[20,20,20])) then
374 ! call dzdn_20 (ds,a,b)
375 ! else if (all(n==[24,24,24])) then
376 ! call dzdn_24 (ds,a,b)
377 ! else if (all(n==[32,32,32])) then
378 ! call dzdn_32 (ds,a,b)
379 ! else if (all(n==[36,36,36])) then
380 ! call dzdn_36 (ds,a,b)
381 ! end if
382 ! else if (n(3) > 1) then
383  if (n(3) > 1) then
384  b(:,:,1) = 0.0
385  c = 1./ds(3)
386  do iz=2,n(3)
387  do iy=1,n(2)
388  !$vector always assert
389  do ix=1,n(1)
390  b(ix,iy,iz) = c*(a(ix,iy,iz ) - a(ix,iy,iz-1))
391  end do
392  end do
393  end do
394  else
395  b = 0.0
396  endif
397 ! call trace%end (itimer)
398  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
399  stagger%count = stagger%count + 1
400 END FUNCTION ddzdn
401 
402 !*******************************************************************************
403 FUNCTION ddxup (ds, a) RESULT (b)
404  real, dimension(:,:,:), intent(in):: a
405  real, dimension(size(a,1),size(a,2),size(a,3)):: b
406  integer :: ix, iy, iz, n(3)
407  real(8):: ds(3), c
408  integer, save:: itimer=0
409 !...............................................................................
410 ! call trace%begin ('stagger_mod::ddxup', itimer=itimer)
411  n = shape(a)
412 ! if (hardwire) then
413 ! if (all(n==[16,16,16])) then
414 ! call dxup_16 (ds,a,b)
415 ! else if (all(n==[20,20,20])) then
416 ! call dxup_20 (ds,a,b)
417 ! else if (all(n==[24,24,24])) then
418 ! call dxup_24 (ds,a,b)
419 ! else if (all(n==[32,32,32])) then
420 ! call dxup_32 (ds,a,b)
421 ! else if (all(n==[36,36,36])) then
422 ! call dxup_36 (ds,a,b)
423 ! end if
424 ! else if (n(1) > 1) then
425  if (n(1) > 1) then
426  b(n(1),:,:) = 0.0
427  c = 1./ds(1)
428  !print*,'mk'
429  do iz=1,n(3)
430  do iy=1,n(2)
431  !$vector always assert
432  do ix=1,n(1)-1
433  b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix ,iy,iz))
434  end do
435  end do
436  end do
437  else
438  b = 0.0
439  endif
440 ! call trace%end (itimer)
441  stagger%flops = stagger%flops + 2*n(3)*n(2)*(n(1)-1)
442  stagger%count = stagger%count + 1
443 END FUNCTION ddxup
444 
445 !*******************************************************************************
446 FUNCTION ddyup (ds, a) RESULT (b)
447  real, dimension(:,:,:), intent(in):: a
448  real, dimension(size(a,1),size(a,2),size(a,3)):: b
449  integer :: ix, iy, iz, n(3)
450  real(8):: ds(3), c
451  integer, save:: itimer=0
452 !...............................................................................
453 ! call trace%begin ('stagger_mod::ddyup', itimer=itimer)
454  n = shape(a)
455 ! if (hardwire) then
456 ! if (all(n==[16,16,16])) then
457 ! call dyup_16 (ds,a,b)
458 ! else if (all(n==[20,20,20])) then
459 ! call dyup_20 (ds,a,b)
460 ! else if (all(n==[24,24,24])) then
461 ! call dyup_24 (ds,a,b)
462 ! else if (all(n==[32,32,32])) then
463 ! call dyup_32 (ds,a,b)
464 ! else if (all(n==[36,36,36])) then
465 ! call dyup_36 (ds,a,b)
466 ! end if
467 ! else if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
468  if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
469  b(:,n(2),:) = 0.0
470  c = 1./ds(2)
471  do iz=1,n(3)
472  do iy=1,n(2)-1
473  !$vector always assert
474  do ix=1,n(1)
475  b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy ,iz))
476  end do
477  end do
478  end do
479  else
480  b = 0.0
481  endif
482 ! call trace%end (itimer)
483  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
484  stagger%count = stagger%count + 1
485 END FUNCTION ddyup
486 
487 !*******************************************************************************
488 FUNCTION ddzup (ds, a) RESULT (b)
489  real, dimension(:,:,:), intent(in):: a
490  real, dimension(size(a,1),size(a,2),size(a,3)):: b
491  integer :: ix, iy, iz, n(3)
492  real(8):: ds(3), c
493  integer, save:: itimer=0
494 !...............................................................................
495 ! call trace%begin ('stagger_mod::ddzup', itimer=itimer)
496  n = shape(a)
497 ! if (hardwire) then
498 ! if (all(n==[16,16,16])) then
499 ! call dzup_16 (ds,a,b)
500 ! else if (all(n==[20,20,20])) then
501 ! call dzup_20 (ds,a,b)
502 ! else if (all(n==[24,24,24])) then
503 ! call dzup_24 (ds,a,b)
504 ! else if (all(n==[32,32,32])) then
505 ! call dzup_32 (ds,a,b)
506 ! else if (all(n==[36,36,36])) then
507 ! call dzup_36 (ds,a,b)
508 ! end if
509 ! else if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
510  if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
511  b(:,:,n(3)) = 0.0
512  c = 1./ds(3)
513  do iz=1,n(3)-1
514  do iy=1,n(2)
515  !$vector always assert
516  do ix=1,n(1)
517  b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz ))
518  end do
519  end do
520  end do
521  else
522  b = 0.0
523  endif
524 ! call trace%end (itimer)
525  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
526  stagger%count = stagger%count + 1
527 END FUNCTION ddzup
528 
529 !*******************************************************************************
530 FUNCTION xdn (a) RESULT (b)
531  real, dimension(:,:,:), intent(in):: a
532  real, dimension(size(a,1),size(a,2),size(a,3)):: b
533  real :: c
534  integer :: ix, iy, iz, n(3)
535  integer, save:: itimer=0
536 !...............................................................................
537 ! call trace%begin ('stagger_mod::xdn', itimer=itimer)
538  n = shape(a)
539  if (hardwire) then
540  if (all(n==[16,16,16])) then
541  call xdn_16 (a, b)
542  else if (all(n==[20,20,20])) then
543  call xdn_20 (a, b)
544  else if (all(n==[24,24,24])) then
545  call xdn_24 (a, b)
546  else if (all(n==[32,32,32])) then
547  call xdn_32 (a, b)
548  else if (all(n==[36,36,36])) then
549  call xdn_36 (a, b)
550  end if
551  else if (n(1) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
552  c = 0.5
553  b(1,:,:) = a(1,:,:)
554  do iz=1,n(3)
555  do iy=1,n(2)
556  !$vector always assert
557  do ix=2,n(1)
558  b(ix,iy,iz) = c*(a(ix ,iy,iz) + a(ix-1,iy,iz))
559  end do
560  end do
561  end do
562  else
563  b = a
564  end if
565 ! call trace%end (itimer)
566  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
567  stagger%count = stagger%count + 1
568 END FUNCTION xdn
569 
570 !*******************************************************************************
571 FUNCTION ydn (a) RESULT (b)
572  real, dimension(:,:,:), intent(in):: a
573  real, dimension(size(a,1),size(a,2),size(a,3)):: b
574  real :: c
575  integer :: ix, iy, iz, n(3)
576  integer, save:: itimer=0
577 !...............................................................................
578 ! call trace%begin ('stagger_mod::ydn', itimer=itimer)
579  n = shape(a)
580  if (hardwire) then
581  if (all(n==[16,16,16])) then
582  call ydn_16 (a, b)
583  else if (all(n==[20,20,20])) then
584  call ydn_20 (a, b)
585  else if (all(n==[24,24,24])) then
586  call ydn_24 (a, b)
587  else if (all(n==[32,32,32])) then
588  call ydn_32 (a, b)
589  else if (all(n==[36,36,36])) then
590  call ydn_36 (a, b)
591  end if
592  else if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
593  c = 0.5
594  b(:,1,:) = a(:,1,:)
595  do iz=1,n(3)
596  do iy=2,n(2)
597  !$vector always assert
598  do ix=1,n(1)
599  b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy-1,iz))
600  end do
601  end do
602  end do
603  else
604  b = a
605  end if
606 ! call trace%end (itimer)
607  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
608  stagger%count = stagger%count + 1
609 END FUNCTION ydn
610 
611 !*******************************************************************************
612 FUNCTION zdn (a) RESULT (b)
613  real, dimension(:,:,:), intent(in):: a
614  real, dimension(size(a,1),size(a,2),size(a,3)):: b
615  real :: c
616  integer :: ix, iy, iz, n(3)
617  integer, save:: itimer=0
618 !...............................................................................
619 ! call trace%begin ('stagger_mod::zdn', itimer=itimer)
620  n = shape(a)
621  if (hardwire) then
622  if (all(n==[16,16,16])) then
623  call zdn_16 (a, b)
624  else if (all(n==[20,20,20])) then
625  call zdn_20 (a, b)
626  else if (all(n==[24,24,24])) then
627  call zdn_24 (a, b)
628  else if (all(n==[32,32,32])) then
629  call zdn_32 (a, b)
630  else if (all(n==[36,36,36])) then
631  call zdn_36 (a, b)
632  end if
633  else if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
634  c = 0.5
635  b(:,:,1) = a(:,:,1)
636  do iz=2,n(3)
637  do iy=1,n(2)
638  !$vector always assert
639  do ix=1,n(1)
640  b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy,iz-1))
641  end do
642  end do
643  end do
644  else
645  b = a
646  end if
647 ! call trace%end (itimer)
648  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
649  stagger%count = stagger%count + 1
650 END FUNCTION zdn
651 
652 !*******************************************************************************
653 FUNCTION xup (a) RESULT (b)
654  real, dimension(:,:,:), intent(in):: a
655  real, dimension(size(a,1),size(a,2),size(a,3)):: b
656  real :: c
657  integer :: ix, iy, iz, n(3)
658  integer, save:: itimer=0
659 !...............................................................................
660 ! call trace%begin ('stagger_mod::xup', itimer=itimer)
661  n = shape(a)
662  if (hardwire) then
663  if (all(n==[16,16,16])) then
664  call xup_16 (a, b)
665  else if (all(n==[20,20,20])) then
666  call xup_20 (a, b)
667  else if (all(n==[24,24,24])) then
668  call xup_24 (a, b)
669  else if (all(n==[32,32,32])) then
670  call xup_32 (a, b)
671  else if (all(n==[36,36,36])) then
672  call xup_36 (a, b)
673  end if
674  else if (n(1) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
675  c = 0.5
676  b(n(1),:,:) = a(n(1),:,:)
677  do iz=1,n(3)
678  do iy=1,n(2)
679  !$vector always assert
680  do ix=1,n(1)-1
681  b(ix,iy,iz) = c*(a(ix+1,iy,iz) + a(ix ,iy,iz))
682  end do
683  end do
684  end do
685  else
686  b = a
687  endif
688 ! call trace%end (itimer)
689  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
690  stagger%count = stagger%count + 1
691 END FUNCTION
692 
693 !*******************************************************************************
694 FUNCTION yup (a) RESULT (b)
695  real, dimension(:,:,:), intent(in):: a
696  real, dimension(size(a,1),size(a,2),size(a,3)):: b
697  real :: c
698  integer :: ix, iy, iz, n(3)
699  integer, save:: itimer=0
700 !...............................................................................
701 ! call trace%begin ('stagger_mod::yup', itimer=itimer)
702  n = shape(a)
703  if (hardwire) then
704  if (all(n==[16,16,16])) then
705  call yup_16 (a, b)
706  else if (all(n==[20,20,20])) then
707  call yup_20 (a, b)
708  else if (all(n==[24,24,24])) then
709  call yup_24 (a, b)
710  else if (all(n==[32,32,32])) then
711  call yup_32 (a, b)
712  else if (all(n==[36,36,36])) then
713  call yup_36 (a, b)
714  end if
715  else if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
716  c = 0.5
717  b(:,n(2),:) = a(:,n(2),:)
718  do iz=1,n(3)
719  do iy=1,n(2)-1
720  !$vector always assert
721  do ix=1,n(1)
722  b(ix,iy,iz) = c*(a(ix,iy+1,iz) + a(ix,iy ,iz))
723  end do
724  end do
725  end do
726  else
727  b = a
728  endif
729 ! call trace%end (itimer)
730  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
731  stagger%count = stagger%count + 1
732 END FUNCTION
733 
734 !*******************************************************************************
735 FUNCTION zup (a) RESULT (b)
736  real, dimension(:,:,:), intent(in):: a
737  real, dimension(size(a,1),size(a,2),size(a,3)):: b
738  real :: c
739  integer :: ix, iy, iz, n(3)
740  integer, save:: itimer=0
741 !...............................................................................
742 ! call trace%begin ('stagger_mod::zup', itimer=itimer)
743  n = shape(a)
744  if (hardwire) then
745  if (all(n==[16,16,16])) then
746  call zup_16 (a, b)
747  else if (all(n==[20,20,20])) then
748  call zup_20 (a, b)
749  else if (all(n==[24,24,24])) then
750  call zup_24 (a, b)
751  else if (all(n==[32,32,32])) then
752  call zup_32 (a, b)
753  else if (all(n==[36,36,36])) then
754  call zup_36 (a, b)
755  end if
756  else if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
757  c = 0.5
758  b(:,:,n(3)) = a(:,:,n(3))
759  do iz=1,n(3)-1
760  do iy=1,n(2)
761  !$vector always assert
762  do ix=1,n(1)
763  b(ix,iy,iz) = c*(a(ix,iy,iz+1) + a(ix,iy,iz ))
764  end do
765  end do
766  end do
767  else
768  b = a
769  endif
770 ! call trace%end (itimer)
771  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
772  stagger%count = stagger%count + 1
773 END FUNCTION
774 
775 !*******************************************************************************
776 FUNCTION xsm (a) RESULT (b)
777  real, dimension(:,:,:), intent(in):: a
778  real, dimension(size(a,1),size(a,2),size(a,3)):: b
779  real :: c
780  integer :: i, n
781  integer, save:: itimer=0
782 !...............................................................................
783 ! call trace%begin ('stagger_mod::xsm', itimer=itimer)
784  n = size(a,1)
785  if (n > 1) then
786  c = 0.25
787  do i=2,n-1
788  b(i,:,:) = (a(i-1,:,:) + a(i+1,:,:))*c + a(i,:,:)*(1.-2.*c)
789  end do
790  do i=1,n,n-1
791  b(i,:,:) = a(i,:,:)
792  end do
793  else
794  b = a
795  endif
796 ! call trace%end (itimer)
797  stagger%flops = stagger%flops + 4*n*n*(n-2)
798  stagger%count = stagger%count + 1
799 END FUNCTION
800 
801 !*******************************************************************************
802 FUNCTION ysm (a) RESULT (b)
803  real, dimension(:,:,:), intent(in):: a
804  real, dimension(size(a,1),size(a,2),size(a,3)):: b
805  real :: c
806  integer :: i, n
807  integer, save:: itimer=0
808 !...............................................................................
809 ! call trace%begin ('stagger_mod::ysm', itimer=itimer)
810  n = size(a,2)
811  if (n > 1) then
812  c = 0.25
813  do i=2,n-1
814  b(:,i,:) = (a(:,i-1,:) + a(:,i+1,:))*c + a(:,i,:)*(1.-2.*c)
815  end do
816  do i=1,n,n-1
817  b(:,i,:) = a(:,i,:)
818  end do
819  else
820  b = a
821  endif
822 ! call trace%end (itimer)
823  stagger%flops = stagger%flops + 4*n*n*(n-2)
824  stagger%count = stagger%count + 1
825 END FUNCTION
826 
827 !*******************************************************************************
828 FUNCTION zsm (a) RESULT (b)
829  real, dimension(:,:,:), intent(in):: a
830  real, dimension(size(a,1),size(a,2),size(a,3)):: b
831  real :: c
832  integer :: i, n
833  integer, save:: itimer=0
834 !...............................................................................
835 ! call trace%begin ('stagger_mod::zsm', itimer=itimer)
836  n = size(a,3)
837  if (n > 1) then
838  c = 0.25
839  do i=2,n-1
840  b(:,:,i) = (a(:,:,i-1) + a(:,:,i+1))*c + a(:,:,i)*(1.-2.*c)
841  end do
842  do i=1,n,n-1
843  b(:,:,i) = a(:,:,i)
844  end do
845  else
846  b = a
847  endif
848 ! call trace%end (itimer)
849  stagger%flops = stagger%flops + 4*n*n*(n-2)
850  stagger%count = stagger%count + 1
851 END FUNCTION
852 
853 !*******************************************************************************
854 FUNCTION xyzsm(a) RESULT (b)
855  real, dimension(:,:,:), intent(in):: a
856  real, dimension(size(a,1),size(a,2),size(a,3)):: b
857  integer, save:: itimer=0
858  call trace%begin ('stagger_mod::xyzsm', itimer=itimer)
859  b = zsm(ysm(xsm(a)))
860  call trace%end (itimer)
861 END FUNCTION
862 
863 !===============================================================================
864 !> Unit test of 6th/5th order operators
865 !===============================================================================
866 SUBROUTINE stagger_test
867  integer :: m(3)
868  real, pointer :: a(:,:,:)
869  real, pointer :: b(:,:,:)
870  real, pointer :: c(:,:,:)
871  real(8) :: ds(3), x, y, z, d
872  integer :: i
873  logical :: ok, allok
874  real, parameter :: eps = 2e-6
875  !-----------------------------------------------------------------------------
876  m = [10,11,12]
877  ds = 1d0/m
878  allocate (a(m(1),m(2),m(3)))
879  allocate (b(m(1),m(2),m(3)))
880  allok = .true.
881  !-----------------------------------------------------------------------------
882  ! x-operators
883  !-----------------------------------------------------------------------------
884  do i=1,m(1)
885  x = i*ds(1)
886  a(i,:,:) = x
887  end do
888  b = xdn(a)
889  ok = .true.
890  do i=4,m(1)-2
891  x = (i-0.5d0)*ds(1)
892  d = x
893  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
894  end do
895  allok = allok .and. ok
896  if (.not.ok) print *,'xdn:', ok
897  b = xup(a)
898  ok = .true.
899  do i=3,m(1)-3
900  x = (i+0.5d0)*ds(1)
901  d = x
902  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
903  end do
904  allok = allok .and. ok
905  if (.not.ok) print *,'xup:', ok
906  !-----------------------------------------------------------------------------
907  ! y-operators
908  !-----------------------------------------------------------------------------
909  do i=1,m(2)
910  x = i*ds(2)
911  a(:,i,:) = x
912  end do
913  b = ydn(a)
914  ok = .true.
915  do i=4,m(2)-2
916  y = (i-0.5d0)*ds(2)
917  d = y
918  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
919  end do
920  allok = allok .and. ok
921  if (.not.ok) print *,'ydn:', ok
922  b = yup(a)
923  ok = .true.
924  do i=3,m(2)-3
925  y = (i+.5d0)*ds(2)
926  d = y
927  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
928  end do
929  allok = allok .and. ok
930  if (.not.ok) print *,'yup:', ok
931  !-----------------------------------------------------------------------------
932  ! z-operators
933  !-----------------------------------------------------------------------------
934  do i=1,m(3)
935  x = i*ds(3)
936  a(:,:,i) = x
937  end do
938  b = zdn(a)
939  ok = .true.
940  do i=4,m(3)-2
941  z = (i-0.5d0)*ds(3)
942  d = z
943  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
944  end do
945  allok = allok .and. ok
946  if (.not.ok) print *,'zdn:', ok
947  b = zup(a)
948  ok = .true.
949  do i=3,m(3)-3
950  z = (i+0.5d0)*ds(3)
951  d = z
952  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
953  end do
954  allok = allok .and. ok
955  if (.not.ok) print *,'zup:', ok
956  !-----------------------------------------------------------------------------
957  ! x-operators
958  !-----------------------------------------------------------------------------
959  do i=1,m(1)
960  x = i*ds(1)
961  a(i,:,:) = 1d0 + x**2
962  end do
963  b = ddxdn(ds,a)
964  ok = .true.
965  do i=4,m(1)-2
966  x = (i-0.5d0)*ds(1)
967  d = 2d0*x
968  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
969  end do
970  allok = allok .and. ok
971  if (.not.ok) print *,'ddxdn:', ok
972  b = ddxup(ds,a)
973  ok = .true.
974  do i=3,m(1)-3
975  x = (i+0.5d0)*ds(1)
976  d = 2d0*x
977  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
978  end do
979  allok = allok .and. ok
980  if (.not.ok) print *,'ddxup:', ok
981  !-----------------------------------------------------------------------------
982  ! y-operators
983  !-----------------------------------------------------------------------------
984  do i=1,m(2)
985  x = i*ds(2)
986  a(:,i,:) = 1d0 + x**2
987  end do
988  b = ddydn(ds,a)
989  ok = .true.
990  do i=4,m(2)-2
991  y = (i-0.5d0)*ds(2)
992  d = 2d0*y
993  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
994  end do
995  allok = allok .and. ok
996  if (.not.ok) print *,'ddydn:', ok
997  b = ddyup(ds,a)
998  ok = .true.
999  do i=3,m(2)-3
1000  y = (i+.5d0)*ds(2)
1001  d = 2d0*y
1002  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1003  end do
1004  allok = allok .and. ok
1005  if (.not.ok) print *,'ddyup:', ok
1006  !-----------------------------------------------------------------------------
1007  ! z-operators
1008  !-----------------------------------------------------------------------------
1009  do i=1,m(3)
1010  x = i*ds(3)
1011  a(:,:,i) = 1d0 + x**2
1012  end do
1013  b = ddzdn(ds,a)
1014  ok = .true.
1015  do i=4,m(3)-2
1016  z = (i-0.5d0)*ds(3)
1017  d = 2d0*z
1018  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1019  end do
1020  allok = allok .and. ok
1021  if (.not.ok) print *,'ddzdn:', ok
1022  b = ddzup(ds,a)
1023  ok = .true.
1024  do i=3,m(3)-3
1025  z = (i+0.5d0)*ds(3)
1026  d = 2d0*z
1027  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1028  end do
1029  allok = allok .and. ok
1030  if (.not.ok) print *,'ddzup:', ok
1031  print *,io%hl
1032  if (allok) then
1033  print *,'2nd/1st order stagger operator test passed'
1034  print *,io%hl
1035  else
1036  print *,'2nd/1st order stagger operator test failed!'
1037  print *,io%hl
1038  error stop 'stagger_test'
1039  end if
1040  !-----------------------------------------------------------------------------
1041  deallocate (a, b)
1042 END SUBROUTINE stagger_test
1043 
1044 END MODULE
Definition: io_mod.f90:4
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4