pure function deye_indices(i,j) result(a) integer, intent(in) :: i,j real(8) :: a a = 0d0 if(i==j)a=1d0 end function deye_indices