! Finding isosurface(s).
! $Header$

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"

subroutine EHFinder_IsoSurface(CCTK_ARGUMENTS)

  use EHFinder_mod
  use IsoSurface_mod

  implicit none

  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS

  CCTK_INT :: n_coords, n_triangles
  CCTK_REAL, dimension(:), allocatable :: x_coord, y_coord, z_coord
  CCTK_INT, dimension(:), allocatable :: t_index
  CCTK_REAL, parameter :: iso_value = zero
  CCTK_INT :: l

! Find index ranges for points excluding ghost zones and symmetry points for
! the 3D grid functions.
#include "include/physical_part.h"

  call CCTK_INFO ('Entered IsoSurface')
  
  print*,nx,ny,nz
  print*,size(f)
  print*,Xf0,Xf1,Xf2,f_length
  print*,size(x)
!  allocate(x_coord(5),y_coord(5),z_coord(5),t_index(20))
  do l = 1, eh_number_level_sets
!    call triangularization ( n_coords, n_triangles, nx, ny, nz, x_coord, &
!                             y_coord, z_coord, t_index, f(:,:,:,l), x, y, z )
    call triangularization ( n_coords, n_triangles, nx, ny, nz, f(:,:,:,l), x, y, z )

  print*,'The surface consists of ',n_triangles/3,' triangles and ',n_coords,' coords'
  end do

end subroutine EHFinder_IsoSurface


!subroutine triangularization(l,contor,Nx,Ny,Nz,coord_x,coord_z,coord_y,triangle_index,grid_value,grid_x,grid_y,grid_z)
subroutine triangularization(l,contor,Nx,Ny,Nz,grid_value,grid_x,grid_y,grid_z)

  use IsoSurface_mod

  implicit none

!declaring variables that will be used in the program
  CCTK_INT,intent(in)::Nx,Ny,Nz
  CCTK_REAL:: u,ABx,ABy,ABz,BCx,BCy,BCz,dx
  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_value
  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_x
  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_y
  CCTK_REAL,dimension(Nx,Ny,Nz),intent(in)::grid_z
  CCTK_INT::i,j,k,jj,kk,index,sum,nsize,ii
  CCTK_INT,intent(out)::l,contor
!  CCTK_REAL,dimension(:),allocatable,intent(out):: coord_x,coord_y,coord_z
!  CCTK_INT,dimension (:),allocatable,intent(out)::triangle_index
!  CCTK_INT,dimension(:),allocatable::resize_triangle_index
!  CCTK_REAL,dimension(:),allocatable::resize_coord_x,resize_coord_y,resize_coord_z
!we store table with triangles in a 2 dimensional constant array using reshape function; Fortran stores column-like and we need row-like
  CCTK_INT, dimension (256,16):: a
  CCTK_INT, dimension (16,30):: b1(16,30)=reshape&
((/-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1,&
3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1,&
3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1,&
3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1,&
9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1,&
9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1,&
2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1,&
8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1,&
9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1,&
4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1,&
3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1,&
1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1/),(/16,30/))





  CCTK_INT, dimension (16,30):: b2(16,30)=reshape&
((/4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1,&
4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1,&
9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1,&
5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1,&
2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1,&
9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1,&
0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1,&
2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1,&
10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1,&
4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1,&
5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1,&
5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1,&
9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1,&
0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1,&
1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1,&
10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1,&
8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1,&
2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1,&
7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1,&
9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1,&
2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1,&
11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1/),(/16,30/))



  CCTK_INT, dimension (16,30):: b3(16,30)=reshape&
((/9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1,&
5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1,&
11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1,&
11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1,&
1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1,&
9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1,&
5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1,&
2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1,&
0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1,&
5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1,&
6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1,&
3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1,&
6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1,&
5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1,&
1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1,&
10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1,&
6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1,&
8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1,&
7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1,&
3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1,&
5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1/),(/16,30/))



  CCTK_INT, dimension (16,30):: b4(16,30)=reshape&
((/0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1,&
9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1,&
8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1,&
5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1,&
0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1,&
6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1,&
10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1,&
10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1,&
8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1,&
1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1,&
3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1,&
0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1,&
10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1,&
3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1,&
6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1,&
9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1,&
8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1,&
3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1,&
6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1,&
0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1,&
10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1,&
10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1,&
2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1,&
7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1,&
7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/),(/16,30/))




  CCTK_INT, dimension (16,30):: b5(16,30)=reshape&
((/2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1,&
2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1,&
1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1,&
11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1,&
8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1,&
0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1,&
7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1,&
10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1,&
2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1,&
6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1,&
7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1,&
2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1,&
1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1,&
10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1,&
10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1,&
0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1,&
7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1,&
6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1,&
8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1,&
9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1,&
6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1/),(/16,30/))



  CCTK_INT, dimension (16,30):: b6(16,30)=reshape&
