subroutine grid_output(this, filename, other)

    class(mesh) :: this
    class(cell), optional, intent(in) :: other
    character(len=20), intent(in) :: filename
    integer :: tag, ip, L(3), nodes, s(3), icom(3), ip_new, gcom(3), lp(3), i

    tag = 33423
    L = this%L
    nodes = this%nodes
    if(present(other)) then
      gcom = (/ int(anint( other%gcom(1) - other%L(1)/2)), int(anint( other%gcom(2) - other%L(2)/2)), int(anint(other%gcom(3) - other%L(3)/2)) /)
      call check_boundary(gcom(1),this%L(1),this%b)
      call check_boundary(gcom(2),this%L(2),this%b)
      call check_boundary(gcom(3),this%L(3),this%b)
    end if

    OPEN(UNIT=tag, FILE=trim(filename)//".vti" )
    write(tag,'(A)')'<?xml version="1.0"?>'
    write(tag,'(A)')'<VTKFile type="ImageData" version="0.1" byte_order="LittleEndian">'
    write(tag,'(A,I1,A,I3,A,I1,A,I3,A,I1,A,I3,A)')'  <ImageData WholeExtent="',0,' ',L(1),&
    ' ',0,' ',L(2),' ',0,' ',L(3),'"   Origin="0 0 0" Spacing="1 1 1">'
    write(tag,'(A,I1,A,I3,A,I1,A,I3,A,I1,A,I3,A)')'    <Piece Extent="',0,' ',L(1),' ',0,' ',L(2),' ',0,' ',L(3),'">'
    write(tag,*)'      <CellData> '
    write(tag,*)'        <DataArray Name="scalar_data" type="Float64" format="ascii">'
    do ip=0, nodes-1
      if(present(other)) then
        call vec_global2local(s, this%position(ip), gcom, this%L)
        ip_new = other%ip(s)
        write(tag,'(F10.2)', ADVANCE='no') other%gt(ip_new)
      else
        write(tag,'(F10.2)', ADVANCE='no') this%gt(ip)
      end if
    end do
    write(tag,'(A)') ""
    write(tag,'(A)')"         </DataArray>"
    write(tag,'(A)')"      </CellData>"
    write(tag,'(A)')"    </Piece>"
    write(tag,'(A)')"</ImageData>"
    write(tag,'(A)')"</VTKFile>"
    CLOSE(tag)
  end subroutine grid_output