DISPATCH
stagger_6th.f90
1 !===============================================================================
2 !> 6th order stagger operators, with self-test procedure
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 FUNCTION ddxdn (ds, a) RESULT (b)
384  real(8) :: ds(3)
385  real :: a(:,:,:)
386  real :: b(size(a,1),size(a,2),size(a,3))
387  !.............................................................................
388  integer :: i, j, k, n(3)
389  real :: c1, c2, c3, d
390  !-----------------------------------------------------------------------------
391  n = shape(a)
392  if (n(1) > 5) then
393  d = ds(1)
394  c3 = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
395  c2 = (-1.-120.*c3)/24.
396  c1 = (1.-3.*c2-5.*c3)/d
397  c2 = c2/d
398  c3 = c3/d
399  b( 1:3,:,:) = 0.0
400  b(n(1)-1:,:,:) = 0.0
401  do k=1,n(3)
402  do j=1,n(2)
403  do i=4,n(1)-2
404  b(i ,j,k) = (c3*(a(i+2,j,k)-a(i-3,j,k)) &
405  + c2*(a(i+1,j,k)-a(i-2,j,k))) &
406  + c1*(a(i ,j,k)-a(i-1,j,k))
407  end do
408  end do
409  end do
410  else
411  b = 0.0
412  endif
413  stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
414  stagger%count = stagger%count + 1
415 END FUNCTION ddxdn
416 
417 !*******************************************************************************
418 FUNCTION ddydn (ds, a) RESULT (b)
419  real(8) :: ds(3)
420  real :: a(:,:,:)
421  real :: b(size(a,1),size(a,2),size(a,3))
422  !.............................................................................
423  integer :: i, j, k, n(3)
424  real :: c1, c2, c3, d
425  !-----------------------------------------------------------------------------
426  n = shape(a)
427  if (n(2) > 5) then
428  d = ds(2)
429  c3 = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
430  c2 = (-1.-120.*c3)/24.
431  c1 = (1.-3.*c2-5.*c3)/d
432  c2 = c2/d
433  c3 = c3/d
434  b(:, 1:3,:) = 0.0
435  b(:,n(2)-1:,:) = 0.0
436  do k=1,n(3)
437  do j=4,n(2)-2
438  do i=1,n(1)
439  b(i,j,k) = (c3*(a(i,j+2,k)-a(i,j-3,k)) &
440  + c2*(a(i,j+1,k)-a(i,j-2,k))) &
441  + c1*(a(i,j ,k)-a(i,j-1,k))
442  end do
443  end do
444  end do
445  else
446  b = 0.0
447  endif
448  stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
449  stagger%count = stagger%count + 1
450 END FUNCTION ddydn
451 
452 !*******************************************************************************
453 FUNCTION ddzdn (ds, a) RESULT (b)
454  real(8) :: ds(3)
455  real :: a(:,:,:)
456  real :: b(size(a,1),size(a,2),size(a,3))
457  !.............................................................................
458  integer :: i, j, k, n(3)
459  real :: c1, c2, c3, d
460  !-----------------------------------------------------------------------------
461  n = shape(a)
462  if (n(3) > 5) then
463  d = ds(3)
464  c3 = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
465  c2 = (-1.-120.*c3)/24.
466  c1 = (1.-3.*c2-5.*c3)/d
467  c2 = c2/d
468  c3 = c3/d
469  b(:,:, 1:3) = 0.0
470  b(:,:,n(3)-1:) = 0.0
471  do k=4,n(3)-2
472  do j=1,n(2)
473  do i=1,n(1)
474  b(i,j,k) = (c3*(a(i,j,k+2)-a(i,j,k-3)) &
475  + c2*(a(i,j,k+1)-a(i,j,k-2))) &
476  + c1*(a(i,j,k )-a(i,j,k-1))
477  end do
478  end do
479  end do
480  else
481  b = 0.0
482  endif
483  stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
484  stagger%count = stagger%count + 1
485 END FUNCTION ddzdn
486 
487 !*******************************************************************************
488 FUNCTION ddxup (ds, a) RESULT (b)
489  real(8) :: ds(3)
490  real :: a(:,:,:)
491  real :: b(size(a,1),size(a,2),size(a,3))
492  !.............................................................................
493  integer :: i, j, k, n(3)
494  real :: c1, c2, c3, d
495  !-----------------------------------------------------------------------------
496  n = shape(a)
497  d = ds(1)
498  if (n(1) > 5) then
499  c3 = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
500  c2 = (-1.-120.*c3)/24.
501  c1 = (1.-3.*c2-5.*c3)/d
502  c2 = c2/d
503  c3 = c3/d
504  b( 1:2,:,:) = 0.0
505  b(n(1)-2:,:,:) = 0.0
506  do k=1,n(3)
507  do j=1,n(2)
508  do i=3,n(1)-3
509  b(i ,j,k) = (c3*(a(i+3,j,k)-a(i-2,j,k)) &
510  + c2*(a(i+2,j,k)-a(i-1,j,k))) &
511  + c1*(a(i+1,j,k)-a(i ,j,k))
512  end do
513  end do
514  end do
515  else
516  b = 0.0
517  endif
518  stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
519  stagger%count = stagger%count + 1
520 END FUNCTION ddxup
521 
522 !*******************************************************************************
523 FUNCTION ddyup (ds, a) RESULT (b)
524  real(8) :: ds(3)
525  real :: a(:,:,:)
526  real :: b(size(a,1),size(a,2),size(a,3))
527  !.............................................................................
528  integer :: i, j, k, n(3)
529  real :: c1, c2, c3, d
530  !-----------------------------------------------------------------------------
531  n = shape(a)
532  d = ds(2)
533  if (n(2) > 5) then
534  c3 = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
535  c2 = (-1.-120.*c3)/24.
536  c1 = (1.-3.*c2-5.*c3)/d
537  c2 = c2/d
538  c3 = c3/d
539  b(:, 1:2,:) = 0.0
540  b(:,n(2)-2:,:) = 0.0
541  do k=1,n(3)
542  do j=3,n(2)-3
543  do i=1,n(1)
544  b(i,j,k) = (c3*(a(i,j+3,k)-a(i,j-2,k)) &
545  + c2*(a(i,j+2,k)-a(i,j-1,k))) &
546  + c1*(a(i,j+1,k)-a(i,j ,k))
547  end do
548  end do
549  end do
550  else
551  b = 0.0
552  endif
553  stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
554  stagger%count = stagger%count + 1
555 END FUNCTION ddyup
556 
557 !*******************************************************************************
558 FUNCTION ddzup (ds, a) RESULT (b)
559  real(8) :: ds(3)
560  real :: a(:,:,:)
561  real :: b(size(a,1),size(a,2),size(a,3))
562  !.............................................................................
563  integer :: i, j, k, n(3)
564  real :: c1, c2, c3, d
565  !-----------------------------------------------------------------------------
566  n = shape(a)
567  d = ds(3)
568  if (n(3) > 5) then
569  c3 = (-1.+(3.**5-3.)/(3.**3-3.))/(5.**5-5.-5.*(3.**5-3))
570  c2 = (-1.-120.*c3)/24.
571  c1 = (1.-3.*c2-5.*c3)/d
572  c2 = c2/d
573  c3 = c3/d
574  b(:,:, 1:2) = 0.0
575  b(:,:,n(1)-2:) = 0.0
576  do k=3,n(3)-3
577  do j=1,n(2)
578  do i=1,n(1)
579  b(i,j,k) = (c3*(a(i,j,k+3)-a(i,j,k-2)) &
580  + c2*(a(i,j,k+2)-a(i,j,k-1))) &
581  + c1*(a(i,j,k+1)-a(i,j,k ))
582  end do
583  end do
584  end do
585  else
586  b = 0.0
587  endif
588  stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
589  stagger%count = stagger%count + 1
590 END FUNCTION ddzup
591 
592 !*******************************************************************************
593 FUNCTION xdn (a) RESULT (b)
594  real :: a(:,:,:)
595  real :: b(size(a,1),size(a,2),size(a,3))
596  !.............................................................................
597  integer :: i, j, k, n(3)
598  real :: c1, c2, c3, d
599  !-----------------------------------------------------------------------------
600  n = shape(a)
601  if (n(1) > 5) then
602  c3 = 3./256.
603  c2 = -25./256.
604  c1 = .5-c2-c3
605  b( 1:3,:,:) = a( 1:3,:,:)
606  b(n(1)-1:,:,:) = a(n(1)-1:,:,:)
607  do k=1,n(3)
608  do j=1,n(2)
609  do i=4,n(1)-2
610  b(i ,j,k) = (c3*(a(i+2,j,k)+a(i-3,j,k)) &
611  + c2*(a(i+1,j,k)+a(i-2,j,k))) &
612  + c1*(a(i ,j,k)+a(i-1,j,k))
613  end do
614  end do
615  end do
616  else
617  b = 0.0
618  endif
619  stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
620  stagger%count = stagger%count + 1
621 END FUNCTION xdn
622 
623 !*******************************************************************************
624 FUNCTION ydn (a) RESULT (b)
625  real :: a(:,:,:)
626  real :: b(size(a,1),size(a,2),size(a,3))
627  !.............................................................................
628  integer :: i, j, k, n(3)
629  real :: c1, c2, c3, d
630  !-----------------------------------------------------------------------------
631  n = shape(a)
632  if (n(2) > 5) then
633  c3 = 3./256.
634  c2 = -25./256.
635  c1 = .5-c2-c3
636  b(:, 1:3,:) = a(:, 1:3,:)
637  b(:,n(2)-1:,:) = a(:,n(2)-1:,:)
638  do k=1,n(3)
639  do j=4,n(2)-2
640  do i=1,n(1)
641  b(i,j,k) = (c3*(a(i,j+2,k)+a(i,j-3,k)) &
642  + c2*(a(i,j+1,k)+a(i,j-2,k))) &
643  + c1*(a(i,j ,k)+a(i,j-1,k))
644  end do
645  end do
646  end do
647  else
648  b = 0.0
649  endif
650  stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
651  stagger%count = stagger%count + 1
652 END FUNCTION ydn
653 
654 !*******************************************************************************
655 FUNCTION zdn (a) RESULT (b)
656  real :: a(:,:,:)
657  real :: b(size(a,1),size(a,2),size(a,3))
658  !.............................................................................
659  integer :: i, j, k, n(3)
660  real :: c1, c2, c3, d
661  !-----------------------------------------------------------------------------
662  n = shape(a)
663  if (n(3) > 5) then
664  c3 = 3./256.
665  c2 = -25./256.
666  c1 = .5-c2-c3
667  b(:,:, 1:3) = a(:,:, 1:3)
668  b(:,:,n(1)-1:) = a(:,:,n(1)-1:)
669  do k=4,n(3)-2
670  do j=1,n(2)
671  do i=1,n(1)
672  b(i,j,k) = (c3*(a(i,j,k+2)+a(i,j,k-3)) &
673  + c2*(a(i,j,k+1)+a(i,j,k-2))) &
674  + c1*(a(i,j,k )+a(i,j,k-1))
675  end do
676  end do
677  end do
678  else
679  b = 0.0
680  endif
681  stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
682  stagger%count = stagger%count + 1
683 END FUNCTION zdn
684 
685 !*******************************************************************************
686 FUNCTION xup (a) RESULT (b)
687  real :: a(:,:,:)
688  real :: b(size(a,1),size(a,2),size(a,3))
689  !.............................................................................
690  integer :: i, j, k, n(3)
691  real :: c1, c2, c3, d
692  !-----------------------------------------------------------------------------
693  n = shape(a)
694  if (n(1) > 5) then
695  c3 = 3./256.
696  c2 = -25./256.
697  c1 = .5-c2-c3
698  b( 1:2,:,:) = a( 1:2,:,:)
699  b(n(1)-2:,:,:) = a(n(1)-2:,:,:)
700  do k=1,n(3)
701  do j=1,n(2)
702  do i=3,n(1)-3
703  b(i ,j,k) = (c3*(a(i+3,j,k)+a(i-2,j,k)) &
704  + c2*(a(i+2,j,k)+a(i-1,j,k))) &
705  + c1*(a(i+1,j,k)+a(i ,j,k))
706  end do
707  end do
708  end do
709  else
710  b = 0.0
711  endif
712  stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
713  stagger%count = stagger%count + 1
714 END FUNCTION xup
715 
716 !*******************************************************************************
717 FUNCTION yup (a) RESULT (b)
718  real :: a(:,:,:)
719  real :: b(size(a,1),size(a,2),size(a,3))
720  !.............................................................................
721  integer :: i, j, k, n(3)
722  real :: c1, c2, c3, d
723  !-----------------------------------------------------------------------------
724  n = shape(a)
725  if (n(2) > 5) then
726  c3 = 3./256.
727  c2 = -25./256.
728  c1 = .5-c2-c3
729  b(:, 1:2,:) = a(:, 1:2,:)
730  b(:,n(2)-2:,:) = a(:,n(2)-2:,:)
731  do k=1,n(3)
732  do j=3,n(2)-3
733  do i=1,n(1)
734  b(i,j,k) = (c3*(a(i,j+3,k)+a(i,j-2,k)) &
735  + c2*(a(i,j+2,k)+a(i,j-1,k))) &
736  + c1*(a(i,j+1,k)+a(i,j ,k))
737  end do
738  end do
739  end do
740  else
741  b = 0.0
742  endif
743  stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
744  stagger%count = stagger%count + 1
745 END FUNCTION yup
746 
747 !*******************************************************************************
748 FUNCTION zup (a) RESULT (b)
749  real :: a(:,:,:)
750  real :: b(size(a,1),size(a,2),size(a,3))
751  !.............................................................................
752  integer :: i, j, k, n(3)
753  real :: c1, c2, c3, d
754  !-----------------------------------------------------------------------------
755  n = shape(a)
756  if (n(3) > 5) then
757  c3 = 3./256.
758  c2 = -25./256.
759  c1 = .5-c2-c3
760  b(:,:, 1:2) = a(:,:, 1:2)
761  b(:,:,n(1)-2:) = a(:,:,n(1)-2:)
762  do k=3,n(3)-3
763  do j=1,n(2)
764  do i=1,n(1)
765  b(i,j,k) = (c3*(a(i,j,k+3)+a(i,j,k-2)) &
766  + c2*(a(i,j,k+2)+a(i,j,k-1))) &
767  + c1*(a(i,j,k+1)+a(i,j,k ))
768  end do
769  end do
770  end do
771  else
772  b = 0.0
773  endif
774  stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
775  stagger%count = stagger%count + 1
776 END FUNCTION zup
777 
778 !*******************************************************************************
779 FUNCTION xsm (a) RESULT (b)
780  real, dimension(:,:,:), intent(in):: a
781  real, dimension(size(a,1),size(a,2),size(a,3)):: b
782  real :: c
783  integer :: i, n
784  integer, save:: itimer=0
785 !...............................................................................
786 ! call trace%begin ('stagger_mod::xsm', itimer=itimer)
787  n = size(a,1)
788  if (n > 1) then
789  c = 0.25
790  do i=2,n-1
791  b(i,:,:) = (a(i-1,:,:) + a(i+1,:,:))*c + a(i,:,:)*(1.-2.*c)
792  end do
793  do i=1,n,n-1
794  b(i,:,:) = a(i,:,:)
795  end do
796  else
797  b = a
798  endif
799 ! call trace%end (itimer)
800  stagger%flops = stagger%flops + 4*n*n*(n-2)
801  stagger%count = stagger%count + 1
802 END FUNCTION
803 
804 !*******************************************************************************
805 FUNCTION ysm (a) RESULT (b)
806  real, dimension(:,:,:), intent(in):: a
807  real, dimension(size(a,1),size(a,2),size(a,3)):: b
808  real :: c
809  integer :: i, n
810  integer, save:: itimer=0
811 !...............................................................................
812 ! call trace%begin ('stagger_mod::ysm', itimer=itimer)
813  n = size(a,2)
814  if (n > 1) then
815  c = 0.25
816  do i=2,n-1
817  b(:,i,:) = (a(:,i-1,:) + a(:,i+1,:))*c + a(:,i,:)*(1.-2.*c)
818  end do
819  do i=1,n,n-1
820  b(:,i,:) = a(:,i,:)
821  end do
822  else
823  b = a
824  endif
825 ! call trace%end (itimer)
826  stagger%flops = stagger%flops + 4*n*n*(n-2)
827  stagger%count = stagger%count + 1
828 END FUNCTION
829 
830 !*******************************************************************************
831 FUNCTION zsm (a) RESULT (b)
832  real, dimension(:,:,:), intent(in):: a
833  real, dimension(size(a,1),size(a,2),size(a,3)):: b
834  real :: c
835  integer :: i, n
836  integer, save:: itimer=0
837 !...............................................................................
838 ! call trace%begin ('stagger_mod::zsm', itimer=itimer)
839  n = size(a,3)
840  if (n > 1) then
841  c = 0.25
842  do i=2,n-1
843  b(:,:,i) = (a(:,:,i-1) + a(:,:,i+1))*c + a(:,:,i)*(1.-2.*c)
844  end do
845  do i=1,n,n-1
846  b(:,:,i) = a(:,:,i)
847  end do
848  else
849  b = a
850  endif
851 ! call trace%end (itimer)
852  stagger%flops = stagger%flops + 4*n*n*(n-2)
853  stagger%count = stagger%count + 1
854 END FUNCTION
855 
856 !*******************************************************************************
857 FUNCTION xyzsm(a) RESULT (b)
858  real, dimension(:,:,:), intent(in):: a
859  real, dimension(size(a,1),size(a,2),size(a,3)):: b
860  integer, save:: itimer=0
861  call trace%begin ('stagger_mod::xyzsm', itimer=itimer)
862  b = zsm(ysm(xsm(a)))
863  call trace%end (itimer)
864 END FUNCTION
865 
866 !===============================================================================
867 !> Unit test of 6th/5th order operators
868 !===============================================================================
869 SUBROUTINE stagger_test
870  USE mpi_mod
871  ! --------------------------------------------------------------------
872  integer :: n(3)
873  real, pointer :: a(:,:,:)
874  real, pointer :: b(:,:,:)
875  real, pointer :: c(:,:,:)
876  real(8) :: ds(3), x, y, z, d
877  integer :: i
878  logical :: ok, allok
879  real, parameter :: eps = 2e-6
880  !-----------------------------------------------------------------------------
881  n = [10,11,12]
882  ds = 1d0/n
883  allocate (a(n(1),n(2),n(3)))
884  allocate (b(n(1),n(2),n(3)))
885  allok = .true.
886  !-----------------------------------------------------------------------------
887  ! x-operators
888  !-----------------------------------------------------------------------------
889  do i=1,n(1)
890  x = i*ds(1)
891  a(i,:,:) = x**2
892  end do
893  b = ddxdn1(ds,a)
894  ok = .true.
895  do i=4,n(1)-2
896  x = (i-0.5d0)*ds(1)
897  d = 2.*x
898  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
899  end do
900  allok = allok .and. ok
901  if (.not.ok) print *,'ddxdn1:', ok
902  b = ddxup1(ds,a)
903  ok = .true.
904  do i=3,n(1)-3
905  x = (i+0.5d0)*ds(1)
906  d = 2.*x
907  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
908  end do
909  allok = allok .and. ok
910  if (.not.ok) print *,'ddxup1:', ok
911  !-----------------------------------------------------------------------------
912  ! y-operators
913  !-----------------------------------------------------------------------------
914  do i=1,n(2)
915  y = i*ds(2)
916  a(:,i,:) = y**2
917  end do
918  b = ddydn1(ds,a)
919  ok = .true.
920  do i=4,n(2)-2
921  y = (i-0.5d0)*ds(2)
922  d = 2.*y
923  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
924  end do
925  allok = allok .and. ok
926  if (.not.ok) print *,'ddydn1:', ok
927  b = ddyup1(ds,a)
928  ok = .true.
929  do i=3,n(2)-3
930  y = (i+.5d0)*ds(2)
931  d = 2.*y
932  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
933  end do
934  allok = allok .and. ok
935  if (.not.ok) print *,'ddyup1:', ok
936  !-----------------------------------------------------------------------------
937  ! z-operators
938  !-----------------------------------------------------------------------------
939  do i=1,n(3)
940  z = i*ds(3)
941  a(:,:,i) = z**2
942  end do
943  b = ddzdn1(ds,a)
944  ok = .true.
945  do i=4,n(3)-2
946  z = (i-0.5d0)*ds(3)
947  d = 2.*z
948  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
949  end do
950  allok = allok .and. ok
951  if (.not.ok) print *,'ddzdn1:', ok
952  b = ddzup1(ds,a)
953  ok = .true.
954  do i=3,n(3)-3
955  z = (i+0.5d0)*ds(3)
956  d = 2.*z
957  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
958  end do
959  allok = allok .and. ok
960  if (.not.ok) print *,'ddzup1:', ok
961  !-----------------------------------------------------------------------------
962  ! x-operators
963  !-----------------------------------------------------------------------------
964  do i=1,n(1)
965  x = i*ds(1)
966  a(i,:,:) = x
967  end do
968  b = xdn1(a)
969  ok = .true.
970  do i=4,n(1)-2
971  x = (i-0.5d0)*ds(1)
972  d = x
973  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
974  end do
975  allok = allok .and. ok
976  if (.not.ok) print *,'xdn1:', ok
977  b = xup1(a)
978  ok = .true.
979  do i=3,n(1)-3
980  x = (i+0.5d0)*ds(1)
981  d = x
982  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
983  end do
984  allok = allok .and. ok
985  if (.not.ok) print *,'xup1:', ok
986  !-----------------------------------------------------------------------------
987  ! y-operators
988  !-----------------------------------------------------------------------------
989  do i=1,n(2)
990  y = i*ds(2)
991  a(:,i,:) = y
992  end do
993  b = ydn1(a)
994  ok = .true.
995  do i=4,n(2)-2
996  y = (i-0.5d0)*ds(2)
997  d = y
998  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
999  end do
1000  allok = allok .and. ok
1001  if (.not.ok) print *,'ydn1:', ok
1002  b = yup1(a)
1003  ok = .true.
1004  do i=3,n(2)-3
1005  y = (i+.5d0)*ds(2)
1006  d = y
1007  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1008  end do
1009  allok = allok .and. ok
1010  if (.not.ok) print *,'yup1:', ok
1011  !-----------------------------------------------------------------------------
1012  ! z-operators
1013  !-----------------------------------------------------------------------------
1014  do i=1,n(3)
1015  z = i*ds(3)
1016  a(:,:,i) = z
1017  end do
1018  b = zdn1(a)
1019  ok = .true.
1020  do i=4,n(3)-2
1021  z = (i-0.5d0)*ds(3)
1022  d = z
1023  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1024  end do
1025  allok = allok .and. ok
1026  if (.not.ok) print *,'zdn1:', ok
1027  b = zup1(a)
1028  ok = .true.
1029  do i=3,n(3)-3
1030  z = (i+0.5d0)*ds(3)
1031  d = z
1032  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1033  end do
1034  allok = allok .and. ok
1035  if (.not.ok) print *,'zup1:', ok
1036  !-----------------------------------------------------------------------------
1037  ! x-operators
1038  !-----------------------------------------------------------------------------
1039  do i=1,n(1)
1040  x = i*ds(1)
1041  a(i,:,:) = x**5
1042  end do
1043  b = xdn(a)
1044  ok = .true.
1045  do i=4,n(1)-2
1046  x = (i-0.5d0)*ds(1)
1047  d = x**5
1048  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1049  end do
1050  allok = allok .and. ok
1051  if (.not.ok) print *,'xdn:', ok
1052  b = xup(a)
1053  ok = .true.
1054  do i=3,n(1)-3
1055  x = (i+0.5d0)*ds(1)
1056  d = x**5
1057  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1058  end do
1059  allok = allok .and. ok
1060  if (.not.ok) print *,'xup:', ok
1061  !-----------------------------------------------------------------------------
1062  ! y-operators
1063  !-----------------------------------------------------------------------------
1064  do i=1,n(2)
1065  x = i*ds(2)
1066  a(:,i,:) = x**5
1067  end do
1068  b = ydn(a)
1069  ok = .true.
1070  do i=4,n(2)-2
1071  y = (i-0.5d0)*ds(2)
1072  d = y**5
1073  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1074  end do
1075  allok = allok .and. ok
1076  if (.not.ok) print *,'ydn:', ok
1077  b = yup(a)
1078  ok = .true.
1079  do i=3,n(2)-3
1080  y = (i+.5d0)*ds(2)
1081  d = y**5
1082  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1083  end do
1084  allok = allok .and. ok
1085  if (.not.ok) print *,'yup:', ok
1086  !-----------------------------------------------------------------------------
1087  ! z-operators
1088  !-----------------------------------------------------------------------------
1089  do i=1,n(3)
1090  x = i*ds(3)
1091  a(:,:,i) = x**5
1092  end do
1093  b = zdn(a)
1094  ok = .true.
1095  do i=4,n(3)-2
1096  z = (i-0.5d0)*ds(3)
1097  d = z**5
1098  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1099  end do
1100  allok = allok .and. ok
1101  if (.not.ok) print *,'zdn:', ok
1102  b = zup(a)
1103  ok = .true.
1104  do i=3,n(3)-3
1105  z = (i+0.5d0)*ds(3)
1106  d = z**5
1107  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1108  end do
1109  allok = allok .and. ok
1110  if (.not.ok) print *,'zup:', ok
1111  !-----------------------------------------------------------------------------
1112  ! x-operators
1113  !-----------------------------------------------------------------------------
1114  do i=1,n(1)
1115  x = i*ds(1)
1116  a(i,:,:) = 1d0 + x**6
1117  end do
1118  b = ddxdn(ds,a)
1119  ok = .true.
1120  do i=4,n(1)-2
1121  x = (i-0.5d0)*ds(1)
1122  d = 6d0*x**5
1123  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1124  end do
1125  allok = allok .and. ok
1126  if (.not.ok) print *,'ddxdn:', ok
1127  b = ddxup(ds,a)
1128  ok = .true.
1129  do i=3,n(1)-3
1130  x = (i+0.5d0)*ds(1)
1131  d = 6d0*x**5
1132  ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1133  end do
1134  allok = allok .and. ok
1135  if (.not.ok) print *,'ddxup:', ok
1136  !-----------------------------------------------------------------------------
1137  ! y-operators
1138  !-----------------------------------------------------------------------------
1139  do i=1,n(2)
1140  x = i*ds(2)
1141  a(:,i,:) = 1d0 + x**6
1142  end do
1143  b = ddydn(ds,a)
1144  ok = .true.
1145  do i=4,n(2)-2
1146  y = (i-0.5d0)*ds(2)
1147  d = 6d0*y**5
1148  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1149  end do
1150  allok = allok .and. ok
1151  if (.not.ok) print *,'ddydn:', ok
1152  b = ddyup(ds,a)
1153  ok = .true.
1154  do i=3,n(2)-3
1155  y = (i+.5d0)*ds(2)
1156  d = 6d0*y**5
1157  ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1158  end do
1159  allok = allok .and. ok
1160  if (.not.ok) print *,'ddyup:', ok
1161  !-----------------------------------------------------------------------------
1162  ! z-operators
1163  !-----------------------------------------------------------------------------
1164  do i=1,n(3)
1165  x = i*ds(3)
1166  a(:,:,i) = 1d0 + x**6
1167  end do
1168  b = ddzdn(ds,a)
1169  ok = .true.
1170  do i=4,n(3)-2
1171  z = (i-0.5d0)*ds(3)
1172  d = 6d0*z**5
1173  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1174  end do
1175  allok = allok .and. ok
1176  if (.not.ok) print *,'ddzdn:', ok
1177  b = ddzup(ds,a)
1178  ok = .true.
1179  do i=3,n(3)-3
1180  z = (i+0.5d0)*ds(3)
1181  d = 6d0*z**5
1182  ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1183  end do
1184  allok = allok .and. ok
1185  if (.not.ok) print *,'ddzup:', ok
1186  print *,io%hl
1187  if (allok) then
1188  print *,'6th/5th order stagger operator test passed'
1189  print *,io%hl
1190  else
1191  print *,'6th/5th order stagger operator test failed!'
1192  print *,io%hl
1193  call mpi%abort('stagger_test')
1194  end if
1195  !-----------------------------------------------------------------------------
1196  deallocate (a, b)
1197 END SUBROUTINE stagger_test
1198 
1199 END MODULE stagger_mod
Definition: io_mod.f90:4
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4