diff --git a/src_output/mapVTKOutput.F90 b/src_output/mapVTKOutput.F90
new file mode 100644
index 00000000..95d7b0f7
--- /dev/null
+++ b/src_output/mapVTKOutput.F90
@@ -0,0 +1,343 @@
+module mapVTKOutput_m
+   use FDETYPES_m
+   use outputTypes_m
+   use outputUtils_m
+   use directoryUtils_m
+   use allocationUtils_m
+   use vtkAPI_m
+   use volumicProbeUtils_m
+   use report_m
+
+   implicit none
+contains
+   subroutine init_mapvtk_output(this, lowerBound, upperBound, field, outputTypeExtension, mpidir, problemInfo)
+      type(mapvtk_output_t), intent(out) :: this
+      type(cell_coordinate_t), intent(in) :: lowerBound, upperBound
+      type(problem_info_t), target ,intent(in) :: problemInfo
+      integer(kind=SINGLE), intent(in) :: mpidir, field
+      character(len=BUFSIZE), intent(in) :: outputTypeExtension
+
+      integer(kind=SINGLE) :: i
+
+      this%mainCoords = lowerBound
+      this%auxCoords = upperBound
+      this%component = field
+
+      this%path = get_output_path()
+      call store_relevant_coordinates(this, problemInfo)
+
+   contains
+
+      function get_output_path() result(outputPath)
+         character(len=BUFSIZE)  :: probeBoundsExtension, prefixFieldExtension
+         character(len=BUFSIZE) :: outputPath
+         probeBoundsExtension = get_coordinates_extension(this%mainCoords, this%auxCoords, mpidir)
+         prefixFieldExtension = get_prefix_extension(field, mpidir)
+         outputPath = &
+            trim(adjustl(outputTypeExtension))//'_'//trim(adjustl(prefixFieldExtension))//'_'//trim(adjustl(probeBoundsExtension))
+         return
+      end function
+
+      subroutine store_media_tag()
+      end subroutine store_media_tag
+   end subroutine init_mapvtk_output
+
+   subroutine store_relevant_coordinates(this, problemInfo)
+      type(mapvtk_output_t), intent(inout) :: this
+      type(problem_info_t), pointer, intent(in) :: problemInfo
+
+      integer :: i, j, k, field, counter
+
+      counter = 0
+      do k = this%mainCoords%Z, this%auxCoords%Z
+      do j = this%mainCoords%Y, this%auxCoords%Y
+      do i = this%mainCoords%X, this%auxCoords%X
+         do field = iEx, iEz
+            if (isEdge(field, i, j, k, problemInfo)) then
+               counter = counter + 1
+            end if
+         end do
+         do field = iHx, iHz
+            if (isWithinBounds(field, i, j, k, problemInfo)) then
+               if (isMaterialExceptPML(field, i, j, k, problemInfo)) then
+                  counter = counter + 1
+               end if
+               if (problemInfo%materialTag%getFaceTag(field, i, j, k) < 0 &
+                   .and. (btest(iabs(problemInfo%materialTag%getFaceTag(field, i, j, k)), field - 1)) &
+                   .and. (.not. isPML(field, i, j, k, problemInfo))) then
+                  counter = counter + 1
+               end if
+            end if
+         end do
+      end do
+      end do
+      end do
+
+      this%nPoints = counter
+      call alloc_and_init(this%coords, 3, this%nPoints, -99)
+      call alloc_and_init(this%materialTag, this%nPoints, -1)
+      call alloc_and_init(this%currentType, this%nPoints, -1)
+
+      counter = 0
+      do k = this%mainCoords%Z, this%auxCoords%Z
+      do j = this%mainCoords%Y, this%auxCoords%Y
+      do i = this%mainCoords%X, this%auxCoords%X
+         do field = iEx, iEz
+            if (isEdge(field, i, j, k, problemInfo)) then
+               counter = counter + 1
+               call writeFaceTagInfo(this, counter, i, j, k, field, problemInfo%materialTag%getFaceTag(field, i, j, k))
+            end if
+         end do
+         do field = iHx, iHz
+            if (isWithinBounds(field, i, j, k, problemInfo)) then
+               if (isMaterialExceptPML(field, i, j, k, problemInfo)) then
+                  counter = counter + 1
+                  call writeFaceTagInfo(this, counter, i, j, k, field, problemInfo%materialTag%getFaceTag(field, i, j, k))
+               end if
+               if (problemInfo%materialTag%getFaceTag(field, i, j, k) < 0 &
+                   .and. (btest(iabs(problemInfo%materialTag%getFaceTag(field, i, j, k)), field - 1)) &
+                   .and. .not. isPML(field, i, j, k, problemInfo)) then
+                  counter = counter + 1
+                  call writeFaceTagInfo(this, counter, i, j, k, field, problemInfo%materialTag%getFaceTag(field, i, j, k))
+               end if
+            end if
+         end do
+      end do
+      end do
+      end do
+   end subroutine store_relevant_coordinates
+
+   logical function isMaterialExceptPML(field, i, j, k, problemInfo)
+      integer, intent(in)             :: field, i, j, k
+      type(problem_info_t), intent(in):: problemInfo
+      isMaterialExceptPML = .not. isMediaVacuum(field, i, j, k, problemInfo)
+      isMaterialExceptPML = isMaterialExceptPML .and. (.not. isPML(field, i, j, k, problemInfo))
+   end function isMaterialExceptPML
+
+   subroutine writeFaceTagInfo(this, counter, i, j, k, field, tag)
+      type(mapvtk_output_t), intent(inout) :: this
+      integer, intent(in) :: i, j, k, counter, tag, field
+      this%coords(1, counter) = i
+      this%coords(2, counter) = j
+      this%coords(3, counter) = k
+      this%currentType(counter) = currentType(field)
+      this%materialTag(counter) = tag
+   end subroutine
+
+   subroutine create_geometry_simulation_vtu(this, control, realXGrid, realYGrid, realZGrid)
+      implicit none
+
+      type(mapvtk_output_t), intent(in) :: this
+      type(sim_control_t), intent(in) :: control
+      real(KIND=RKIND), pointer, dimension(:), intent(in) :: realXGrid, realYGrid, realZGrid
+
+      !type(vtk_file) :: vtkOutput
+      type(vtk_unstructured_grid), target :: ugrid
+
+      integer :: ierr, i, npts, unit
+      real(RKIND), allocatable :: x(:), y(:), z(:), materialTag(:)
+      character(len=BUFSIZE) :: info_str
+      character(len=BUFSIZE) :: metadata_filename, vtuPath
+
+      integer, allocatable :: conn(:), offsets(:), types(:)
+      integer :: numNodes, numEdges, numQuads
+      real(kind=RKIND), allocatable :: nodes(:, :)
+      integer(kind=SINGLE), allocatable :: edges(:, :), quads(:, :)
+
+      call create_folder(this%path, ierr)
+      vtuPath = join_path(this%path, get_last_component(this%path))//vtuFileExtension
+
+      call createUnstructuredDataForVTU(this%nPoints, this%coords, this%currentType, nodes, edges, quads, numNodes, numEdges, numQuads, control%vtkindex, realXGrid, realYGrid, realZGrid)
+      call ugrid%add_points(real(nodes, 4))
+
+      allocate (conn(2*numEdges + 4*numQuads))
+      conn(1:2*numEdges) = reshape(edges, [2*numEdges])
+      conn(2*numEdges + 1:2*numEdges + 4*numQuads) = reshape(quads, [4*numQuads])
+
+      allocate (offsets(numEdges + numQuads))
+      do i = 1, numEdges + numQuads
+         if (i <= numEdges) then
+            if (i == 1) then
+               offsets(i) = 2
+            else
+               offsets(i) = offsets(i - 1) + 2
+            end if
+         else
+            if (i == 1) then
+               offsets(i) = 4
+            else
+               offsets(i) = offsets(i - 1) + 4
+            end if
+         end if
+      end do
+
+      allocate (types(numEdges + numQuads))
+      types(1:numEdges) = 3
+      types(numEdges + 1:numEdges + numQuads) = 9
+
+      call ugrid%add_cell_connectivity(conn, offsets, types)
+      if (size(offsets) /= numQuads) then
+         print *, "Problema con offsets"
+      end if
+      if (size(types) /= numQuads) then
+         print *, "Problema con types"
+      end if
+      if (offsets(numQuads) /= size(conn)) then
+         print *, "Tenemos un problema con conn y offset"
+      end if
+
+      call ugrid%write_file(vtuPath)
+
+      !---------------------------------------------
+      ! Metadata: write to .txt file
+      !---------------------------------------------
+      info_str = 'PEC=0, already_YEEadvanced_byconformal=5, NOTOUCHNOUSE=6, '// &
+                 'WIRE=7, WIRE-COLISION=8, COMPO=3, DISPER=1, DIEL=2, SLOT=4, '// &
+                 'CONF=5/6, OTHER=-1 (ADD +0.5 for borders)'
+
+      metadata_filename = trim(this%path)//'.txt'
+      open (newunit=unit, file=metadata_filename, status='replace', action='write', iostat=ierr)
+      if (ierr /= 0) then
+         print *, 'Error opening metadata file: ', metadata_filename
+      else
+         write (unit, '(A)') trim(info_str)
+         close (unit)
+      end if
+
+   end subroutine create_geometry_simulation_vtu
+
+   logical function isEdge(campo, iii, jjj, kkk, problemInfo)
+      integer(4), intent(in) :: campo, iii, jjj, kkk
+      type(problem_info_t), pointer, intent(in) :: problemInfo
+
+      type(MediaData_t), pointer, dimension(:) :: mData
+      type(limit_t), pointer, dimension(:) :: problemDimension
+
+      integer(4) :: imed, imed1, imed2, imed3, imed4, contaborde
+
+      mData => problemInfo%materialList
+      problemDimension => problemInfo%problemDimension
+      isEdge = .false.
+      contaborde = 0
+
+      call get_media_from_coord_and_h_neighbours(campo, iii, jjj, kkk,  problemInfo%geometryToMaterialData, imed, imed1, imed2, imed3, imed4)
+
+      if (imed /= 1) then
+
+         if (mData(imed)%is%SGBC) then
+
+            if (mData(imed1)%is%SGBC) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed1)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed1 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed2)%is%SGBC) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed2)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed2 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed3)%is%SGBC) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed3)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed3 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed4)%is%SGBC) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed4)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed4 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+         elseif (mData(imed)%is%Multiport) then
+
+            if (mData(imed1)%is%Multiport) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed1)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed1 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed2)%is%Multiport) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed2)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed2 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed3)%is%Multiport) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed3)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed3 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed4)%is%Multiport) then
+               if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed4)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed4 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+         elseif (mData(imed)%is%AnisMultiport) then
+
+            if (mData(imed1)%is%AnisMultiport) then
+               if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed1)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed1 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed2)%is%AnisMultiport) then
+               if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed2)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed2 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed3)%is%AnisMultiport) then
+               if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed3)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed3 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+            if (mData(imed4)%is%AnisMultiport) then
+               if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= &
+                   trim(adjustl(mData(imed4)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1
+            elseif (imed4 /= 1) then
+               contaborde = contaborde + 1
+            end if
+
+         else
+            if ((imed /= imed1) .and. (imed1 /= 1)) contaborde = contaborde + 1
+            if ((imed /= imed2) .and. (imed2 /= 1)) contaborde = contaborde + 1
+            if ((imed /= imed3) .and. (imed3 /= 1)) contaborde = contaborde + 1
+            if ((imed /= imed4) .and. (imed4 /= 1)) contaborde = contaborde + 1
+         end if
+
+         if ((imed1 == 1 .and. imed2 == 1 .and. imed3 == 1 .and. imed4 /= 1) .or. &
+             (imed2 == 1 .and. imed3 == 1 .and. imed4 == 1 .and. imed1 /= 1) .or. &
+             (imed3 == 1 .and. imed4 == 1 .and. imed1 == 1 .and. imed2 /= 1) .or. &
+             (imed4 == 1 .and. imed1 == 1 .and. imed2 == 1 .and. imed3 /= 1) .or. &
+             (imed1 == 1 .and. imed2 == 1 .and. imed3 == 1 .and. imed4 == 1) .or. &
+             (contaborde > 0)) isEdge = .true.
+
+         if ((iii > problemDimension(campo)%XE) .or. (jjj > problemDimension(campo)%YE) .or. &
+             (kkk > problemDimension(campo)%ZE)) isEdge = .false.
+
+         if ((iii < problemDimension(campo)%XI) .or. (jjj < problemDimension(campo)%YI) .or. &
+             (kkk < problemDimension(campo)%ZI)) isEdge = .false.
+
+      else
+         isEdge = .false.
+      end if
+
+   end function
+
+end module mapVTKOutput_m
diff --git a/src_output/vtkAPI.F90 b/src_output/vtkAPI.F90
new file mode 100644
index 00000000..bca99f0b
--- /dev/null
+++ b/src_output/vtkAPI.F90
@@ -0,0 +1,365 @@
+module vtkAPI_m
+   implicit none
+   private
+   public :: vtk_data_array, vtk_grid, vtk_structured_grid, vtk_unstructured_grid
+
+   !==========================
+   ! Data array type
+   !==========================
+   type :: vtk_data_array
+      character(len=:), allocatable :: name
+      character(len=:), allocatable :: type
+      integer :: num_components
+      real, allocatable :: data(:)
+   end type vtk_data_array
+
+   !==========================
+   ! Base abstract class
+   !==========================
+   type, abstract :: vtk_grid
+      real, allocatable :: points(:, :)
+      type(vtk_data_array), allocatable :: scalars(:)
+      type(vtk_data_array), allocatable :: vectors(:)
+      type(vtk_data_array), allocatable :: cell_scalars(:)
+      type(vtk_data_array), allocatable :: cell_vectors(:)
+   contains
+      procedure(add_points_generic), deferred :: add_points
+      procedure(add_scalar_generic), deferred :: add_scalar
+      procedure(add_vector_generic), deferred :: add_vector
+      procedure(add_cell_scalar_generic), deferred :: add_cell_scalar
+      procedure(add_cell_vector_generic), deferred :: add_cell_vector
+      procedure(write_file_generic), deferred :: write_file
+   end type vtk_grid
+
+   !==========================
+   ! Structured Grid (VTS)
+   !==========================
+   type, extends(vtk_grid) :: vtk_structured_grid
+      integer :: nx = 0, ny = 0, nz = 0
+   contains
+      procedure :: add_points => add_points_structured
+      procedure :: add_scalar => add_scalar_structured
+      procedure :: add_vector => add_vector_structured
+      procedure :: add_cell_scalar => add_cell_scalar_structured
+      procedure :: add_cell_vector => add_cell_vector_structured
+      procedure :: write_file => write_vts_file
+   end type vtk_structured_grid
+
+   !==========================
+   ! Unstructured Grid (VTU)
+   !==========================
+   type, extends(vtk_grid) :: vtk_unstructured_grid
+      integer :: num_points = 0
+      integer :: num_cells = 0
+      integer, allocatable :: connectivity(:)
+      integer, allocatable :: offsets(:)
+      integer, allocatable :: types(:)
+   contains
+      procedure :: add_points => add_points_unstructured
+      procedure :: add_scalar => add_scalar_unstructured
+      procedure :: add_vector => add_vector_unstructured
+      procedure :: add_cell_scalar => add_cell_scalar_unstructured
+      procedure :: add_cell_vector => add_cell_vector_unstructured
+      procedure :: add_cell_connectivity
+      procedure :: write_file => write_vtu_file
+   end type vtk_unstructured_grid
+
+   !==========================
+   ! Generic deferred interfaces
+   !==========================
+   abstract interface
+      subroutine add_points_generic(this, pts)
+         import :: vtk_grid
+         class(vtk_grid), intent(inout) :: this
+         real, intent(in) :: pts(:, :)
+      end subroutine add_points_generic
+
+      subroutine add_scalar_generic(this, name, data)
+         import :: vtk_grid
+         class(vtk_grid), intent(inout) :: this
+         character(len=*), intent(in) :: name
+         real, intent(in) :: data(:)
+      end subroutine add_scalar_generic
+
+      subroutine add_vector_generic(this, name, data)
+         import :: vtk_grid
+         class(vtk_grid), intent(inout) :: this
+         character(len=*), intent(in) :: name
+         real, intent(in) :: data(:)
+      end subroutine add_vector_generic
+
+      subroutine add_cell_scalar_generic(this, name, data)
+         import :: vtk_grid
+         class(vtk_grid), intent(inout) :: this
+         character(len=*), intent(in) :: name
+         real, intent(in) :: data(:)
+      end subroutine add_cell_scalar_generic
+
+      subroutine add_cell_vector_generic(this, name, data)
+         import :: vtk_grid
+         class(vtk_grid), intent(inout) :: this
+         character(len=*), intent(in) :: name
+         real, intent(in) :: data(:)
+      end subroutine add_cell_vector_generic
+
+      subroutine write_file_generic(this, filename)
+         import :: vtk_grid
+         class(vtk_grid), intent(in) :: this
+         character(len=*), intent(in) :: filename
+      end subroutine write_file_generic
+   end interface
+
+contains
+   !==========================
+   !==== Structured Grid Methods ====
+   !==========================
+
+   subroutine add_points_structured(this, pts)
+      class(vtk_structured_grid), intent(inout) :: this
+      real, intent(in) :: pts(:, :)
+      integer :: npts
+      npts = this%nx*this%ny*this%nz
+      if (size(pts, 1) /= 3) error stop 'add_points_structured: first dim must be 3'
+      if (size(pts, 2) /= npts) error stop 'add_points_structured: wrong number of points'
+      if (allocated(this%points)) deallocate (this%points)
+      allocate (this%points(3, npts))
+      this%points = pts
+   end subroutine add_points_structured
+
+   subroutine add_scalar_structured(this, name, data)
+      class(vtk_structured_grid), intent(inout) :: this
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      call add_array_generic(this%scalars, name, data, 1)
+   end subroutine add_scalar_structured
+
+   subroutine add_vector_structured(this, name, data)
+      class(vtk_structured_grid), intent(inout) :: this
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      call add_array_generic(this%vectors, name, data, 3)
+   end subroutine add_vector_structured
+
+   subroutine add_cell_scalar_structured(this, name, data)
+      class(vtk_structured_grid), intent(inout) :: this
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      call add_array_generic(this%cell_scalars, name, data, 1)
+   end subroutine add_cell_scalar_structured
+
+   subroutine add_cell_vector_structured(this, name, data)
+      class(vtk_structured_grid), intent(inout) :: this
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      call add_array_generic(this%cell_vectors, name, data, 3)
+   end subroutine add_cell_vector_structured
+
+   !==========================
+   ! Write VTS
+   !==========================
+   subroutine write_vts_file(this, filename)
+      class(vtk_structured_grid), intent(in) :: this
+      character(len=*), intent(in) :: filename
+      integer :: iunit
+      open (newunit=iunit, file=filename, status='replace', action='write', form='formatted')
+      write (iunit, *) '<VTKFile type="StructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">'
+      write (iunit, *) '  <StructuredGrid WholeExtent="0 ', this%nx - 1, ' 0 ', this%ny - 1, ' 0 ', this%nz - 1, '">'
+      write (iunit, *) '    <Piece Extent="0 ', this%nx - 1, ' 0 ', this%ny - 1, ' 0 ', this%nz - 1, '">'
+      call write_pointdata(iunit, this%scalars, this%vectors)
+      call write_celldata(iunit, this%cell_scalars, this%cell_vectors)
+      call write_points(iunit, this%points)
+      write (iunit, *) '    </Piece>'
+      write (iunit, *) '  </StructuredGrid>'
+      write (iunit, *) '</VTKFile>'
+      close (iunit)
+   end subroutine write_vts_file
+
+   !==========================
+   !==== Unstructured Grid Methods ====
+   !==========================
+   subroutine add_points_unstructured(this, pts)
+      real, intent(in) :: pts(:, :)
+      class(vtk_unstructured_grid), intent(inout) :: this
+      if (size(pts, 1) /= 3) error stop 'add_points_unstructured: first dim must be 3'
+      this%num_points = size(pts, 2)
+      if (allocated(this%points)) deallocate (this%points)
+      allocate (this%points(3, this%num_points))
+      this%points = pts
+   end subroutine add_points_unstructured
+
+   subroutine add_scalar_unstructured(this, name, data)
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      class(vtk_unstructured_grid), intent(inout) :: this
+      call add_array_generic(this%scalars, name, data, 1)
+   end subroutine add_scalar_unstructured
+
+   subroutine add_vector_unstructured(this, name, data)
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      class(vtk_unstructured_grid), intent(inout) :: this
+      call add_array_generic(this%vectors, name, data, 3)
+   end subroutine add_vector_unstructured
+
+   subroutine add_cell_scalar_unstructured(this, name, data)
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      class(vtk_unstructured_grid), intent(inout) :: this
+      call add_array_generic(this%cell_scalars, name, data, 1)
+   end subroutine add_cell_scalar_unstructured
+
+   subroutine add_cell_vector_unstructured(this, name, data)
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      class(vtk_unstructured_grid), intent(inout) :: this
+      call add_array_generic(this%cell_vectors, name, data, 3)
+   end subroutine add_cell_vector_unstructured
+
+   subroutine add_cell_connectivity(this, conn, offsets, types)
+      class(vtk_unstructured_grid), intent(inout) :: this
+      integer, intent(in) :: conn(:), offsets(:), types(:)
+      this%num_cells = size(offsets)
+      if (allocated(this%connectivity)) deallocate (this%connectivity)
+      if (allocated(this%offsets)) deallocate (this%offsets)
+      if (allocated(this%types)) deallocate (this%types)
+      allocate (this%connectivity(size(conn)))
+      allocate (this%offsets(size(offsets)))
+      allocate (this%types(size(types)))
+      this%connectivity = conn
+      this%offsets = offsets
+      this%types = types
+   end subroutine add_cell_connectivity
+
+   !==========================
+   ! Write VTU
+   !==========================
+   subroutine write_vtu_file(this, filename)
+      character(len=*), intent(in) :: filename
+      class(vtk_unstructured_grid), intent(in) :: this
+      integer :: iunit
+      open (newunit=iunit, file=filename, status='replace', action='write', form='formatted')
+      write (iunit, *) '<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">'
+      write (iunit, *) '  <UnstructuredGrid>'
+      write (iunit, *) '    <Piece NumberOfPoints="', this%num_points, '" NumberOfCells="', this%num_cells, '">'
+      call write_pointdata(iunit, this%scalars, this%vectors)
+      call write_celldata(iunit, this%cell_scalars, this%cell_vectors)
+      call write_cells(iunit, this%connectivity, this%offsets, this%types)
+      call write_points(iunit, this%points)
+      write (iunit, *) '    </Piece>'
+      write (iunit, *) '  </UnstructuredGrid>'
+      write (iunit, *) '</VTKFile>'
+      close (iunit)
+   end subroutine write_vtu_file
+
+!==========================
+!==== Shared helper routines ====
+!==========================
+   subroutine add_array_generic(array_list, name, data, ncomp)
+      type(vtk_data_array), allocatable, intent(inout) :: array_list(:)
+      character(len=*), intent(in) :: name
+      real, intent(in) :: data(:)
+      integer, intent(in) :: ncomp
+      type(vtk_data_array), allocatable :: tmp(:)
+      integer :: n
+      if (.not. allocated(array_list)) then
+         allocate (array_list(1))
+         n = 1
+      else
+         n = size(array_list) + 1
+         allocate (tmp(n))
+         tmp(1:n - 1) = array_list
+         call move_alloc(tmp, array_list)
+      end if
+      array_list(n)%name = name
+      array_list(n)%type = 'Float32'
+      array_list(n)%num_components = ncomp
+      array_list(n)%data = data
+   end subroutine add_array_generic
+
+   subroutine write_pointdata(iunit, scalars, vectors)
+      integer, intent(in) :: iunit
+      type(vtk_data_array), allocatable, intent(in) :: scalars(:)
+      type(vtk_data_array), allocatable, intent(in) :: vectors(:)
+      integer :: i
+      character(len=1024) :: tag
+      tag = '      <PointData'
+      if (allocated(scalars)) then
+         if (size(scalars) > 0) tag = trim(tag)//' Scalars="'//trim(scalars(1)%name)//'"'
+      end if
+      if (allocated(vectors)) then
+         if (size(vectors) > 0) tag = trim(tag)//' Vectors="'//trim(vectors(1)%name)//'"'
+      end if
+      tag = trim(tag)//'>'
+      write (iunit, '(A)') trim(tag)
+      if (allocated(scalars)) then
+         do i = 1, size(scalars)
+            write (iunit, '(A)') '        <DataArray type="Float32" Name="'//trim(scalars(i)%name)//'" format="ascii">'
+            write (iunit, '(1000(F12.6,1X))') scalars(i)%data
+            write (iunit, '(A)') '        </DataArray>'
+         end do
+      end if
+      if (allocated(vectors)) then
+         do i = 1, size(vectors)
+  write (iunit, '(A)') '        <DataArray type="Float32" Name="'//trim(vectors(i)%name)//'" NumberOfComponents="3" format="ascii">'
+            write (iunit, '(1000(F12.6,1X))') vectors(i)%data
+            write (iunit, '(A)') '        </DataArray>'
+         end do
+      end if
+      write (iunit, *) '      </PointData>'
+   end subroutine write_pointdata
+
+   subroutine write_celldata(iunit, scalars, vectors)
+      integer, intent(in) :: iunit
+      type(vtk_data_array), allocatable, intent(in) :: scalars(:)
+      type(vtk_data_array), allocatable, intent(in) :: vectors(:)
+      integer :: i
+      if (.not. allocated(scalars) .and. .not. allocated(vectors)) then
+         write (iunit, *) '      <CellData>'
+         write (iunit, *) '      </CellData>'
+         return
+      end if
+      write (iunit, *) '      <CellData>'
+      if (allocated(scalars)) then
+         do i = 1, size(scalars)
+            write (iunit, '(A)') '        <DataArray type="Float32" Name="'//trim(scalars(i)%name)//'" format="ascii">'
+            write (iunit, '(1000(F12.6,1X))') scalars(i)%data
+            write (iunit, '(A)') '        </DataArray>'
+         end do
+      end if
+      if (allocated(vectors)) then
+         do i = 1, size(vectors)
+  write (iunit, '(A)') '        <DataArray type="Float32" Name="'//trim(vectors(i)%name)//'" NumberOfComponents="3" format="ascii">'
+            write (iunit, '(1000(F12.6,1X))') vectors(i)%data
+            write (iunit, '(A)') '        </DataArray>'
+         end do
+      end if
+      write (iunit, *) '      </CellData>'
+   end subroutine write_celldata
+
+   subroutine write_points(iunit, pts)
+      integer, intent(in) :: iunit
+      real, intent(in) :: pts(:, :)
+      write (iunit, *) '      <Points>'
+      write (iunit, *) '        <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii">'
+      write (iunit, '(1000(F12.6,1X))') pts
+      write (iunit, *) '        </DataArray>'
+      write (iunit, *) '      </Points>'
+   end subroutine write_points
+
+   subroutine write_cells(iunit, conn, offsets, types)
+      integer, intent(in) :: iunit
+      integer, intent(in) :: conn(:), offsets(:), types(:)
+      write (iunit, *) '      <Cells>'
+      write (iunit, *) '        <DataArray type="Int32" Name="connectivity" format="ascii">'
+      write (iunit, '(1000(I8,1X))') conn
+      write (iunit, *) '        </DataArray>'
+      write (iunit, *) '        <DataArray type="Int32" Name="offsets" format="ascii">'
+      write (iunit, '(1000(I8,1X))') offsets
+      write (iunit, *) '        </DataArray>'
+      write (iunit, *) '        <DataArray type="UInt8" Name="types" format="ascii">'
+      write (iunit, '(1000(I3,1X))') types
+      write (iunit, *) '        </DataArray>'
+      write (iunit, *) '      </Cells>'
+   end subroutine write_cells
+
+end module vtkAPI_m