function tetrahedron_jacobian(t, x, type)
!=======================================!
! !
! Jacobian part of surface integrations !
! !
!=======================================!
use w90_constants, only: dp
implicit none
real(kind=dp) :: tetrahedron_jacobian
real(kind=dp), DIMENSION(3, 3), INTENT(IN) :: t
real(kind=dp), DIMENSION(3), INTENT(IN) :: x
real(kind=dp), DIMENSION(3, 2) :: J
real(kind=dp) :: y
real(kind=dp) :: Ans
INTEGER, INTENT(IN) :: type
INTEGER :: j_, k, a, b, c, d
IF (type == 1) THEN
J(1, 1) = -x(1); J(1, 2) = -x(1)
J(2, 1) = x(2); J(2, 2) = 0.0_dp
J(3, 1) = 0.0_dp; J(3, 2) = x(3)
ELSE IF (type == 2) THEN
J(1, 1) = x(1); J(1, 2) = 0.0_dp
J(2, 1) = 0.0_dp; J(2, 2) = x(2)
J(3, 1) = 1.0_dp - x(1) - x(3); J(3, 2) = -x(3)
ELSE IF (type == 3) THEN
y = x(1)*(x(2) - 1.0_dp)*x(3)/(-x(2) + x(1)*x(2) + x(2)*x(3) - x(1)*x(3))
J(1, 1) = x(1) - y; J(1, 2) = -y
J(2, 1) = y - 1.0_dp; J(2, 2) = y - 1.0_dp + x(2)
J(3, 1) = 1.0_dp - x(1); J(3, 2) = 0.0_dp
ELSE
J(1, 1) = x(1); J(1, 2) = 0.0_dp
J(2, 1) = 0.0_dp; J(2, 2) = x(2)
J(3, 1) = 1.0_dp - x(1) - x(3); J(3, 2) = 1.0_dp - x(2) - x(3)
END IF
Ans = 0.0_dp
DO j_ = 1, 3
DO k = 1, 3
DO a = 1, 3
DO b = 1, 3
DO c = 1, 3
DO d = 1, 3
Ans = Ans + t(j_, a)*t(j_, b)*t(k, c)*t(k, d)*J(a, 1)*J(c, 2)*(J(b, 1)*J(d, 2) - J(b, 2)*J(d, 1))
END DO
END DO
END DO
END DO
END DO
END DO
tetrahedron_jacobian = SQRT(ABS(Ans))
end function tetrahedron_jacobian