((/4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1,&
10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1,&
8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1,&
0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1,&
1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1,&
8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1,&
10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1,&
4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1,&
10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1,&
5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1,&
11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1,&
9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1,&
6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1,&
7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1,&
3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1,&
7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1,&
9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1,&
3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1,&
6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1,&
9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1,&
1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1,&
4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1,&
7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1,&
6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1,&
3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1,&
0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1,&
6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1/),(/16,30/))


  CCTK_INT, dimension (16,30):: b7(16,30)=reshape&
((/1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1,&
0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1,&
11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1,&
6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1,&
5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1,&
9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1,&
1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1,&
1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1,&
10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1,&
0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1,&
5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1,&
10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1,&
11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1,&
9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1,&
7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1,&
2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1,&
8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1,&
9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1,&
9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1,&
1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1,&
9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1,&
9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1,&
5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1/),(/16,30/))


  CCTK_INT, dimension (16,30):: b8(16,30)=reshape&
((/0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1,&
10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1,&
2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1,&
0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1,&
0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1,&
9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1,&
5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1,&
3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1,&
5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1,&
8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1,&
0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1,&
9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1,&
0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1,&
1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1,&
3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1,&
4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1,&
9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1,&
11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1,&
11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1,&
2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1,&
9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1,&
3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1,&
1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1,&
4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1,&
4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/),(/16,30/))



  CCTK_INT, dimension (16,16):: b9(16,16)=reshape&
