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