!
! utility/help subroutines
!

	Module Util_Mod
	Use GlobalVars_Mod
        Use Mesh_Mod
	Use MemAlloc_Mod
	Use MPI_Mod
	implicit none
        save

	Contains

	! Proc 0 reads in a global mesh from a file and distributes
        ! mesh elements to all procs
	!
	Subroutine Read_dist_mesh(filename, myPid, totalProcs, &
	                loc_elmts, elmt_node_list, elmt_edge_list, &
		        node_coord_list, elmt_ngb_list)
	character(len=*) :: filename
	integer, intent(in) :: myPid, totalProcs
	integer :: loc_elmts
        integer, dimension(:), pointer :: elmt_node_list
        integer, dimension(:), pointer :: elmt_edge_list
        integer, dimension(:), pointer :: elmt_ngb_list
	real, dimension(:), pointer :: node_coord_list

	character(len=100) :: ptitle
	integer, dimension(:), pointer :: edge_node_list
	integer :: tnodes, telmts, tedges, iholder, idx, jdx, kdx, istart
        integer :: nelmts_send, nelmts_recvd, nelmts_local, hdl(5), msg_tag
	integer, dimension(:), pointer :: tmp_array
	integer, dimension(:), pointer :: iptr
	integer, target :: ibuf(2)
	integer :: kk

	if(myPid == 0) then
	  open(unit=6, file=filename, action="READ")

	  ! Skip first string line
   	  read(unit=6,FMT=100) ptitle

	  ! Read in total nodes, elemts, and edges
	  read(unit=6,fmt=*) tnodes, telmts, tedges, iholder
	  
	  if(debug_) then
            print*,'tnodes = ', tnodes
            print*,'telmts = ', telmts
	  endif

	  ! Skip another string line
	  read(unit=6,FMT=100) ptitle

	  ! Temp node_coord, elmt_node and elmt_edge arrays 
	  allocate(node_coord_list(2*tnodes), stat=err_code)
	  allocate(elmt_node_list(3*telmts), stat=err_code)
	  allocate(elmt_edge_list(3*telmts), stat=err_code)
	  allocate(edge_node_list(3*tedges), stat=err_code)

	  ! Read in all node coords
	  Do idx=1, tnodes
	    read(unit=6,fmt=*) iholder ! skip another line
	    read(unit=6,fmt=*) node_coord_list(2*idx-1), node_coord_list(2*idx)
	  EndDo

          ! Read in all elmt nodes
	  Do idx=1, telmts
	    read(unit=6,fmt=*) iholder, iholder  ! skip another line
	    read(unit=6,fmt=*) elmt_node_list(3*idx-2), elmt_node_list(3*idx-1), &
		           elmt_node_list(3*idx) 
	  EndDo

	  ! Read in all elmt edges
	  Do idx=1, telmts 
            read(unit=6,fmt=*) iholder  ! skip another line
	    read(unit=6,fmt=*) elmt_edge_list(3*idx-2), elmt_edge_list(3*idx-1), & 
                           elmt_edge_list(3*idx)
	  EndDo

	  ! Read in all edge nodes
          Do idx=1, tedges
            read(unit=6,fmt=*) iholder  ! skip another line
            read(unit=6,fmt=*) edge_node_list(3*idx-2), edge_node_list(3*idx-1), &
                           edge_node_list(3*idx)
          EndDo


	  ! Enforce counterclockwise convention between elem_edge and elmt_node 
	  call Check_mesh(tnodes, tedges, telmts, & 
	         elmt_node_list, node_coord_list, elmt_edge_list, edge_node_list)


	  ! Mesh partition among all procs
	  nelmts_send = telmts/totalProcs
          nelmts_local = telmts - (totalProcs - 1)*nelmts_send

	  allocate(elmt_ngb_list(3*(nelmts_local+1)), stat=err_code)
	  allocate(tmp_array(3*(nelmts_local+1)), stat=err_code)	

	  if(debug_) then
	    print*,'nelmts_local = ', nelmts_local
	    print*,'nelmts_send = ', nelmts_send
	  endif

	  ! Send mesh info to other procs
	  Do idx = 1, totalProcs-1
	    istart = 3*nelmts_local + 3*(idx-1)*nelmts_send

	    ! Compute elmt_ngb_list
	    elmt_ngb_list(1:3*nelmts_send) = -1
	    kk = 0
	    Do jdx = 1, nelmts_send
	      Do kdx = 1, telmts
	        if(nelmts_local+(idx-1)*nelmts_send+jdx .ne. kdx) then
	          if(elmt_edge_list(istart+3*jdx-2) == &
		     elmt_edge_list(3*kdx-2) .or. &
		     elmt_edge_list(istart+3*jdx-2) == &
		     elmt_edge_list(3*kdx-1) .or. &
		     elmt_edge_list(istart+3*jdx-2) == &
		     elmt_edge_list(3*kdx)) then
		    elmt_ngb_list(3*jdx-2) = kdx; kk = kk + 1
		  endif 
	 	  if(elmt_edge_list(istart+3*jdx-1) == &
		     elmt_edge_list(3*kdx-2) .or. &
                     elmt_edge_list(istart+3*jdx-1) == &
		     elmt_edge_list(3*kdx-1) .or. & 
                     elmt_edge_list(istart+3*jdx-1) == &
		     elmt_edge_list(3*kdx)) then
                    elmt_ngb_list(3*jdx-1) = kdx; kk = kk + 1
                  endif	
		  if(elmt_edge_list(istart+3*jdx) == &
		     elmt_edge_list(3*kdx-2) .or. &    
                     elmt_edge_list(istart+3*jdx) == &
		     elmt_edge_list(3*kdx-1) .or. &  
                     elmt_edge_list(istart+3*jdx) == &
		     elmt_edge_list(3*kdx)) then     
                    elmt_ngb_list(3*jdx) = kdx; kk = kk + 1  
                  endif
	        endif
	      EndDo
	    EndDo

	    !
	    ! Send sub-mesh to proc idx 
	    !
            msg_tag = 0; ibuf(1) = nelmts_send; ibuf(2) = tnodes
            call Msg_iSend_int(ibuf, 2, idx, msg_tag, hdl(1)) 

            msg_tag = msg_tag + 1   
	    tmp_array(1:3*nelmts_send) = &
              elmt_node_list(istart+1:istart+3*nelmts_send)
            call Msg_iSend_int(tmp_array, 3*nelmts_send, idx, msg_tag, hdl(2))
            msg_tag = msg_tag + 1
	    tmp_array(1:3*nelmts_send) = &
	      elmt_edge_list(istart+1:istart+3*nelmts_send)
	    call Msg_iSend_int(tmp_array, 3*nelmts_send, idx, msg_tag, hdl(3)) 
            msg_tag = msg_tag + 1
	    call Msg_iSend_int(elmt_ngb_list, 3*nelmts_send, idx, msg_tag, hdl(4)) 
	    msg_tag = msg_tag + 1
	    call Msg_iSend_real(node_coord_list, 2*tnodes, idx, msg_tag, hdl(5))
	    !
	    call Msg_wait(hdl, 5)
          EndDo
	
          kk = 0 
	  elmt_ngb_list(1:3*nelmts_local) = -1 
          Do jdx = 1, nelmts_local
            Do kdx = 1, telmts   
              if(jdx .ne. kdx) then 
                if(elmt_edge_list(3*jdx-2) == elmt_edge_list(3*kdx-2) .or. &
                   elmt_edge_list(3*jdx-2) == elmt_edge_list(3*kdx-1) .or. &
                   elmt_edge_list(3*jdx-2) == elmt_edge_list(3*kdx)) then
                  elmt_ngb_list(3*jdx-2) = kdx; kk = kk + 1
                endif
                if(elmt_edge_list(3*jdx-1) == elmt_edge_list(3*kdx-2) .or. &
                   elmt_edge_list(3*jdx-1) == elmt_edge_list(3*kdx-1) .or. &
                   elmt_edge_list(3*jdx-1) == elmt_edge_list(3*kdx)) then
                  elmt_ngb_list(3*jdx-1) = kdx; kk = kk + 1
                endif
                if(elmt_edge_list(3*jdx) == elmt_edge_list(3*kdx-2) .or. &
                   elmt_edge_list(3*jdx) == elmt_edge_list(3*kdx-1) .or. &
                   elmt_edge_list(3*jdx) == elmt_edge_list(3*kdx)) then
                  elmt_ngb_list(3*jdx) = kdx; kk = kk + 1
                endif
              endif
            EndDo 
          EndDo

	  loc_elmts = nelmts_local 
 100	  format(A)

	else 
  	  ! Procs other than 0
	  ! Receive sub-mesh from proc 0
	  msg_tag = 0 
          call Msg_iRecv_int(ibuf, 2, 0, msg_tag, hdl(1))
	  call Msg_wait(hdl, 1)
	  nelmts_recvd = ibuf(1); tnodes = ibuf(2)
          msg_tag = msg_tag + 1 
	
	  allocate(elmt_node_list(3*nelmts_recvd), stat=err_code) 
	  allocate(elmt_edge_list(3*nelmts_recvd), stat=err_code) 
	  allocate(elmt_ngb_list(3*nelmts_recvd), stat=err_code) 
	  allocate(node_coord_list(2*tnodes), stat=err_code)

	  call Msg_iRecv_int(elmt_node_list, 3*nelmts_recvd, 0, &
			     msg_tag, hdl(1)) 
          msg_tag = msg_tag + 1
	  call Msg_iRecv_int(elmt_edge_list, 3*nelmts_recvd, 0, &
                             msg_tag, hdl(2)) 
          msg_tag = msg_tag + 1
	  call Msg_iRecv_int(elmt_ngb_list, 3*nelmts_recvd, 0, & 
                             msg_tag, hdl(3))  
          msg_tag = msg_tag + 1
	  call Msg_iRecv_real(node_coord_list, 2*tnodes, 0, msg_tag, hdl(4))

	  call Msg_wait(hdl, 4)

	  loc_elmts = nelmts_recvd
	endif
	End Subroutine Read_dist_mesh


	Subroutine Build_mesh_struct(in_mesh, loc_elmts, &
			 elmt_node_array, elmt_edge_array, &
		         node_coord_array, elmt_neigh, &
		         my_pid, tot_procs)
        Use GlobalVars_Mod 
        Use Mesh_Mod 
        implicit none
          
        Type(Mesh) in_mesh 
	integer, intent (in) :: loc_elmts, my_pid, tot_procs
        integer, dimension(:), pointer :: elmt_node_array, &
                               elmt_edge_array, elmt_neigh 
        real, dimension(:), pointer :: node_coord_array
	integer, dimension(:), pointer :: tmp_nodes, tmp_edges
	integer, dimension(:), pointer :: sum_array
	integer :: idx, jdx, offset
	logical :: found

	in_mesh%loc_elements = loc_elmts
	allocate(tmp_nodes(ele_nodes*in_mesh%loc_elements), stat=err_code)	
	allocate(tmp_edges(ele_edges*in_mesh%loc_elements), stat=err_code)	
	allocate(sum_array(tot_procs), stat=err_code)

	! Compute offset for global element id
	sum_array = 0
        if(myPid < tot_procs-1) then
	  sum_array(my_pid+1:tot_procs) = in_mesh%loc_elements
	endif
	call MPI_ALLREDUCE(sum_array, sum_array, tot_procs, MPI_INTEGER, &
			   MPI_SUM, MPI_COMM_WORLD, err_code)

	in_mesh%loc_nodes = 0; in_mesh%loc_edges = 0
	!
	! Build node arrays
	!
	Do idx = 1, ele_nodes*in_mesh%loc_elements
          found = .false.
          Do jdx = 1, idx-1
            if(elmt_node_array(idx) == elmt_node_array(jdx)) then
	       tmp_nodes(idx) = tmp_nodes(jdx) ! local index
              found = .true.; exit
            endif
          EndDo
          if(.not. found) then
	    in_mesh%loc_nodes = in_mesh%loc_nodes + 1
	    tmp_nodes(idx) = in_mesh%loc_nodes ! local indx
	  endif  
        EndDo

	!
        ! Build edge arrays
        !          
        Do idx = 1, ele_edges*in_mesh%loc_elements
          found = .false.
          Do jdx = 1, idx-1
            if(elmt_edge_array(idx) == elmt_edge_array(jdx)) then
	      tmp_edges(idx) = tmp_edges(jdx) ! local index
              found = .true.; exit
            endif
          EndDo
          if(.not. found) then
	    in_mesh%loc_edges = in_mesh%loc_edges + 1
	    tmp_edges(idx) = in_mesh%loc_edges ! local indx
	  endif
        EndDo

        ! 
        ! Create local mesh structure 
	!
        in_mesh = Alloc_Mesh(in_mesh%loc_elements, &
			     in_mesh%loc_edges, in_mesh%loc_nodes)

	! Assign local ids in edge and element structures
  	Do idx = 1, in_mesh%loc_elements
	  in_mesh%elements(idx)%node_indx(1:3) = tmp_nodes(3*idx-2:3*idx) 
	  in_mesh%elements(idx)%edge_indx(1:3) = tmp_edges(3*idx-2:3*idx) 
        EndDo
	Do idx = 1, in_mesh%loc_elements
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(1))%node_indx(1) &
	    = in_mesh%elements(idx)%node_indx(1)
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(1))%node_indx(2) &
	    = in_mesh%elements(idx)%node_indx(2)
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(2))%node_indx(1) & 
            = in_mesh%elements(idx)%node_indx(2) 
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(2))%node_indx(2) &
            = in_mesh%elements(idx)%node_indx(3)
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(3))%node_indx(1) & 
            = in_mesh%elements(idx)%node_indx(3) 
          in_mesh%edges(in_mesh%elements(idx)%edge_indx(3))%node_indx(2) &
            = in_mesh%elements(idx)%node_indx(1)
	EndDo

	! Assign global ids for nodes and edges
	Do idx = 1, in_mesh%loc_elements
	  in_mesh%nodes(in_mesh%elements(idx)%node_indx(1))%global_id = &
	    	   elmt_node_array(3*idx-2)
	  in_mesh%nodes(in_mesh%elements(idx)%node_indx(2))%global_id = &
		   elmt_node_array(3*idx-1)
	  in_mesh%nodes(in_mesh%elements(idx)%node_indx(3))%global_id = &
		   elmt_node_array(3*idx)
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(1))%global_id = & 
                   elmt_edge_array(3*idx-2)
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(2))%global_id = & 
                   elmt_edge_array(3*idx-1)
	  in_mesh%edges(in_mesh%elements(idx)%edge_indx(3))%global_id = & 
                   elmt_edge_array(3*idx)
	EndDo
	
	! Assign global id for elements
        if(my_pid == 0) then
	  offset = 0
	else 
	  offset = sum_array(my_pid)
	endif
	!
	Do idx = 1, in_mesh%loc_elements
	  in_mesh%elements(idx)%global_id = offset + 1
	  offset = offset + 1
	EndDo 
	
	! Build node coord 
	Do idx = 1, in_mesh%loc_nodes
	  offset = in_mesh%nodes(idx)%global_id
	  in_mesh%nodes(idx)%coord(1:2) = node_coord_array(2*offset-1:2*offset)
	EndDo

	Do idx = 1, in_mesh%loc_elements
          in_mesh%elements(idx)%neigh_indx(1:3) = elmt_neigh(3*idx-2:3*idx)
        EndDo

        deallocate(elmt_node_array, elmt_edge_array, node_coord_array, &
		   elmt_neigh)
	!
        End Subroutine Build_Mesh_Struct


	Subroutine Write_distrMesh(trimesh, filename)
	Type(Mesh) trimesh
	character(len=*) :: filename
	Type(Element), dimension(:), pointer :: elmt_array
        Type(Node), dimension(:), pointer :: node_array
	integer :: funit, idx, tag
	integer :: msg(1)

	call Global_Barrier()

	elmt_array => trimesh%elements
        node_array => trimesh%nodes
	funit = 77; msg(1) = 0; tag = 0

	if(debug_) &
	  print*,'In proc ', myPid, 'elmts, edges, nodes = ', &
            trimesh%loc_elements, trimesh%loc_edges, trimesh%loc_nodes

	trimesh%loc_edges = trimesh%loc_edges
	trimesh%loc_nodes = trimesh%loc_nodes

	if(myPid == 0) then
	  open(Unit=funit, File=filename, Status="replace")
	  Do idx = 1, trimesh%loc_elements
            write(funit,300) node_array(elmt_array(idx)%node_indx(1))%coord(1), &
                             node_array(elmt_array(idx)%node_indx(1))%coord(2), &
                             node_array(elmt_array(idx)%node_indx(2))%coord(1), &
                             node_array(elmt_array(idx)%node_indx(2))%coord(2), &
                             node_array(elmt_array(idx)%node_indx(3))%coord(1), &
                             node_array(elmt_array(idx)%node_indx(3))%coord(2)
          EndDo	  
	  close(funit)
	  call Msg_Send_int(msg, 1, myPid+1, tag)
	else
	  call Msg_Recv_int(msg, 1, myPid-1, tag)
	  open(Unit=funit, File="mesh_data", Status="old", Position="APPEND")
	  Do idx=1, trimesh%loc_elements           
	    write(funit,300) node_array(elmt_array(idx)%node_indx(1))%coord(1), &
                             node_array(elmt_array(idx)%node_indx(1))%coord(2), &
                             node_array(elmt_array(idx)%node_indx(2))%coord(1), &
                             node_array(elmt_array(idx)%node_indx(2))%coord(2), &
                             node_array(elmt_array(idx)%node_indx(3))%coord(1), &
                             node_array(elmt_array(idx)%node_indx(3))%coord(2)
          EndDo 
	  close(funit)
	  if(myPid < totalProcs-1) call Msg_Send_int(msg, 1, myPid+1, tag) 
	endif 
 300    format(6(f8.3, 2x)) 
	!
	End Subroutine Write_distrMesh


	! Check counterclockwise property for the mesh
        Subroutine Check_mesh(tot_nodes, tot_edges, tot_elmts, &
                     elmt_node_list, node_coord_list, elmt_edge_list, edge_nodes)
        implicit none
        integer :: tot_nodes, tot_edges, tot_elmts
        integer :: elmt_node_list(3,tot_elmts)
        integer :: elmt_edge_list(3,tot_elmts)
        real :: node_coord_list(2,tot_nodes)
        integer :: edge_nodes(3,tot_edges)
        integer :: idx, itmp, e(3)
        real :: v1(2), v2(2), cp

        Do idx = 1, tot_elmts
          v1(1) = node_coord_list(1,elmt_node_list(2,idx)) - &
                    node_coord_list(1,elmt_node_list(1,idx))
          v1(2) = node_coord_list(2,elmt_node_list(2,idx)) - &
                    node_coord_list(2,elmt_node_list(1,idx))
          v2(1) = node_coord_list(1,elmt_node_list(3,idx)) - &
                    node_coord_list(1,elmt_node_list(1,idx))
          v2(2) = node_coord_list(2,elmt_node_list(3,idx)) - &
                    node_coord_list(2,elmt_node_list(1,idx))
          cp = v1(1)*v2(2) - v1(2)*v2(1)
          if(cp < 0.0) then
            itmp = elmt_node_list(3,idx)
            elmt_node_list(3,idx) = elmt_node_list(1,idx)
            elmt_node_list(1,idx) = itmp
          endif
        EndDo

        Do idx = 1, tot_elmts
          e(1) = elmt_edge_list(1,idx)
          e(2) = elmt_edge_list(2,idx)
          e(3) = elmt_edge_list(3,idx)
          if(edge_nodes(1, e(1)) == elmt_node_list(2,idx) .and. &
             edge_nodes(2, e(1)) == elmt_node_list(3,idx) &
             .or. &
             edge_nodes(2, e(1)) == elmt_node_list(2,idx) .and. &
             edge_nodes(1, e(1)) == elmt_node_list(3,idx)) then
             elmt_edge_list(2,idx) = e(1)
          else if(edge_nodes(1, e(1)) == elmt_node_list(3,idx) .and. &
                  edge_nodes(2, e(1)) == elmt_node_list(1,idx) &
                  .or. &
                  edge_nodes(2, e(1)) == elmt_node_list(3,idx) .and. &
                  edge_nodes(1, e(1)) == elmt_node_list(1,idx)) then
                  elmt_edge_list(3,idx) = e(1)
          endif

	  if(edge_nodes(1, e(2)) == elmt_node_list(1,idx) .and. &
             edge_nodes(2, e(2)) == elmt_node_list(2,idx) &
             .or. &
             edge_nodes(2, e(2)) == elmt_node_list(1,idx) .and. &
             edge_nodes(1, e(2)) == elmt_node_list(2,idx)) then
             elmt_edge_list(1,idx) = e(2)
          else if(edge_nodes(1, e(2)) == elmt_node_list(3,idx) .and. &
                  edge_nodes(2, e(2)) == elmt_node_list(1,idx) &
                  .or. &
                  edge_nodes(2, e(2)) == elmt_node_list(3,idx) .and. &
                  edge_nodes(1, e(2)) == elmt_node_list(1,idx)) then
                  elmt_edge_list(3,idx) = e(2)
          endif
             
          if(edge_nodes(1, e(3)) == elmt_node_list(1,idx) .and. &
             edge_nodes(2, e(3)) == elmt_node_list(2,idx) &
             .or. &
             edge_nodes(2, e(3)) == elmt_node_list(1,idx) .and. &
             edge_nodes(1, e(3)) == elmt_node_list(2,idx)) then
             elmt_edge_list(1,idx) = e(3)
          else if(edge_nodes(1, e(3)) == elmt_node_list(2,idx) .and. &
                  edge_nodes(2, e(3)) == elmt_node_list(3,idx) &
                  .or. &
                  edge_nodes(2, e(3)) == elmt_node_list(2,idx) .and. &
                  edge_nodes(1, e(3)) == elmt_node_list(3,idx)) then
                  elmt_edge_list(2,idx) = e(3)
          endif
        EndDo
        End Subroutine Check_mesh


	! ---------------------------------------------------------------
        ! Mesh Visualize
	!
        ! Create a visualization file using TecPlot format
	! Create zones for each processor - very fast for tecplot
	! ---------------------------------------------------------------
         SUBROUTINE PAMR_VISUALIZE(this, mesh_file)
           IMPLICIT NONE
           type (mesh), intent(in) :: this
           character(len=*), intent(in), optional :: mesh_file

           integer, parameter :: OS_FILE = 6
           integer :: i, iter, ierror, ios
           character(len=32) :: in_elems, in_nodes

	! ---------------------------------------------------------------
	! Processor 0 creates the initial file
	! ---------------------------------------------------------------

           if ( myPid == 0 ) then

             if (present(mesh_file)) then
                 OPEN(UNIT=OS_FILE, FILE=mesh_file, STATUS="REPLACE",   &
     &                ACTION="WRITE", IOSTAT=ios)
             else
                 OPEN(UNIT=OS_FILE, FILE="pamr_mesh.plt",                &
     &                STATUS="REPLACE", ACTION="WRITE", IOSTAT=ios)
             end if

             if ( ios /= 0 ) then
                  print *,"mesh_visualize: Can't create mesh" ; stop
             end if
 
             WRITE(UNIT=OS_FILE, FMT='(a)') 'VARIABLES = "X", "Y"'
 
             WRITE(in_elems,'(i32)') this%loc_elements
             WRITE(in_nodes,'(i32)') this%loc_nodes
 
             WRITE(UNIT=OS_FILE, FMT='(a,a,a,a,a)')                     &
     &             "ZONE N=",trim(adjustl(in_nodes)),                   &
     &             ", E=",trim(adjustl(in_elems)),                      &
     &             ", F=FEPOINT, ET=TRIANGLE"
 
             do i = 1, this%loc_nodes
                WRITE(UNIT=OS_FILE, FMT='(F6.2,F6.2)')                  &
     &               this%nodes(i)%coord(1), this%nodes(i)%coord(2)
             end do

	     WRITE(UNIT=OS_FILE, FMT=*)

                  do i = 1, this%loc_elements
                      WRITE(UNIT=OS_FILE, FMT='(I6,I6,I6)')             &
     &                      this%elements(i)%node_indx(1),              &
     &                      this%elements(i)%node_indx(2),              &
     &                      this%elements(i)%node_indx(3)
                  end do

                  CLOSE(UNIT=OS_FILE, STATUS="KEEP", IOSTAT=ios)
 
                end if
 
                call MPI_BARRIER(MPI_COMM_WORLD, ierror)

	   ! ---------------------------------------------------------
           ! Remaining processors append to existing file
           ! ---------------------------------------------------------
  
           do iter = 2, totalProcs 
 
                if ( iter == (myPid+1) ) then
 
                  if (present(mesh_file)) then
                     OPEN(UNIT=OS_FILE, FILE=mesh_file, STATUS="OLD",   &
     &                    POSITION="APPEND", ACTION="WRITE", IOSTAT=ios)
                  else
                     OPEN(UNIT=OS_FILE,FILE="pamr_mesh.plt",             &
     &                    STATUS="OLD", POSITION="APPEND",              &
     &                    ACTION="WRITE", IOSTAT=ios)
                  end if
 
                  if ( ios /= 0 ) then
                       print *,"mesh_visualize: Can't create mesh"
                       stop
                  end if
 
                  WRITE(in_elems,'(i32)') this%loc_elements
                  WRITE(in_nodes,'(i32)') this%loc_nodes
 
                  WRITE(UNIT=OS_FILE, FMT=*)
 
                  WRITE(UNIT=OS_FILE, FMT='(a,a,a,a,a)')                &
     &                  "ZONE N=",trim(adjustl(in_nodes)),              &
     &                  ", E=",trim(adjustl(in_elems)),                 &
     &                  ", F=FEPOINT, ET=TRIANGLE"
 
                  do i = 1, this%loc_nodes
                     WRITE(UNIT=OS_FILE, FMT='(F6.2,F6.2)')             &
     &                    this%nodes(i)%coord(1), this%nodes(i)%coord(2)
                  end do
 
                  WRITE(UNIT=OS_FILE, FMT=*)
 
                  do i = 1, this%loc_elements
                      WRITE(UNIT=OS_FILE, FMT='(I6,I6,I6)')             &
     &                      this%elements(i)%node_indx(1),              &
     &                      this%elements(i)%node_indx(2),              &
     &                      this%elements(i)%node_indx(3)
                  end do
 
                  CLOSE(UNIT=OS_FILE, STATUS="KEEP", IOSTAT=ios)
 
                end if
 
                call MPI_BARRIER(MPI_COMM_WORLD, ierror)
 
           end do     
        END SUBROUTINE PAMR_VISUALIZE


	integer Function Insert_item(list, item_key, n_crnt_items)
        Type(r_indx), dimension(:), pointer :: list
        integer :: item_key, n_crnt_items
        integer :: idx
        !
        if(n_crnt_items > 0) then
          Do idx = 1, n_crnt_items
            if(list(idx)%global_id == item_key) then
              Insert_item = idx; return
            endif
          EndDo
        endif
        n_crnt_items = n_crnt_items + 1
        Insert_item = n_crnt_items
        list(n_crnt_items)%global_id = item_key
        End Function Insert_item


	integer Function InsertToList(list, item_key, n_crnt_items)
        integer, dimension(:), pointer :: list
        integer :: item_key, n_crnt_items
        integer :: idx
        !
        if(n_crnt_items > 0) then
          Do idx = 1, n_crnt_items
            if(list(idx) == item_key) then
              InsertToList = idx; return
            endif
          EndDo
        endif
        n_crnt_items = n_crnt_items + 1
        InsertToList = n_crnt_items
        list(n_crnt_items) = item_key
        End Function InsertToList
	!
        End Module Util_Mod
