subroutine tetrahedron_array_init(tet_array)
implicit none
integer, dimension(6, 20), intent(inout) :: tet_array
! with correction, Kawamura PRB 89 094515
! z !
! ! 61--62--63-64 ! tetrahedron #1: 22-38-39-43
tet_array(1, 1:4) = (/22, 38, 39, 43/) ! /| /| layer 4! tetrahedron #2: 22-42-38-43
tet_array(2, 1:4) = (/22, 42, 38, 43/) ! 49--50--51-52| ! tetrahedron #3: 22-26-42-43
tet_array(3, 1:4) = (/22, 26, 42, 43/) ! | | y | | ! tetrahedron #4: 22-27-26-43
tet_array(4, 1:4) = (/22, 27, 26, 43/) ! | 13--14--1|-16 ! tetrahedron #5: 22-23-27-43
tet_array(5, 1:4) = (/22, 23, 27, 43/) ! |/ |/ layer 1! tetrahedron #6: 22-39-23-43
tet_array(6, 1:4) = (/22, 39, 23, 43/) ! 1---2---3--4 x !
! 22-43 : main diagonal !
!=============================!
tet_array(1, 5:20) = (/6, 37, 35, 64, 5, 33, 56, 48, 1, 54, 40, 47, 59, 23, 42, 18/)
tet_array(2, 5:20) = (/2, 46, 33, 64, 6, 41, 54, 44, 1, 62, 34, 48, 63, 18, 47, 17/)
tet_array(3, 5:20) = (/18, 10, 41, 64, 2, 9, 62, 60, 1, 30, 58, 44, 47, 38, 27, 21/)
tet_array(4, 5:20) = (/17, 28, 9, 64, 18, 11, 30, 59, 1, 32, 25, 60, 48, 21, 44, 5/)
tet_array(5, 5:20) = (/21, 19, 11, 64, 17, 3, 32, 63, 1, 24, 31, 59, 44, 26, 39, 6/)
tet_array(6, 5:20) = (/5, 55, 3, 64, 21, 35, 24, 47, 1, 56, 7, 63, 60, 6, 59, 2/)
!note: 18 points(4, 8, 12, 13, 14, 15, 16, 20, 29, 36, 45, 49, 50, 51, 52, 53, 57, 61) are not used
end subroutine tetrahedron_array_init
