9 integer,
private:: verbose=0
10 logical,
save:: hardwire=.false.
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
46 FUNCTION ddxdn1 (ds, a)
RESULT (b)
47 real,
dimension(:,:,:):: a
48 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
50 integer :: ix, iy, iz, n(3)
51 integer,
save:: itimer=0
61 b(ix,iy,iz) = c*(a(ix ,iy,iz) - a(ix-1,iy,iz))
69 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
70 stagger%count = stagger%count + 1
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)
79 integer,
save:: itimer=0
89 b(ix,iy,iz) = c*(a(ix,iy ,iz) - a(ix,iy-1,iz))
97 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
98 stagger%count = stagger%count + 1
102 FUNCTION ddzdn1 (ds, a)
RESULT (b)
103 real,
dimension(:,:,:):: a
104 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
106 integer :: ix, iy, iz, n(3)
107 integer,
save:: itimer=0
117 b(ix,iy,iz) = c*(a(ix,iy,iz ) - a(ix,iy,iz-1))
125 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
126 stagger%count = stagger%count + 1
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)
135 integer,
save:: itimer=0
146 b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix ,iy,iz))
154 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
155 stagger%count = stagger%count + 1
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)
164 integer,
save:: itimer=0
174 b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy ,iz))
182 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
183 stagger%count = stagger%count + 1
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)
192 integer,
save:: itimer=0
202 b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz ))
210 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
211 stagger%count = stagger%count + 1
215 FUNCTION xdn1 (a)
RESULT (b)
216 real,
dimension(:,:,:):: a
217 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
219 integer :: ix, iy, iz, n(3)
220 integer,
save:: itimer=0
230 b(ix,iy,iz) = c*(a(ix ,iy,iz) + a(ix-1,iy,iz))
238 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
239 stagger%count = stagger%count + 1
243 FUNCTION ydn1 (a)
RESULT (b)
244 real,
dimension(:,:,:):: a
245 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
247 integer :: ix, iy, iz, n(3)
248 integer,
save:: itimer=0
258 b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy-1,iz))
266 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
267 stagger%count = stagger%count + 1
271 FUNCTION zdn1 (a)
RESULT (b)
272 real,
dimension(:,:,:):: a
273 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
275 integer :: ix, iy, iz, n(3)
276 integer,
save:: itimer=0
286 b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy,iz-1))
294 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
295 stagger%count = stagger%count + 1
299 FUNCTION xup1 (a)
RESULT (b)
300 real,
dimension(:,:,:):: a
301 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
303 integer :: ix, iy, iz, n(3)
304 integer,
save:: itimer=0
310 b(n(1),:,:) = a(n(1),:,:)
314 b(ix,iy,iz) = c*(a(ix+1,iy,iz) + a(ix ,iy,iz))
322 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
323 stagger%count = stagger%count + 1
327 FUNCTION yup1 (a)
RESULT (b)
328 real,
dimension(:,:,:):: a
329 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
331 integer :: ix, iy, iz, n(3)
332 integer,
save:: itimer=0
338 b(:,n(2),:) = a(:,n(2),:)
342 b(ix,iy,iz) = c*(a(ix,iy+1,iz) + a(ix,iy ,iz))
350 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
351 stagger%count = stagger%count + 1
355 FUNCTION zup1 (a)
RESULT (b)
356 real,
dimension(:,:,:):: a
357 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
359 integer :: ix, iy, iz, n(3)
360 integer,
save:: itimer=0
366 b(:,:,n(3)) = a(:,:,n(3))
370 b(ix,iy,iz) = c*(a(ix,iy,iz+1) + a(ix,iy,iz ))
378 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
379 stagger%count = stagger%count + 1
384 FUNCTION ddxdn (ds, a)
RESULT (b)
385 real,
dimension(:,:,:):: a
386 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
389 integer :: ix, iy, iz, n(3)
392 if (
size(a,1) > 1)
then 394 b(
size(a,1),:,:) = 0.0
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))
408 stagger%flops = stagger%flops + 5*n(2)*n(3)*(n(1)-1)
409 stagger%count = stagger%count + 1
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)
421 if (
size(a,2) > 1)
then 423 b(:,
size(a,2),:) = 0.0
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))
437 stagger%flops = stagger%flops + 5*n(1)*n(3)*(n(2)-3)
438 stagger%count = stagger%count + 1
442 FUNCTION ddzdn (ds, a)
RESULT (b)
443 real,
dimension(:,:,:):: a
444 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
447 integer :: ix, iy, iz, n(3)
450 if (
size(a,3) > 1)
then 452 b(:,:,
size(a,3)) = 0.0
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))
466 stagger%flops = stagger%flops + 5*n(1)*n(2)*(n(3)-3)
467 stagger%count = stagger%count + 1
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)
478 if (
size(a,1) > 1)
then 480 b(
size(a,1)-1:,:,:) = 0.0
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))
494 stagger%flops = stagger%flops + 5*n(2)*n(3)*(n(1)-1)
495 stagger%count = stagger%count + 1
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)
507 if (
size(a,2) > 1)
then 509 b(:,
size(a,2)-1:,:) = 0.0
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))
523 stagger%flops = stagger%flops + 5*n(1)*n(3)*(n(2)-3)
524 stagger%count = stagger%count + 1
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)
536 if (
size(a,3) > 1)
then 538 b(:,:,
size(a,3)-1:) = 0.0
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))
552 stagger%flops = stagger%flops + 5*n(1)*n(2)*(n(3)-3)
553 stagger%count = stagger%count + 1
557 FUNCTION xdn (a)
RESULT (b)
558 real,
dimension(:,:,:):: a
559 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
561 integer :: ix, iy, iz, n(3)
564 if (
size(a,1) > 1)
then 567 b(1:2,:,:) = a(1:2,:,:)
568 b(
size(a,1),:,:) = a(
size(a,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))
580 stagger%flops = stagger%flops + 5*n(2)*n(3)*(n(1)-1)
581 stagger%count = stagger%count + 1
585 FUNCTION ydn (a)
RESULT (b)
586 real,
dimension(:,:,:):: a
587 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
589 integer :: ix, iy, iz, n(3)
592 if (
size(a,2) > 1)
then 595 b(:,1:2,:) = a(:,1:2,:)
596 b(:,
size(a,2),:) = a(:,
size(a,2),:)
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))
608 stagger%flops = stagger%flops + 5*n(1)*n(3)*(n(2)-3)
609 stagger%count = stagger%count + 1
613 FUNCTION zdn (a)
RESULT (b)
614 real,
dimension(:,:,:):: a
615 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
617 integer :: ix, iy, iz, n(3)
620 if (
size(a,3) > 1)
then 623 b(:,:,1:2) = a(:,:,1:2)
624 b(:,:,
size(a,3)) = a(:,:,
size(a,3))
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))
636 stagger%flops = stagger%flops + 5*n(1)*n(2)*(n(3)-3)
637 stagger%count = stagger%count + 1
641 FUNCTION xup (a)
RESULT (b)
642 real,
dimension(:,:,:):: a
643 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
645 integer :: ix, iy, iz, n(3)
648 if (
size(a,1) > 1)
then 652 b(
size(a,1)-1:,:,:) = a(
size(a,1)-1:,:,:)
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))
664 stagger%flops = stagger%flops + 5*n(2)*n(3)*(n(1)-1)
665 stagger%count = stagger%count + 1
669 FUNCTION yup (a)
RESULT (b)
670 real,
dimension(:,:,:):: a
671 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
673 integer :: ix, iy, iz, n(3)
676 if (
size(a,2) > 1)
then 680 b(:,
size(a,2)-1:,:) = a(:,
size(a,2)-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))
692 stagger%flops = stagger%flops + 5*n(1)*n(3)*(n(2)-3)
693 stagger%count = stagger%count + 1
697 FUNCTION zup (a)
RESULT (b)
698 real,
dimension(:,:,:):: a
699 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
701 integer :: ix, iy, iz, n(3)
704 if (
size(a,3) > 1)
then 708 b(:,:,
size(a,3)-1:) = a(:,:,
size(a,3)-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))
720 stagger%flops = stagger%flops + 5*n(1)*n(2)*(n(3)-3)
721 stagger%count = stagger%count + 1
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
730 integer,
save:: itimer=0
737 b(i,:,:) = (a(i-1,:,:) + a(i+1,:,:))*c + a(i,:,:)*(1.-2.*c)
746 stagger%flops = stagger%flops + 4*n*n*(n-2)
747 stagger%count = stagger%count + 1
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
756 integer,
save:: itimer=0
763 b(:,i,:) = (a(:,i-1,:) + a(:,i+1,:))*c + a(:,i,:)*(1.-2.*c)
772 stagger%flops = stagger%flops + 4*n*n*(n-2)
773 stagger%count = stagger%count + 1
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
782 integer,
save:: itimer=0
789 b(:,:,i) = (a(:,:,i-1) + a(:,:,i+1))*c + a(:,:,i)*(1.-2.*c)
798 stagger%flops = stagger%flops + 4*n*n*(n-2)
799 stagger%count = stagger%count + 1
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)
809 call trace%end (itimer)
815 SUBROUTINE stagger_test
817 real,
pointer :: a(:,:,:)
818 real,
pointer :: b(:,:,:)
819 real,
pointer :: c(:,:,:)
820 real(8) :: ds(3), x, y, z, d
823 real,
parameter :: eps = 2e-6
827 allocate (a(m(1),m(2),m(3)))
828 allocate (b(m(1),m(2),m(3)))
842 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
844 if (.not.ok) print *,
'xdn:', ok
851 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
853 if (.not.ok) print *,
'xup:', ok
867 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
869 if (.not.ok) print *,
'ydn:', ok
876 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
878 if (.not.ok) print *,
'yup:', ok
892 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
894 if (.not.ok) print *,
'zdn:', ok
901 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
903 if (.not.ok) print *,
'zup:', ok
910 a(i,:,:) = 1d0 + x**4
917 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
919 if (.not.ok) print *,
'ddxdn:', ok
926 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
928 if (.not.ok) print *,
'ddxup:', ok
935 a(:,i,:) = 1d0 + x**4
942 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
944 if (.not.ok) print *,
'ddydn:', ok
951 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
953 if (.not.ok) print *,
'ddyup:', ok
960 a(:,:,i) = 1d0 + x**4
967 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
969 if (.not.ok) print *,
'ddzdn:', ok
976 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
978 if (.not.ok) print *,
'ddzup:', ok
982 print *,
'4th/3rd order stagger operator test passed' 985 print *,
'4th/3rd order stagger operator test failed!' 987 error stop
'stagger_test' 991 END SUBROUTINE stagger_test
6th order stagger operators, with self-test procedure