((/9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1,&
0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1,&
3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1,&
3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1,&
0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1,&
9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1,&
1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,&
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/),(/16,16/))

  do j=1,30
    do i=1,16
      a(j,i)=b1(i,j)
    end do 
  end do

  jj=0
  do j=31,60
    jj=jj+1
    do i=1,16
      a(j,i)=b2(i,jj)
    end do 
  end do

  jj=0
  do j=61,90
    jj=jj+1
    do i=1,16
      a(j,i)=b3(i,jj)
    end do 
  end do

  jj=0
  do j=91,120
    jj=jj+1
    do i=1,16
      a(j,i)=b4(i,jj)
    end do 
  end do

  jj=0
  do j=121,150
    jj=jj+1
    do i=1,16
      a(j,i)=b5(i,jj)
    end do 
  end do

  jj=0
  do j=151,180
    jj=jj+1
    do i=1,16
      a(j,i)=b6(i,jj)
    end do 
  end do

  jj=0
  do j=181,210
    jj=jj+1
    do i=1,16
      a(j,i)=b7(i,jj)
    end do 
  end do

  jj=0
  do j=211,240
    jj=jj+1
    do i=1,16
      a(j,i)=b8(i,jj)
    end do 
  end do

  jj=0
  do j=241,257
    jj=jj+1
    do i=1,16
      a(j,i)=b9(i,jj)
    end do 
  end do


  print*,'1'

  !allocating initial space for coordinate and triangle index arrays.These arrays will be extended if needed during the execution of the program using allocate and deallocate functions.
  nsize=5
  allocate(coord_x(nsize),coord_y(nsize),coord_z(nsize),triangle_index(nsize))

  !marching through the grid and computing "sum" for each individual cube. "l" is the contor for the vectors that keeps the coordinates of intersection points. we will also produce a vector called triangle_index that tells which points build which triangles.
  l=0
  contor=0
  do k=1,Nz-1
    do j=1,Ny-1
      do i=1,Nx-1                                 
        sum=1
          if (grid_value(i,j,k)<0) then 
            sum=sum+2**0
          end if
          if (grid_value(i+1,j,k)<0) then 
            sum=sum+2**3
          end if
          if (grid_value(i+1,j+1,k)<0) then 
            sum=sum+2**2
          end if
          if (grid_value(i,j+1,k)<0) then 
            sum=sum+2**1
          end if
          if (grid_value(i+1,j+1,k+1)<0) then 
            sum=sum+2**6
          end if
          if (grid_value(i+1,j,k+1)<0) then 
            sum=sum+2**7
          end if
          if (grid_value(i,j,k+1)<0) then 
            sum=sum+2**4
          end if
          if (grid_value(i,j+1,k+1)<0) then 
            sum=sum+2**5
          end if

          if (a(sum,1)>-1) then
            do jj=1,15,3
              if (a(sum,jj)>-1) then
                do kk=jj,jj+2

                  if (a(sum,kk)==0) then 
                    u=abs(grid_value(i,j,k))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i,j,k))+abs(grid_value(i,j+1,k)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)
                    coord_y(l)=grid_y(i,j,k)+u
                    coord_z(l)=grid_z(i,j,k)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==1) then 
                    u=abs(grid_value(i,j+1,k))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i,j+1,k))+abs(grid_value(i+1,j+1,k)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)+u
                    coord_y(l)=grid_y(i,j+1,k)
                    coord_z(l)=grid_z(i,j,k)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==2) then 
                    u=abs(grid_value(i+1,j,k))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i+1,j+1,k))+abs(grid_value(i+1,j,k)))
                    l=l+1
                    coord_x(l)=grid_x(i+1,j,k)
                    coord_y(l)=grid_y(i,j,k)+u
                    coord_z(l)=grid_z(i,j,k)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==3) then
                    u=abs(grid_value(i,j,k))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i+1,j,k))+abs(grid_value(i,j,k)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)+u
                    coord_y(l)=grid_y(i,j,k)
                    coord_z(l)=grid_z(i,j,k)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==4) then 
                    u=abs(grid_value(i,j,k+1))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i,j,k+1))+abs(grid_value(i,j+1,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)
                    coord_y(l)=grid_y(i,j,k)+u
                    coord_z(l)=grid_z(i,j,k+1)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==5) then 
                    u=abs(grid_value(i,j+1,k+1))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i,j+1,k+1))+abs(grid_value(i+1,j+1,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)+u
                    coord_y(l)=grid_y(i,j+1,k)
                    coord_z(l)=grid_z(i,j,k+1)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==6) then 
                    u=abs(grid_value(i+1,j,k+1))*(grid_y(i,j+1,k)-grid_y(i,j,k))/(abs(grid_value(i+1,j+1,k+1))+abs(grid_value(i+1,j,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i+1,j,k)
                    coord_y(l)=grid_y(i,j,k)+u
                    coord_z(l)=grid_z(i,j,k+1)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if
 
                  if (a(sum,kk)==7) then
                    u=abs(grid_value(i,j,k+1))*(grid_x(i+1,j,k)-grid_x(i,j,k))/(abs(grid_value(i,j,k+1))+abs(grid_value(i+1,j,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)+u
                    coord_y(l)=grid_y(i,j,k)
                    coord_z(l)=grid_z(i,j,k+1)
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==8) then 
                    u=abs(grid_value(i,j,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i,j,k))+abs(grid_value(i,j,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)
                    coord_y(l)=grid_y(i,j,k)
                    coord_z(l)=grid_z(i,j,k)+u
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==9) then 
                    u=abs(grid_value(i,j+1,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i,j+1,k))+abs(grid_value(i,j+1,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i,j,k)
                    coord_y(l)=grid_y(i,j+1,k)
                    coord_z(l)=grid_z(i,j,k)+u
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==10) then 
                    u=abs(grid_value(i+1,j+1,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i+1,j+1,k))+abs(grid_value(i+1,j+1,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i+1,j,k)
                    coord_y(l)=grid_y(i,j+1,k)
                    coord_z(l)=grid_z(i,j,k)+u
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if

                  if (a(sum,kk)==11) then 
                    u=abs(grid_value(i+1,j,k))*(grid_z(i,j,k+1)-grid_z(i,j,k))/(abs(grid_value(i+1,j,k))+abs(grid_value(i+1,j,k+1)))
                    l=l+1
                    coord_x(l)=grid_x(i+1,j,k)
                    coord_y(l)=grid_y(i,j,k)
                    coord_z(l)=grid_z(i,j,k)+u
                    contor=contor+1
                    triangle_index(contor)=l-1
                  end if
                end do

                  if (l>3) then
                    ABx=coord_x(l)-coord_x(l-1)
                    BCx=coord_x(l-1)-coord_x(l-2)
                    ABy=coord_y(l)-coord_y(l-1)
                    BCy=coord_y(l-1)-coord_y(l-2)
                    ABz=coord_z(l)-coord_z(l-1)
                    BCz=coord_z(l-1)-coord_z(l-2)
                  if (((ABx*BCy)==(ABy*BCx)).and.((ABx*BCz)==(ABz*BCx)).and.((ABy*BCz)==(ABz*BCy))) then
                    l=l-3 
                    contor=contor-3
                  end if
                end if
              end if

              if ((nsize-l)<10) then
                allocate(resize_coord_x(nsize),resize_coord_y(nsize),resize_coord_z(nsize),resize_triangle_index(nsize))
                resize_coord_x(1:l)=coord_x(1:l)
                resize_coord_y(1:l)=coord_y(1:l)
                resize_coord_z(1:l)=coord_z(1:l)
                resize_triangle_index(1:l)=triangle_index(1:l)
                deallocate(coord_x,coord_y,coord_z,triangle_index)
                nsize=nsize*2
                allocate(coord_x(nsize),coord_y(nsize),coord_z(nsize),triangle_index(nsize))
                coord_x(1:l)=resize_coord_x(1:l)
                coord_y(1:l)=resize_coord_y(1:l)
                coord_z(1:l)=resize_coord_z(1:l)
                triangle_index(1:l)=resize_triangle_index(1:l)
                deallocate(resize_coord_x,resize_coord_y,resize_coord_z,resize_triangle_index)
            end if
          end do
        end if
      end do
    end do
  end do 
end subroutine triangularization