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
383 FUNCTION ddxdn (ds, a)
RESULT (b)
386 real :: b(size(a,1),size(a,2),size(a,3))
388 integer :: i, j, k, n(3)
389 real :: c1, c2, c3, d
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
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))
413 stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
414 stagger%count = stagger%count + 1
418 FUNCTION ddydn (ds, a)
RESULT (b)
421 real :: b(size(a,1),size(a,2),size(a,3))
423 integer :: i, j, k, n(3)
424 real :: c1, c2, c3, d
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
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))
448 stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
449 stagger%count = stagger%count + 1
453 FUNCTION ddzdn (ds, a)
RESULT (b)
456 real :: b(size(a,1),size(a,2),size(a,3))
458 integer :: i, j, k, n(3)
459 real :: c1, c2, c3, d
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
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))
483 stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
484 stagger%count = stagger%count + 1
488 FUNCTION ddxup (ds, a)
RESULT (b)
491 real :: b(size(a,1),size(a,2),size(a,3))
493 integer :: i, j, k, n(3)
494 real :: c1, c2, c3, d
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
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))
518 stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
519 stagger%count = stagger%count + 1
523 FUNCTION ddyup (ds, a)
RESULT (b)
526 real :: b(size(a,1),size(a,2),size(a,3))
528 integer :: i, j, k, n(3)
529 real :: c1, c2, c3, d
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
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))
553 stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
554 stagger%count = stagger%count + 1
558 FUNCTION ddzup (ds, a)
RESULT (b)
561 real :: b(size(a,1),size(a,2),size(a,3))
563 integer :: i, j, k, n(3)
564 real :: c1, c2, c3, d
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
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 ))
588 stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
589 stagger%count = stagger%count + 1
593 FUNCTION xdn (a)
RESULT (b)
595 real :: b(size(a,1),size(a,2),size(a,3))
597 integer :: i, j, k, n(3)
598 real :: c1, c2, c3, d
605 b( 1:3,:,:) = a( 1:3,:,:)
606 b(n(1)-1:,:,:) = a(n(1)-1:,:,:)
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))
619 stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
620 stagger%count = stagger%count + 1
624 FUNCTION ydn (a)
RESULT (b)
626 real :: b(size(a,1),size(a,2),size(a,3))
628 integer :: i, j, k, n(3)
629 real :: c1, c2, c3, d
636 b(:, 1:3,:) = a(:, 1:3,:)
637 b(:,n(2)-1:,:) = a(:,n(2)-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))
650 stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
651 stagger%count = stagger%count + 1
655 FUNCTION zdn (a)
RESULT (b)
657 real :: b(size(a,1),size(a,2),size(a,3))
659 integer :: i, j, k, n(3)
660 real :: c1, c2, c3, d
667 b(:,:, 1:3) = a(:,:, 1:3)
668 b(:,:,n(1)-1:) = a(:,:,n(1)-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))
681 stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
682 stagger%count = stagger%count + 1
686 FUNCTION xup (a)
RESULT (b)
688 real :: b(size(a,1),size(a,2),size(a,3))
690 integer :: i, j, k, n(3)
691 real :: c1, c2, c3, d
698 b( 1:2,:,:) = a( 1:2,:,:)
699 b(n(1)-2:,:,:) = a(n(1)-2:,:,:)
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))
712 stagger%flops = stagger%flops + 8*n(2)*n(3)*(n(1)-5)
713 stagger%count = stagger%count + 1
717 FUNCTION yup (a)
RESULT (b)
719 real :: b(size(a,1),size(a,2),size(a,3))
721 integer :: i, j, k, n(3)
722 real :: c1, c2, c3, d
729 b(:, 1:2,:) = a(:, 1:2,:)
730 b(:,n(2)-2:,:) = a(:,n(2)-2:,:)
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))
743 stagger%flops = stagger%flops + 8*n(1)*n(3)*(n(2)-5)
744 stagger%count = stagger%count + 1
748 FUNCTION zup (a)
RESULT (b)
750 real :: b(size(a,1),size(a,2),size(a,3))
752 integer :: i, j, k, n(3)
753 real :: c1, c2, c3, d
760 b(:,:, 1:2) = a(:,:, 1:2)
761 b(:,:,n(1)-2:) = a(:,:,n(1)-2:)
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 ))
774 stagger%flops = stagger%flops + 8*n(1)*n(2)*(n(3)-5)
775 stagger%count = stagger%count + 1
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
784 integer,
save:: itimer=0
791 b(i,:,:) = (a(i-1,:,:) + a(i+1,:,:))*c + a(i,:,:)*(1.-2.*c)
800 stagger%flops = stagger%flops + 4*n*n*(n-2)
801 stagger%count = stagger%count + 1
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
810 integer,
save:: itimer=0
817 b(:,i,:) = (a(:,i-1,:) + a(:,i+1,:))*c + a(:,i,:)*(1.-2.*c)
826 stagger%flops = stagger%flops + 4*n*n*(n-2)
827 stagger%count = stagger%count + 1
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
836 integer,
save:: itimer=0
843 b(:,:,i) = (a(:,:,i-1) + a(:,:,i+1))*c + a(:,:,i)*(1.-2.*c)
852 stagger%flops = stagger%flops + 4*n*n*(n-2)
853 stagger%count = stagger%count + 1
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)
863 call trace%end (itimer)
869 SUBROUTINE stagger_test
873 real,
pointer :: a(:,:,:)
874 real,
pointer :: b(:,:,:)
875 real,
pointer :: c(:,:,:)
876 real(8) :: ds(3), x, y, z, d
879 real,
parameter :: eps = 2e-6
883 allocate (a(n(1),n(2),n(3)))
884 allocate (b(n(1),n(2),n(3)))
898 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
900 allok = allok .and. ok
901 if (.not.ok) print *,
'ddxdn1:', ok
907 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
909 allok = allok .and. ok
910 if (.not.ok) print *,
'ddxup1:', ok
923 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
925 allok = allok .and. ok
926 if (.not.ok) print *,
'ddydn1:', ok
932 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
934 allok = allok .and. ok
935 if (.not.ok) print *,
'ddyup1:', ok
948 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
950 allok = allok .and. ok
951 if (.not.ok) print *,
'ddzdn1:', ok
957 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
959 allok = allok .and. ok
960 if (.not.ok) print *,
'ddzup1:', ok
973 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
975 allok = allok .and. ok
976 if (.not.ok) print *,
'xdn1:', ok
982 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
984 allok = allok .and. ok
985 if (.not.ok) print *,
'xup1:', ok
998 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1000 allok = allok .and. ok
1001 if (.not.ok) print *,
'ydn1:', ok
1007 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1009 allok = allok .and. ok
1010 if (.not.ok) print *,
'yup1:', ok
1023 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1025 allok = allok .and. ok
1026 if (.not.ok) print *,
'zdn1:', ok
1032 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1034 allok = allok .and. ok
1035 if (.not.ok) print *,
'zup1:', ok
1048 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1050 allok = allok .and. ok
1051 if (.not.ok) print *,
'xdn:', ok
1057 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1059 allok = allok .and. ok
1060 if (.not.ok) print *,
'xup:', ok
1073 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1075 allok = allok .and. ok
1076 if (.not.ok) print *,
'ydn:', ok
1082 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1084 allok = allok .and. ok
1085 if (.not.ok) print *,
'yup:', ok
1098 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1100 allok = allok .and. ok
1101 if (.not.ok) print *,
'zdn:', ok
1107 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1109 allok = allok .and. ok
1110 if (.not.ok) print *,
'zup:', ok
1116 a(i,:,:) = 1d0 + x**6
1123 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1125 allok = allok .and. ok
1126 if (.not.ok) print *,
'ddxdn:', ok
1132 ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
1134 allok = allok .and. ok
1135 if (.not.ok) print *,
'ddxup:', ok
1141 a(:,i,:) = 1d0 + x**6
1148 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1150 allok = allok .and. ok
1151 if (.not.ok) print *,
'ddydn:', ok
1157 ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
1159 allok = allok .and. ok
1160 if (.not.ok) print *,
'ddyup:', ok
1166 a(:,:,i) = 1d0 + x**6
1173 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1175 allok = allok .and. ok
1176 if (.not.ok) print *,
'ddzdn:', ok
1182 ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
1184 allok = allok .and. ok
1185 if (.not.ok) print *,
'ddzup:', ok
1188 print *,
'6th/5th order stagger operator test passed' 1191 print *,
'6th/5th order stagger operator test failed!' 1193 call mpi%abort(
'stagger_test')
1197 END SUBROUTINE stagger_test
6th order stagger operators, with self-test procedure