14 integer,
private:: verbose=0
15 logical,
save:: hardwire=.false.
18 integer(kind=8):: flops=0
19 integer(kind=8):: count=0
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
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
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
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
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
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
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
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
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
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
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
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
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
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
152 integer :: ix, iy, iz, n(3)
153 integer,
save:: itimer=0
178 b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix-1,iy,iz))
186 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
187 stagger%count = stagger%count + 1
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
195 integer :: ix, iy, iz, n(3)
196 integer,
save:: itimer=0
221 b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy-1,iz))
229 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
230 stagger%count = stagger%count + 1
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
238 integer :: ix, iy, iz, n(3)
239 integer,
save:: itimer=0
264 b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz-1))
272 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
273 stagger%count = stagger%count + 1
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
281 integer :: ix, iy, iz, n(3)
282 integer,
save:: itimer=0
306 b(ix,iy,iz) = c*(a(ix ,iy,iz) - a(ix-1,iy,iz))
314 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
315 stagger%count = stagger%count + 1
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)
324 integer,
save:: itimer=0
348 b(ix,iy,iz) = c*(a(ix,iy ,iz) - a(ix,iy-1,iz))
356 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
357 stagger%count = stagger%count + 1
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
365 integer :: ix, iy, iz, n(3)
366 integer,
save:: itimer=0
390 b(ix,iy,iz) = c*(a(ix,iy,iz ) - a(ix,iy,iz-1))
398 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
399 stagger%count = stagger%count + 1
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)
408 integer,
save:: itimer=0
433 b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix ,iy,iz))
441 stagger%flops = stagger%flops + 2*n(3)*n(2)*(n(1)-1)
442 stagger%count = stagger%count + 1
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)
451 integer,
save:: itimer=0
475 b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy ,iz))
483 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
484 stagger%count = stagger%count + 1
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)
493 integer,
save:: itimer=0
517 b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz ))
525 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
526 stagger%count = stagger%count + 1
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
534 integer :: ix, iy, iz, n(3)
535 integer,
save:: itimer=0
540 if (all(n==[16,16,16]))
then 542 else if (all(n==[20,20,20]))
then 544 else if (all(n==[24,24,24]))
then 546 else if (all(n==[32,32,32]))
then 548 else if (all(n==[36,36,36]))
then 551 else if (n(1) > 1)
then 558 b(ix,iy,iz) = c*(a(ix ,iy,iz) + a(ix-1,iy,iz))
566 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
567 stagger%count = stagger%count + 1
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
575 integer :: ix, iy, iz, n(3)
576 integer,
save:: itimer=0
581 if (all(n==[16,16,16]))
then 583 else if (all(n==[20,20,20]))
then 585 else if (all(n==[24,24,24]))
then 587 else if (all(n==[32,32,32]))
then 589 else if (all(n==[36,36,36]))
then 592 else if (n(2) > 1)
then 599 b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy-1,iz))
607 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
608 stagger%count = stagger%count + 1
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
616 integer :: ix, iy, iz, n(3)
617 integer,
save:: itimer=0
622 if (all(n==[16,16,16]))
then 624 else if (all(n==[20,20,20]))
then 626 else if (all(n==[24,24,24]))
then 628 else if (all(n==[32,32,32]))
then 630 else if (all(n==[36,36,36]))
then 633 else if (n(3) > 1)
then 640 b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy,iz-1))
648 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
649 stagger%count = stagger%count + 1
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
657 integer :: ix, iy, iz, n(3)
658 integer,
save:: itimer=0
663 if (all(n==[16,16,16]))
then 665 else if (all(n==[20,20,20]))
then 667 else if (all(n==[24,24,24]))
then 669 else if (all(n==[32,32,32]))
then 671 else if (all(n==[36,36,36]))
then 674 else if (n(1) > 1)
then 676 b(n(1),:,:) = a(n(1),:,:)
681 b(ix,iy,iz) = c*(a(ix+1,iy,iz) + a(ix ,iy,iz))
689 stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
690 stagger%count = stagger%count + 1
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
698 integer :: ix, iy, iz, n(3)
699 integer,
save:: itimer=0
704 if (all(n==[16,16,16]))
then 706 else if (all(n==[20,20,20]))
then 708 else if (all(n==[24,24,24]))
then 710 else if (all(n==[32,32,32]))
then 712 else if (all(n==[36,36,36]))
then 715 else if (n(2) > 1)
then 717 b(:,n(2),:) = a(:,n(2),:)
722 b(ix,iy,iz) = c*(a(ix,iy+1,iz) + a(ix,iy ,iz))
730 stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
731 stagger%count = stagger%count + 1
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
739 integer :: ix, iy, iz, n(3)
740 integer,
save:: itimer=0
745 if (all(n==[16,16,16]))
then 747 else if (all(n==[20,20,20]))
then 749 else if (all(n==[24,24,24]))
then 751 else if (all(n==[32,32,32]))
then 753 else if (all(n==[36,36,36]))
then 756 else if (n(3) > 1)
then 758 b(:,:,n(3)) = a(:,:,n(3))
763 b(ix,iy,iz) = c*(a(ix,iy,iz+1) + a(ix,iy,iz ))
771 stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
772 stagger%count = stagger%count + 1
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
781 integer,
save:: itimer=0
788 b(i,:,:) = (a(i-1,:,:) + a(i+1,:,:))*c + a(i,:,:)*(1.-2.*c)
797 stagger%flops = stagger%flops + 4*n*n*(n-2)
798 stagger%count = stagger%count + 1
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
807 integer,
save:: itimer=0
814 b(:,i,:) = (a(:,i-1,:) + a(:,i+1,:))*c + a(:,i,:)*(1.-2.*c)
823 stagger%flops = stagger%flops + 4*n*n*(n-2)
824 stagger%count = stagger%count + 1
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
833 integer,
save:: itimer=0
840 b(:,:,i) = (a(:,:,i-1) + a(:,:,i+1))*c + a(:,:,i)*(1.-2.*c)
849 stagger%flops = stagger%flops + 4*n*n*(n-2)
850 stagger%count = stagger%count + 1
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)
860 call trace%end (itimer)
866 SUBROUTINE stagger_test
868 real,
pointer :: a(:,:,:)
869 real,
pointer :: b(:,:,:)
870 real,
pointer :: c(:,:,:)
871 real(8) :: ds(3), x, y, z, d
874 real,
parameter :: eps = 2e-6
878 allocate (a(m(1),m(2),m(3)))
879 allocate (b(m(1),m(2),m(3)))
893 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
895 allok = allok .and. ok
896 if (.not.ok) print *,
'xdn:', ok
902 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
904 allok = allok .and. ok
905 if (.not.ok) print *,
'xup:', ok
918 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
920 allok = allok .and. ok
921 if (.not.ok) print *,
'ydn:', ok
927 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
929 allok = allok .and. ok
930 if (.not.ok) print *,
'yup:', ok
943 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
945 allok = allok .and. ok
946 if (.not.ok) print *,
'zdn:', ok
952 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
954 allok = allok .and. ok
955 if (.not.ok) print *,
'zup:', ok
961 a(i,:,:) = 1d0 + x**2
968 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
970 allok = allok .and. ok
971 if (.not.ok) print *,
'ddxdn:', ok
977 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
979 allok = allok .and. ok
980 if (.not.ok) print *,
'ddxup:', ok
986 a(:,i,:) = 1d0 + x**2
993 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
995 allok = allok .and. ok
996 if (.not.ok) print *,
'ddydn:', ok
1002 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1004 allok = allok .and. ok
1005 if (.not.ok) print *,
'ddyup:', ok
1011 a(:,:,i) = 1d0 + x**2
1018 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1020 allok = allok .and. ok
1021 if (.not.ok) print *,
'ddzdn:', ok
1027 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1029 allok = allok .and. ok
1030 if (.not.ok) print *,
'ddzup:', ok
1033 print *,
'2nd/1st order stagger operator test passed' 1036 print *,
'2nd/1st order stagger operator test failed!' 1038 error stop
'stagger_test' 1042 END SUBROUTINE stagger_test
6th order stagger operators, with self-test procedure