diff --git a/doc/pages/user-guide/case-file.md b/doc/pages/user-guide/case-file.md
index 995dba2c409..944fb8d3237 100644
--- a/doc/pages/user-guide/case-file.md
+++ b/doc/pages/user-guide/case-file.md
@@ -215,6 +215,36 @@ file documentation.
    `base_value` keyword, and then assigned a zone value inside a point zone. The
    point zone is specified by the `name` keyword, and should be defined in the
    `case.point_zones` object. See more about point zones @ref point-zones.md.
+5. `field`, where the initial condition is retrieved from a field file.
+   The following keywords can be used:
+   
+| Name             | Description                                                                                        | Admissible values            | Default value |
+| ---------------- | -------------------------------------------------------------------------------------------------- | ---------------------------- | ------------- |
+| `file_name`      | Name of the field file to use (e.g. `myfield0.f00034`).                                            | Strings ending with `f*****` | -             |
+| `interpolate`    | Whether to interpolate the velocity and pressure fields from the field file onto the current mesh. | `true` or `false`            | `false`       |
+| `tolerance`      | Tolerance for the point search.                                                                    | Positive real.               | `1e-6`        |
+| `mesh_file_name` | If interpolation is enabled, the name of the field file that contains the mesh coordinates.        | Strings ending with `f*****` | `file_name`   |
+
+   @attention Interpolating a field from the same mesh but different 
+   polynomial order is performed implicitly and does not require to enable
+   interpolation.
+   
+   @note It is recommended to interpolate from `fld` files that were 
+   written in double precision. 
+   To check if your `fld` file was written in double precision, run
+   the command:
+   ~~~~~~~~~~~~~~~{.sh}
+   head -1 field0.f00000
+   ~~~~~~~~~~~~~~~
+   The output `#std 4 ...` indicates single precision, 
+   whereas `#std 8 ...` indicates double precision.
+   Neko write single precision `fld` files by default. To write your 
+   files in double precision, set `case.output_precision` to
+   `"double"`. 
+  
+   @attention Neko does not detect wether interpolation is needed or not. 
+   Interpolation will always be performed if `"interpolate"` is set 
+   to `true` even if the field file matches with the current simulation. 
 
 ### Blasius profile
 The `blasius` object is used to specify the Blasius profile that can be used for the
@@ -467,8 +497,12 @@ that can be described concisely directly in the table.
 | `output_value`                          | The frequency of sampling in terms of `output_control`.                                           | Positive real or integer                         | -             |
 | `inflow_condition.type`                 | Velocity inflow condition type.                                                                   | `user`, `uniform`, `blasius`                     | -             |
 | `inflow_condition.value`                | Value of the inflow velocity.                                                                     | Vector of 3 reals                                | -             |
-| `initial_condition.type`                | Initial condition type.                                                                           | `user`, `uniform`, `blasius`                     | -             |
+| `initial_condition.type`                | Initial condition type.                                                                           | `user`, `uniform`, `blasius`, `field`            | -             |
 | `initial_condition.value`               | Value of the velocity initial condition.                                                          | Vector of 3 reals                                | -             |
+| `initial_condition.file_name`           | If `"type" = "field"`, the path to the field file to read from.                                   | String ending with `.fld`, `.chkp`, `.nek5000` or `f*****`.  | -             |
+| `initial_condition.sample_index`        | If `"type" = "field"`, and file type is `fld` or `nek5000`, the index of the file to sampled.     | Positive integer.                                | -1            |
+| `initial_condition.previous_mesh`       | If `"type" = "field"`, and file type is `chkp`, the previous mesh from which to interpolate.      | String ending with `.nmsh`.                      | -             |
+| `initial_condition.tolerance`           | If `"type" = "field"`, and file type is `chkp`, tolerance to use for mesh interpolation.          | Positive real.                                   | 1e-6          |
 | `blasius.delta`                         | Boundary layer thickness in the Blasius profile.                                                  | Positive real                                    | -             |
 | `blasius.freestream_velocity`           | Free-stream velocity in the Blasius profile.                                                      | Vector of 3 reals                                | -             |
 | `blasius.approximation`                 | Numerical approximation of the Blasius profile.                                                   | `linear`, `quadratic`, `cubic`, `quartic`, `sin` | -             |
@@ -544,7 +578,10 @@ file documentation.
    `base_value` keyword, and then assigned a zone value inside a point zone. The
    point zone is specified by the `name` keyword, and should be defined in the
    `case.point_zones` object. See more about point zones @ref point-zones.md.
-
+4. `field`, where the initial condition is retrieved from a field file. Works
+   in the same way as for the fluid. See the 
+   [fluid section](@ref case-file_fluid-ic) for detailed explanations.
+  
 ### Source terms
 
 The configuration of source terms is the same as for the fluid. A demonstration
diff --git a/src/.depends b/src/.depends
index d7c005b40fa..2ea758815a5 100644
--- a/src/.depends
+++ b/src/.depends
@@ -113,7 +113,7 @@ io/map_file.o : io/map_file.f90 io/format/map.o comm/comm.o common/utils.o io/ge
 io/re2_file.o : io/re2_file.f90 common/log.o adt/htable.o io/map_file.o io/format/map.o io/format/re2.o common/datadist.o comm/mpi_types.o comm/comm.o mesh/point.o mesh/mesh.o common/utils.o config/num_types.o io/generic_file.o 
 io/rea_file.o : io/rea_file.f90 common/log.o adt/htable.o common/datadist.o comm/comm.o io/map_file.o io/re2_file.o io/format/rea.o io/format/map.o mesh/point.o mesh/mesh.o common/utils.o config/num_types.o io/generic_file.o 
 io/fld_file.o : io/fld_file.f90 comm/mpi_types.o math/math.o common/datadist.o comm/comm.o common/utils.o mesh/mesh.o fluid/mean_sqr_flow.o fluid/mean_flow.o io/fld_file_data.o math/vector.o common/structs.o sem/space.o sem/dofmap.o field/field_list.o field/field.o io/generic_file.o config/num_types.o 
-io/fld_file_data.o : io/fld_file_data.f90 math/vector.o math/math.o config/num_types.o 
+io/fld_file_data.o : io/fld_file_data.f90 mesh/mesh.o common/utils.o common/global_interpolation.o sem/space.o sem/dofmap.o field/field.o math/vector.o math/math.o config/num_types.o 
 io/vtk_file.o : io/vtk_file.f90 comm/comm.o common/log.o mesh/tri_mesh.o mesh/tet_mesh.o field/mesh_field.o sem/dofmap.o field/field.o mesh/mesh.o common/utils.o io/generic_file.o config/num_types.o 
 io/stl_file.o : io/stl_file.f90 io/format/stl.o comm/comm.o common/utils.o comm/mpi_types.o mesh/point.o common/log.o mesh/tri_mesh.o io/generic_file.o config/num_types.o 
 io/nmsh_file.o : io/nmsh_file.f90 common/log.o comm/mpi_types.o common/datadist.o mesh/element.o io/format/nmsh.o adt/tuple.o mesh/point.o common/utils.o mesh/mesh.o config/num_types.o comm/comm.o io/generic_file.o 
@@ -130,7 +130,7 @@ io/fluid_stats_output.o : io/fluid_stats_output.f90 math/matrix.o io/output.o de
 io/mean_sqr_flow_output.o : io/mean_sqr_flow_output.f90 io/output.o config/num_types.o fluid/mean_sqr_flow.o 
 io/data_streamer.o : io/data_streamer.F90 config/neko_config.o comm/mpi_types.o comm/comm.o device/device.o common/utils.o sem/coef.o field/field.o config/num_types.o 
 common/sampler.o : common/sampler.f90 common/time_based_controller.o config/num_types.o common/profiler.o common/utils.o common/log.o comm/comm.o io/fld_file.o io/output.o 
-common/global_interpolation.o : common/global_interpolation.F90 comm/mpi_types.o math/math.o comm/comm.o sem/local_interpolation.o common/utils.o common/log.o mesh/mesh.o sem/dofmap.o sem/space.o config/num_types.o 
+common/global_interpolation.o : common/global_interpolation.F90 common/structs.o comm/mpi_types.o math/math.o comm/comm.o sem/local_interpolation.o common/utils.o common/log.o mesh/mesh.o sem/dofmap.o sem/space.o config/num_types.o 
 common/profiler.o : common/profiler.F90 common/runtime_statistics.o common/craypat.o device/hip/roctx.o device/cuda/nvtx.o device/device.o config/neko_config.o 
 common/craypat.o : common/craypat.F90 
 bc/bc.o : bc/bc.f90 common/utils.o adt/tuple.o adt/stack.o mesh/facet_zone.o mesh/mesh.o sem/space.o sem/coef.o sem/dofmap.o device/device.o config/num_types.o config/neko_config.o 
@@ -208,7 +208,7 @@ fluid/mean_flow.o : fluid/mean_flow.f90 field/field.o field/mean_field.o
 fluid/fluid_stats.o : fluid/fluid_stats.f90 common/utils.o config/neko_config.o device/device.o common/stats_quant.o field/field_list.o field/field.o sem/coef.o math/operators.o math/math.o config/num_types.o math/bcknd/device/device_math.o field/mean_field.o 
 fluid/mean_sqr_flow.o : fluid/mean_sqr_flow.f90 field/field.o field/mean_sqr_field.o 
 fluid/flow_profile.o : fluid/flow_profile.f90 config/num_types.o 
-fluid/flow_ic.o : fluid/flow_ic.f90 mesh/point_zone_registry.o mesh/point_zone.o common/json_utils.o common/user_intf.o math/bcknd/device/device_math.o math/math.o sem/coef.o common/utils.o field/field.o device/device.o fluid/flow_profile.o config/neko_config.o gs/gather_scatter.o config/num_types.o 
+fluid/flow_ic.o : fluid/flow_ic.f90 sem/space.o sem/interpolation.o common/global_interpolation.o io/file.o io/fld_file.o io/fld_file_data.o mesh/point_zone_registry.o mesh/point_zone.o common/json_utils.o common/user_intf.o math/bcknd/device/device_math.o math/math.o sem/coef.o common/utils.o field/field.o device/device.o fluid/flow_profile.o config/neko_config.o gs/gather_scatter.o common/log.o config/num_types.o 
 fluid/advection.o : fluid/advection.f90 time_schemes/time_scheme_controller.o field/field_series.o sem/coef.o field/field.o sem/space.o config/num_types.o 
 fluid/advection_fctry.o : fluid/advection_fctry.f90 fluid/bcknd/advection/adv_oifs.o fluid/bcknd/advection/adv_no_dealias.o fluid/bcknd/advection/adv_dealias.o common/json_utils.o fluid/advection.o 
 fluid/bcknd/advection/adv_dealias.o : fluid/bcknd/advection/adv_dealias.f90 device/device.o sem/interpolation.o math/operators.o config/neko_config.o math/bcknd/device/device_math.o sem/coef.o field/field.o sem/space.o math/math.o config/num_types.o fluid/advection.o 
@@ -268,7 +268,7 @@ scalar/scalar_pnpn.o : scalar/scalar_pnpn.f90 field/scratch_registry.o common/ti
 scalar/scalar_aux.o : scalar/scalar_aux.f90 krylov/krylov.o config/num_types.o common/log.o 
 scalar/scalar_residual.o : scalar/scalar_residual.f90 config/num_types.o mesh/mesh.o sem/space.o bc/facet_normal.o scalar/source_scalar.o sem/coef.o field/field.o math/ax.o gs/gather_scatter.o 
 scalar/scalar_residual_fctry.o : scalar/scalar_residual_fctry.f90 scalar/bcknd/sx/scalar_residual_sx.o scalar/bcknd/cpu/scalar_residual_cpu.o scalar/bcknd/device/scalar_residual_device.o config/neko_config.o scalar/scalar_residual.o 
-scalar/scalar_ic.o : scalar/scalar_ic.f90 mesh/point_zone_registry.o mesh/point_zone.o common/json_utils.o common/user_intf.o math/math.o sem/coef.o common/utils.o field/field.o device/device.o math/bcknd/device/device_math.o config/num_types.o config/neko_config.o gs/gather_scatter.o 
+scalar/scalar_ic.o : scalar/scalar_ic.f90 sem/space.o sem/interpolation.o common/global_interpolation.o io/file.o common/checkpoint.o io/fld_file.o io/fld_file_data.o common/log.o field/field_registry.o mesh/point_zone_registry.o mesh/point_zone.o common/json_utils.o common/user_intf.o math/math.o sem/coef.o common/utils.o field/field.o device/device.o math/bcknd/device/device_math.o config/num_types.o config/neko_config.o gs/gather_scatter.o 
 scalar/scalar_source_term.o : scalar/scalar_source_term.f90 common/user_intf.o sem/coef.o field/field_list.o field/field.o source_terms/source_term_handler.o source_terms/source_term.o scalar/scalar_user_source_term.o 
 scalar/scalar_user_source_term.o : scalar/scalar_user_source_term.f90 sem/dofmap.o math/math.o math/bcknd/device/device_math.o device/device.o sem/coef.o field/field_list.o source_terms/source_term.o common/utils.o config/num_types.o config/neko_config.o 
 scalar/bcknd/cpu/scalar_residual_cpu.o : scalar/bcknd/cpu/scalar_residual_cpu.f90 math/math.o config/num_types.o mesh/mesh.o sem/space.o sem/coef.o field/field.o math/ax.o scalar/scalar_residual.o 
diff --git a/src/case.f90 b/src/case.f90
index 3342bec7670..dd0736c2d6e 100644
--- a/src/case.f90
+++ b/src/case.f90
@@ -103,7 +103,7 @@ subroutine case_init_from_file(C, case_file)
     end if
 
     call MPI_Bcast(integer_val, 1, MPI_INTEGER, 0, NEKO_COMM, ierr)
-    if (pe_rank .ne. 0) allocate(character(len=integer_val)::json_buffer)
+    if (pe_rank .ne. 0) allocate(character(len = integer_val) :: json_buffer)
     call MPI_Bcast(json_buffer, integer_val, MPI_CHARACTER, 0, NEKO_COMM, ierr)
     call C%params%load_from_string(json_buffer)
 
@@ -137,7 +137,7 @@ subroutine case_init_common(C)
     logical :: found, logical_val
     integer :: integer_val
     real(kind=rp) :: real_val
-    character(len=:), allocatable :: string_val
+    character(len = :), allocatable :: string_val
     real(kind=rp) :: stats_start_time, stats_output_val
     integer :: stats_sampling_interval
     integer :: output_dir_len
@@ -270,6 +270,9 @@ subroutine case_init_common(C)
     !
     call json_get(C%params, 'case.fluid.initial_condition.type',&
                   string_val)
+
+    call neko_log%section("Fluid initial condition ")
+
     if (trim(string_val) .ne. 'user') then
        call set_flow_ic(C%fluid%u, C%fluid%v, C%fluid%w, C%fluid%p, &
             C%fluid%c_Xh, C%fluid%gs_Xh, string_val, C%params)
@@ -278,8 +281,14 @@ subroutine case_init_common(C)
             C%fluid%c_Xh, C%fluid%gs_Xh, C%usr%fluid_user_ic, C%params)
     end if
 
+    call neko_log%end_section()
+
     if (scalar) then
+
        call json_get(C%params, 'case.scalar.initial_condition.type', string_val)
+
+       call neko_log%section("Scalar initial condition ")
+
        if (trim(string_val) .ne. 'user') then
           call set_scalar_ic(C%scalar%s, &
             C%scalar%c_Xh, C%scalar%gs_Xh, string_val, C%params)
@@ -287,6 +296,9 @@ subroutine case_init_common(C)
           call set_scalar_ic(C%scalar%s, &
             C%scalar%c_Xh, C%scalar%gs_Xh, C%usr%scalar_user_ic, C%params)
        end if
+
+       call neko_log%end_section()
+
     end if
 
     ! Add initial conditions to BDF scheme (if present)
diff --git a/src/common/global_interpolation.F90 b/src/common/global_interpolation.F90
index 5d9e703182b..60d3710a102 100644
--- a/src/common/global_interpolation.F90
+++ b/src/common/global_interpolation.F90
@@ -46,49 +46,57 @@ module global_interpolation
   use comm
   use math, only: copy
   use neko_mpi_types
+  use structs, only: array_ptr_t
   use, intrinsic :: iso_c_binding
   implicit none
   private
   !> Implements global interpolation for arbitrary points in the domain.
   type, public :: global_interpolation_t
-     !> Dofmap from which we interpolate the points
-     type(dofmap_t), pointer :: dof
-     !> Mesh on which we interpolate
-     type(mesh_t), pointer :: mesh
-     !> Space
+     !> X coordinates from which to interpolate.
+     type(array_ptr_t) :: x
+     !> Y coordinates from which to interpolate.
+     type(array_ptr_t) :: y
+     !> Z coordinates from which to interpolate.
+     type(array_ptr_t) :: z
+     !> Geometric dimension of the simulation.
+     integer :: gdim
+     !> Number of elements.
+     integer :: nelv
+     !> Space.
      type(space_t), pointer :: Xh
-     !> Interpolator for local points
+     !> Interpolator for local points.
      type(local_interpolator_t) :: local_interp
-     !> If all points are local on this PE
+     !> If all points are local on this PE.
      logical :: all_points_local = .false.
-     !! Gslib handle
-     !! @note: Remove when we remove gslib
+     !! Gslib handle.
+     !! @note: Remove when we remove gslib.
      integer :: gs_handle
      logical :: gs_init = .false.
      !> Components related to the points we want to evalute
      !> Number of points we want to evaluate
      integer :: n_points
-     !> x,y,z coordinates, findpts format
+     !> x,y,z coordinates, findpts format.
      !! @note: When replacing gs we can change format
      real(kind=rp), allocatable :: xyz(:,:)
-     !> List of owning processes
+     !> List of owning processes.
      integer, allocatable :: proc_owner(:)
-     !> List of owning elements
+     !> List of owning elements.
      integer, allocatable :: el_owner(:)
      type(c_ptr) :: el_owner_d = c_null_ptr
-     !> r,s,t coordinates findpts format
+     !> r,s,t coordinates findpts format.
      !! @note: When replacing gs we can change format
      real(kind=rp), allocatable :: rst(:,:)
-     !> Distance squared between original and interpolated point
+     !> Distance squared between original and interpolated point.
      !! (in xyz space) (according to gslib)
      real(kind=rp), allocatable :: dist2(:)
-     !> Error code for each point, needed for gslib
+     !> Error code for each point, needed for gslib.
      integer, allocatable :: error_code(:)
-     !> Tolerance for distance squared between original and interpolated point
-     real(kind=rp) :: tol = 5e-13
+     !> Tolerance for distance squared between original and interpolated point.
+     real(kind=rp) :: tol = 5d-13
    contains
      !> Initialize the global interpolation object on a dofmap.
-     procedure, pass(this) :: init => global_interpolation_init
+     procedure, pass(this) :: init_xyz => global_interpolation_init_xyz
+     procedure, pass(this) :: init_dof => global_interpolation_init_dof
      !> Destructor
      procedure, pass(this) :: free => global_interpolation_free
      !> Destructor for arrays related to evaluation points
@@ -98,43 +106,76 @@ module global_interpolation
      !> Finds the process owner, global element number,
      !! and local rst coordinates for each point.
      !! Sets up correct values to be able to evalute the points
-     procedure, pass(this) :: find_points_coords => global_interpolation_find_coords
+     procedure, pass(this) :: find_points_coords => &
+          global_interpolation_find_coords
      procedure, pass(this) :: find_points_xyz => global_interpolation_find_xyz
      generic :: find_points => find_points_xyz, find_points_coords
      !> Evaluate the value of the field in each point.
      procedure, pass(this) :: evaluate => global_interpolation_evaluate
-     !> Evaluate only local points
-  end type global_interpolation_t
 
+     !> Generic constructor
+     generic :: init => init_dof, init_xyz
+
+  end type global_interpolation_t
 
 contains
-  !> Initialize user defined variables.
+
+  !> Initialize the global interpolation object on a dofmap.
   !! @param dof Dofmap on which the interpolation is to be carried out.
   !! @param tol Tolerance for Newton iterations.
-  subroutine global_interpolation_init(this, dof, tol)
+  subroutine global_interpolation_init_dof(this, dof, tol)
     class(global_interpolation_t), intent(inout) :: this
     type(dofmap_t), target :: dof
     real(kind=rp), optional :: tol
-    integer :: lx, ly, lz, nelv, max_pts_per_iter
 
-    this%dof => dof
-    this%Xh => dof%Xh
-    this%mesh => dof%msh
-    if(present(tol)) this%tol = tol
+    ! NOTE: Passing dof%x(:,1,1,1), etc in init_xyz passes down the entire
+    ! dof%x array and not a slice. It is done this way for
+    ! this%x%ptr to point to dof%x (see global_interpolation_init_xyz).
+    call this%init_xyz(dof%x(:,1,1,1), dof%y(:,1,1,1), dof%z(:,1,1,1), &
+         dof%msh%gdim, dof%msh%nelv, dof%Xh, tol = tol)
+
+  end subroutine global_interpolation_init_dof
+
+  !> Initialize the global interpolation object on a set of coordinates.
+  !! @param x x-coordinates.
+  !! @param y y-coordinates.
+  !! @param z z-coordinates.
+  !! @param gdim Geometric dimension.
+  !! @param nelv Number of elements of the mesh in which to search for the
+  !! points.
+  !! @param Xh Space on which to interpolate.
+  !! @param tol Tolerance for Newton iterations.
+  subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, Xh, tol)
+    class(global_interpolation_t), intent(inout) :: this
+    real(kind=rp), intent(in), target :: x(:)
+    real(kind=rp), intent(in), target :: y(:)
+    real(kind=rp), intent(in), target :: z(:)
+    integer, intent(in) :: gdim
+    integer, intent(in) :: nelv
+    type(space_t), intent(in), target :: Xh
+    real(kind=rp), intent(in), optional :: tol
+    integer :: lx, ly, lz, max_pts_per_iter
 
 #ifdef HAVE_GSLIB
 
-    lx = this%Xh%lx
-    ly = this%Xh%ly
-    lz = this%Xh%lz
-    nelv = this%mesh%nelv
-    !Number of points to iterate on simultaneosuly
+    this%x%ptr => x
+    this%y%ptr => y
+    this%z%ptr => z
+    this%gdim = gdim
+    this%nelv = nelv
+    this%Xh => Xh
+    if (present(tol)) this%tol = tol
+
+    ! Number of points to iterate on simultaneosuly
     max_pts_per_iter = 128
+    lx = Xh%lx
+    ly = Xh%ly
+    lz = Xh%lz
 
     call fgslib_findpts_setup(this%gs_handle, &
          NEKO_COMM, pe_size, &
-         this%mesh%gdim, &
-         dof%x, dof%y, dof%z, & ! Physical nodal values
+         this%gdim, &
+         this%x%ptr, this%y%ptr, this%z%ptr, & ! Physical nodal values
          lx, ly, lz, nelv, & ! Mesh dimensions
          2*lx, 2*ly, 2*lz, & ! Mesh size for bounding box computation
          0.01, & ! relative size to expand bounding boxes by
@@ -145,17 +186,21 @@ subroutine global_interpolation_init(this, dof, tol)
     call neko_error('Neko needs to be built with GSLIB support')
 #endif
 
-  end subroutine global_interpolation_init
+  end subroutine global_interpolation_init_xyz
 
 
   !> Destructor
   subroutine global_interpolation_free(this)
     class(global_interpolation_t), intent(inout) :: this
 
-    nullify(this%mesh)
-    nullify(this%dof)
+    nullify(this%x%ptr)
+    nullify(this%y%ptr)
+    nullify(this%z%ptr)
     nullify(this%Xh)
 
+    this%nelv = 0
+    this%gdim = 0
+
     call this%free_points()
 
 #ifdef HAVE_GSLIB
@@ -174,11 +219,11 @@ subroutine global_interpolation_free_points(this)
     this%n_points = 0
     this%all_points_local = .false.
 
-    if (allocated(this%xyz))        deallocate(this%xyz)
-    if (allocated(this%rst))        deallocate(this%rst)
+    if (allocated(this%xyz)) deallocate(this%xyz)
+    if (allocated(this%rst)) deallocate(this%rst)
     if (allocated(this%proc_owner)) deallocate(this%proc_owner)
-    if (allocated(this%el_owner))   deallocate(this%el_owner)
-    if (allocated(this%dist2))       deallocate(this%dist2)
+    if (allocated(this%el_owner)) deallocate(this%el_owner)
+    if (allocated(this%dist2)) deallocate(this%dist2)
     if (allocated(this%error_code)) deallocate(this%error_code)
 
     if (c_associated(this%el_owner_d)) then
@@ -207,13 +252,13 @@ subroutine global_interpolation_find_common(this)
          this%error_code, 1, &
          this%proc_owner, 1, &
          this%el_owner, 1, &
-         this%rst, this%mesh%gdim, &
+         this%rst, this%gdim, &
          this%dist2, 1, &
-         this%xyz(1,1), this%mesh%gdim, &
-         this%xyz(2,1), this%mesh%gdim, &
-         this%xyz(3,1), this%mesh%gdim, this%n_points)
+         this%xyz(1,1), this%gdim, &
+         this%xyz(2,1), this%gdim, &
+         this%xyz(3,1), this%gdim, this%n_points)
 
-    do i=1,this%n_points
+    do i = 1 , this%n_points
 
        !
        ! Check validity of points
@@ -236,7 +281,8 @@ subroutine global_interpolation_find_common(this)
                 ' Interpolation on these points will return 0.0. dist2: ', &
                 this%dist2(i),&
                 'el_owner, rst coords, pe: ',&
-                this%el_owner(i), this%rst(1,i), this%rst(2,i), this%rst(3,i), pe_rank
+                this%el_owner(i), this%rst(1,i), this%rst(2,i), &
+                this%rst(3,i), pe_rank
 
     end do
 
@@ -247,23 +293,22 @@ subroutine global_interpolation_find_common(this)
     call fgslib_findpts_eval(this%gs_handle, x_check, &
                              1, this%error_code, 1, &
                              this%proc_owner, 1, this%el_owner, 1, &
-                             this%rst, this%mesh%gdim, &
-                             this%n_points, this%dof%x)
+                             this%rst, this%gdim, &
+                             this%n_points, this%x%ptr)
 
     call fgslib_findpts_eval(this%gs_handle, y_check, &
                              1, this%error_code, 1, &
                              this%proc_owner, 1, this%el_owner, 1, &
-                             this%rst, this%mesh%gdim, &
-                             this%n_points, this%dof%y)
+                             this%rst, this%gdim, &
+                             this%n_points, this%y%ptr)
 
     call fgslib_findpts_eval(this%gs_handle, z_check, &
                              1, this%error_code, 1, &
                              this%proc_owner, 1, this%el_owner, 1, &
-                             this%rst, this%mesh%gdim, &
-                             this%n_points, this%dof%z)
-
+                             this%rst, this%gdim, &
+                             this%n_points, this%z%ptr)
 
-    do i=1,this%n_points
+    do i = 1 , this%n_points
 
        !
        ! Check validity of points
@@ -277,11 +322,11 @@ subroutine global_interpolation_find_common(this)
        if ( zdiff .gt. this%tol) isdiff = .true.
 
        if (isdiff) then
-          write(*,*) 'Points with coords: ',&
-                this%xyz(1,i),this%xyz(2,i),this%xyz(3,i), &
-                'Differ from interpolated coords: ',&
-                x_check(i), y_check(i), z_check(i),&
-                'Distance squared: ',&
+          write(*,*) 'Points with coords: ', &
+                this%xyz(1, i), this%xyz(2, i), this%xyz(3, i), &
+                'Differ from interpolated coords: ', &
+                x_check(i), y_check(i), z_check(i), &
+                'Distance squared: ', &
                 xdiff, ydiff, zdiff
        end if
 
@@ -304,8 +349,8 @@ end subroutine global_interpolation_find_common
   !! in the correct global element as well as which process that owns the point.
   !! After this the values at these points can be evaluated.
   !! If the locations of the points change this must be called again.
-  !! - `error_code`: returns `0` if point found, `1` if closest point on a border
-  !! (check dist2), `2` if not found
+  !! - `error_code`: returns `0` if point found, `1` if closest point on a
+  !! border (check dist2), `2` if not found
   !! - `dist2`: distance squared (used to compare the points found by each
   !! processor)
   !! @param x The x-coordinates of the points.
@@ -330,9 +375,9 @@ subroutine global_interpolation_find_coords(this, x, y, z, n_points)
     call global_interpolation_init_point_arrays(this)
 
     do i = 1, n_points
-       this%xyz(1,i) = x(i,1,1,1)
-       this%xyz(2,i) = y(i,1,1,1)
-       this%xyz(3,i) = z(i,1,1,1)
+       this%xyz(1, i) = x(i,1,1,1)
+       this%xyz(2, i) = y(i,1,1,1)
+       this%xyz(3, i) = z(i,1,1,1)
     end do
 
     call global_interpolation_find_common(this)
@@ -342,15 +387,15 @@ end subroutine global_interpolation_find_coords
   subroutine global_interpolation_init_point_arrays(this)
     class(global_interpolation_t) :: this
 
-    allocate(this%xyz(3,this%n_points))
-    allocate(this%rst(3,this%n_points))
+    allocate(this%xyz(3, this%n_points))
+    allocate(this%rst(3, this%n_points))
     allocate(this%proc_owner(this%n_points))
     allocate(this%el_owner(this%n_points))
     allocate(this%dist2(this%n_points))
     allocate(this%error_code(this%n_points))
 
     if (NEKO_BCKND_DEVICE .eq. 1) &
-       call device_map(this%el_owner, this%el_owner_d,this%n_points)
+       call device_map(this%el_owner, this%el_owner_d, this%n_points)
 
   end subroutine global_interpolation_init_point_arrays
 
@@ -358,8 +403,8 @@ end subroutine global_interpolation_init_point_arrays
   !! in the correct global element as well as which process that owns the point.
   !! After this the values at these points can be evaluated.
   !! If the locations of the points change this must be called again.
-  !! - `error_code`: returns `0` if point found, `1` if closest point on a border
-  !! (check dist2), `2` if not found
+  !! - `error_code`: returns `0` if point found, `1` if closest point on a
+  !! border (check dist2), `2` if not found
   !! - `dist2`: distance squared (used to compare the points found by each
   !! processor)
   !! @param xyz The coordinates of the points.
@@ -368,7 +413,7 @@ subroutine global_interpolation_find_xyz(this, xyz, n_points)
     class(global_interpolation_t), intent(inout) :: this
     integer, intent(in) :: n_points
     !!Perhaps this should be kind dp
-    real(kind=rp), intent(inout) :: xyz(3,n_points)
+    real(kind=rp), intent(inout) :: xyz(3, n_points)
 
 
     call this%free_points()
@@ -378,18 +423,19 @@ subroutine global_interpolation_find_xyz(this, xyz, n_points)
     call global_interpolation_init_point_arrays(this)
 
     !> make deep copy incase xyz goes out of scope or deallocated
-    call copy(this%xyz,xyz,3*n_points)
+    call copy(this%xyz, xyz, 3 * n_points)
 
     call global_interpolation_find_common(this)
 
   end subroutine global_interpolation_find_xyz
 
-  !> Finds the corresponding r,s,t coordinates and redistributes the points to the owning rank
-  !! in the correct global element as well as which process that owns the point.
+  !> Finds the corresponding r,s,t coordinates and redistributes the points to
+  !! the owning rank in the correct global element as well as which process
+  !! that owns the point.
   !! After this the values at these points can be evaluated.
   !! If the locations of the points change this must be called again.
-  !! - `error_code`: returns `0` if point found, `1` if closest point on a border
-  !! (check dist2), `2` if not found
+  !! - `error_code`: returns `0` if point found, `1` if closest point on a
+  !! border (check dist2), `2` if not found.
   !! - `dist2`: distance squared (used to compare the points found by each
   !! processor)
   !! @param xyz The coordinates of the points.
@@ -408,7 +454,7 @@ subroutine global_interpolation_find_and_redist(this, xyz, n_points)
     call global_interpolation_init_point_arrays(this)
 
     !> make deep copy incase xyz goes out of scope or deallocated
-    call copy(this%xyz,xyz,3*n_points)
+    call copy(this%xyz, xyz, 3 * n_points)
 
     call global_interpolation_find_common(this)
     !> Sets new points and redistributes them
@@ -417,17 +463,17 @@ subroutine global_interpolation_find_and_redist(this, xyz, n_points)
 
     do i = 1, this%n_points
        if (this%proc_owner(i) .ne. pe_rank) then
-          write(*,*) 'Redistribution failed on rank: ', pe_rank,&
+          write(*,*) 'Redistribution failed on rank: ', pe_rank, &
                      'for point with coord: ', &
-                     this%xyz(1,i),this%xyz(2,i),this%xyz(3,i)
+                     this%xyz(1, i), this%xyz(2, i), this%xyz(3, i)
           exit
        end if
     end do
 
     n_points = this%n_points
     deallocate(xyz)
-    allocate(xyz(3,n_points))
-    call copy(xyz,this%xyz,3*n_points)
+    allocate(xyz(3, n_points))
+    call copy(xyz, this%xyz, 3*n_points)
 
     call this%local_interp%init(this%Xh, this%rst(1,:),&
                                 this%rst(2,:), this%rst(3,:), n_points)
@@ -452,28 +498,30 @@ subroutine global_interpolation_redist(this)
     n_points_per_pe = 0
     n_points_from_pe = 0
     !> Calculate which processes this proc has points on
-    do i = 1,this%n_points
-       n_points_per_pe(this%proc_owner(i)) = n_points_per_pe(this%proc_owner(i)) + 1
+    do i = 1, this%n_points
+       n_points_per_pe(this%proc_owner(i)) = &
+            n_points_per_pe(this%proc_owner(i)) + 1
     end do
     !> Sum number of points on all pes to compute n_new_points
     !! Store how many points to receive from each pe
-    do i = 0,(pe_size-1)
-       call MPI_Reduce(n_points_per_pe(i),n_new_points,1,MPI_INTEGER,MPI_SUM, i, NEKO_COMM, ierr)
+    do i = 0, (pe_size - 1)
+       call MPI_Reduce(n_points_per_pe(i), n_new_points, 1, MPI_INTEGER, &
+            MPI_SUM, i, NEKO_COMM, ierr)
        call MPI_Gather(n_points_per_pe(i), 1, MPI_INTEGER,&
                       n_points_from_pe, 1, MPI_INTEGER, i, NEKO_COMM, ierr)
     end do
 
     allocate(n_point_offset_from_pe(0:(pe_size-1)))
     n_point_offset_from_pe(0) = 0
-    do i = 1,(pe_size-1)
+    do i = 1, (pe_size - 1)
        n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
                                  + n_point_offset_from_pe(i-1)
     end do
 
-    allocate(new_xyz(3,n_new_points))
+    allocate(new_xyz(3, n_new_points))
     max_n_points_to_send = maxval(n_points_per_pe)
-    allocate(xyz_send_to_pe(3,max_n_points_to_send))
-    do i = 0, (pe_size-1)
+    allocate(xyz_send_to_pe(3, max_n_points_to_send))
+    do i = 0, (pe_size - 1)
        !> This could be avoided by adding all indices to a list
        k = 0
        do j = 1, this%n_points
@@ -483,14 +531,14 @@ subroutine global_interpolation_redist(this)
           end if
        end do
        if (k .ne. n_points_per_pe(i)) then
-          write(*,*) 'PE: ', pe_rank, ' has k= ', k,&
-                     'points for PE:', i,' but should have: ',&
+          write(*,*) 'PE: ', pe_rank, ' has k= ', k, &
+                     'points for PE:', i,' but should have: ', &
                      n_points_per_pe(i)
           call neko_error('Error in redistribution of points')
        end if
-       call MPI_Gatherv(xyz_send_to_pe,3*n_points_per_pe(i),&
-                        MPI_DOUBLE_PRECISION, new_xyz,3*n_points_from_pe,&
-                        3*n_point_offset_from_pe,&
+       call MPI_Gatherv(xyz_send_to_pe,3*n_points_per_pe(i), &
+                        MPI_DOUBLE_PRECISION, new_xyz,3*n_points_from_pe, &
+                        3*n_point_offset_from_pe, &
                         MPI_DOUBLE_PRECISION, i, NEKO_COMM, ierr)
 
     end do
@@ -499,7 +547,7 @@ subroutine global_interpolation_redist(this)
 
     this%n_points = n_new_points
     call global_interpolation_init_point_arrays(this)
-    call copy(this%xyz,new_xyz,3*n_new_points)
+    call copy(this%xyz, new_xyz, 3 * n_new_points)
 
     deallocate(n_point_offset_from_pe)
     deallocate(n_points_from_pe)
@@ -516,19 +564,20 @@ end subroutine global_interpolation_redist
   subroutine global_interpolation_evaluate(this, interp_values, field)
     class(global_interpolation_t), intent(inout) :: this
     real(kind=rp), intent(inout) :: interp_values(this%n_points)
-    real(kind=rp), intent(inout) :: field(this%dof%size())
+    real(kind=rp), intent(inout) :: field(this%nelv*this%Xh%lxyz)
+
 
 #ifdef HAVE_GSLIB
     if (.not. this%all_points_local) then
        call fgslib_findpts_eval(this%gs_handle, interp_values, &
                                 1, this%error_code, 1, &
                                 this%proc_owner, 1, this%el_owner, 1, &
-                                this%rst, this%mesh%gdim, &
+                                this%rst, this%gdim, &
                                 this%n_points, field)
     else
        if (this%n_points .gt. 0) &
-          call this%local_interp%evaluate(interp_values, this%el_owner,&
-                                          field, this%mesh%nelv)
+          call this%local_interp%evaluate(interp_values, this%el_owner, &
+                                          field, this%nelv)
     end if
 #else
     call neko_error('Neko needs to be built with GSLIB support')
diff --git a/src/common/utils.f90 b/src/common/utils.f90
index 6be482a3959..a2811502f4a 100644
--- a/src/common/utils.f90
+++ b/src/common/utils.f90
@@ -46,7 +46,7 @@ module utils
   public :: neko_error, neko_warning, nonlinear_index, filename_chsuffix, &
        filename_suffix, filename_suffix_pos, filename_tslash_pos, &
        linear_index, split_string, NEKO_FNAME_LEN, index_is_on_facet, &
-       concat_string_array
+       concat_string_array, extract_fld_file_index
 
 
 contains
@@ -84,6 +84,51 @@ subroutine filename_chsuffix(fname, new_fname, new_suffix)
 
   end subroutine filename_chsuffix
 
+  !> Extracts the index of a field file. For example, "myfield.f00045"
+  !! will return `45`. If the suffix of the file name is invalid, returns
+  !! a default index value.
+  !! @param fld_filename Name of the fld file, e.g. `myfield0.f00035`.
+  !! @param default_index The index to return in case the suffix of
+  !! `fld_filename` is invalid.
+  function extract_fld_file_index(fld_filename, default_index) result(index)
+    character(len=*), intent(in) :: fld_filename
+    integer, intent(in) :: default_index
+
+    character(len=80) :: suffix
+    integer :: index, fpos, i
+    logical :: valid
+
+    call filename_suffix(fld_filename, suffix)
+
+    valid = .true.
+
+    ! This value will be modified when reading the file name extension
+    ! e.g. "field0.f00035" will set sample_idx = 35
+    index = default_index
+
+    !
+    ! Try to extract the index of the field file from the suffix "fxxxxx"
+    !
+    fpos = scan(trim(suffix), 'f')
+    if (fpos .eq. 1) then
+       ! Make sure that the suffix only contains integers from 0 to 9
+       do i = 2, len(trim(suffix))
+          if (.not. (iachar(suffix(i:i)) >= iachar('0') &
+               .and. iachar(suffix(i:i)) <= iachar('9'))) then
+             valid = .false.
+          end if
+       end do
+    else
+       valid = .false.
+    end if
+
+    ! Must be exactly 6 characters long, i.e. an 'f' with 5 integers after
+    if (len(trim(suffix)) .ne. 6) valid = .false.
+
+    if (valid) read (suffix(2:), "(I5.5)") index
+
+  end function extract_fld_file_index
+
   !> Split a string based on delimiter (tokenizer)
   !! OBS: very hacky, this should really be improved, it is rather embarrasing
   !! code.
diff --git a/src/fluid/flow_ic.f90 b/src/fluid/flow_ic.f90
index 8d821ab8982..e95a4eea65d 100644
--- a/src/fluid/flow_ic.f90
+++ b/src/fluid/flow_ic.f90
@@ -33,21 +33,29 @@
 !> Initial flow condition
 module flow_ic
   use num_types, only : rp
+  use logger, only: neko_log, LOG_SIZE
   use gather_scatter, only : gs_t, GS_OP_ADD
   use neko_config, only : NEKO_BCKND_DEVICE
   use flow_profile, only : blasius_profile, blasius_linear, blasius_cubic, &
     blasius_quadratic, blasius_quartic, blasius_sin
   use device, only: device_memcpy, HOST_TO_DEVICE
   use field, only : field_t
-  use utils, only : neko_error
+  use utils, only : neko_error, filename_suffix, filename_chsuffix, &
+       neko_warning, NEKO_FNAME_LEN, extract_fld_file_index
   use coefs, only : coef_t
   use math, only : col2, cfill, cfill_mask
   use device_math, only : device_col2, device_cfill, device_cfill_mask
   use user_intf, only : useric
   use json_module, only : json_file
-  use json_utils, only: json_get
+  use json_utils, only: json_get, json_get_or_default
   use point_zone, only: point_zone_t
   use point_zone_registry, only: neko_point_zone_registry
+  use fld_file_data, only: fld_file_data_t
+  use fld_file, only: fld_file_t
+  use file, only: file_t
+  use global_interpolation, only: global_interpolation_t
+  use interpolation, only: interpolator_t
+  use space, only: space_t, GLL
   implicit none
   private
 
@@ -69,28 +77,67 @@ subroutine set_flow_ic_int(u, v, w, p, coef, gs, type, params)
     type(gs_t), intent(inout) :: gs
     character(len=*) :: type
     type(json_file), intent(inout) :: params
-    real(kind=rp) :: delta
+    real(kind=rp) :: delta, tol
     real(kind=rp), allocatable :: uinf(:)
     real(kind=rp), allocatable :: zone_value(:)
-    character(len=:), allocatable :: blasius_approximation
-    character(len=:), allocatable :: zone_name
+    character(len=:), allocatable :: read_str
+    character(len=NEKO_FNAME_LEN) :: fname, mesh_fname
+    logical :: interpolate
 
+
+    !
+    ! Uniform (Uinf, Vinf, Winf)
+    !
     if (trim(type) .eq. 'uniform') then
+
        call json_get(params, 'case.fluid.initial_condition.value', uinf)
        call set_flow_ic_uniform(u, v, w, uinf)
+
+    !
+    ! Blasius boundary layer
+    !
     else if (trim(type) .eq. 'blasius') then
+
        call json_get(params, 'case.fluid.blasius.delta', delta)
        call json_get(params, 'case.fluid.blasius.approximation', &
-                     blasius_approximation)
+                     read_str)
        call json_get(params, 'case.fluid.blasius.freestream_velocity', uinf)
-       call set_flow_ic_blasius(u, v, w, delta, uinf, blasius_approximation)
+
+       call set_flow_ic_blasius(u, v, w, delta, uinf, read_str)
+
+    !
+    ! Point zone initial condition
+    !
     else if (trim(type) .eq. 'point_zone') then
+
        call json_get(params, 'case.fluid.initial_condition.base_value', uinf)
        call json_get(params, 'case.fluid.initial_condition.zone_name', &
-                     zone_name)
+                     read_str)
        call json_get(params, 'case.fluid.initial_condition.zone_value', &
-                     zone_value)
-       call set_flow_ic_point_zone(u, v, w, uinf, zone_name, zone_value)
+            zone_value)
+
+       call set_flow_ic_point_zone(u, v, w, uinf, read_str, zone_value)
+
+    !
+    ! Field initial condition (from fld file)
+    !
+    else if (trim(type) .eq. 'field') then
+
+       call json_get(params, 'case.fluid.initial_condition.file_name', &
+            read_str)
+       fname = trim(read_str)
+       call json_get_or_default(params, &
+            'case.fluid.initial_condition.interpolate', interpolate, &
+            .false.)
+       call json_get_or_default(params, &
+            'case.fluid.initial_condition.tolerance', tol, 0.000001_rp)
+       call json_get_or_default(params, &
+            'case.fluid.initial_condition.mesh_file_name', read_str, &
+            "none")
+       mesh_fname = trim(read_str)
+
+       call set_flow_ic_fld(u, v, w, p, fname, interpolate, tol, mesh_fname)
+
     else
        call neko_error('Invalid initial condition')
     end if
@@ -110,6 +157,8 @@ subroutine set_flow_ic_usr(u, v, w, p, coef, gs, usr_ic, params)
     procedure(useric) :: usr_ic
     type(json_file), intent(inout) :: params
 
+
+    call neko_log%message("Type: user")
     call usr_ic(u, v, w, p, params)
 
     call set_flow_ic_common(u, v, w, p, coef, gs)
@@ -129,11 +178,11 @@ subroutine set_flow_ic_common(u, v, w, p, coef, gs)
 
     if (NEKO_BCKND_DEVICE .eq. 1) then
        call device_memcpy(u%x, u%x_d, n, &
-                          HOST_TO_DEVICE, sync=.false.)
+                          HOST_TO_DEVICE, sync = .false.)
        call device_memcpy(v%x, v%x_d, n, &
-                          HOST_TO_DEVICE, sync=.false.)
+                          HOST_TO_DEVICE, sync = .false.)
        call device_memcpy(w%x, w%x_d, n, &
-                          HOST_TO_DEVICE, sync=.false.)
+                          HOST_TO_DEVICE, sync = .false.)
     end if
 
     ! Ensure continuity across elements for initial conditions
@@ -159,7 +208,14 @@ subroutine set_flow_ic_uniform(u, v, w, uinf)
     type(field_t), intent(inout) :: v
     type(field_t), intent(inout) :: w
     real(kind=rp), intent(in) :: uinf(3)
-    integer :: n
+    integer :: n, i
+    character(len=LOG_SIZE) :: log_buf
+
+    call neko_log%message("Type : uniform")
+    write (log_buf, '(A, 3(ES12.6, A))') "Value: [", (uinf(i), ", ", i=1, 2), &
+         uinf(3), "]"
+    call neko_log%message(log_buf)
+
     u = uinf(1)
     v = uinf(2)
     w = uinf(3)
@@ -183,17 +239,26 @@ subroutine set_flow_ic_blasius(u, v, w, delta, uinf, type)
     character(len=*), intent(in) :: type
     procedure(blasius_profile), pointer :: bla => null()
     integer :: i
+    character(len=LOG_SIZE) :: log_buf
 
-    select case(trim(type))
-    case('linear')
+    call neko_log%message("Type         : blasius")
+    write (log_buf, '(A,ES12.6)') "delta        : ", delta
+    call neko_log%message(log_buf)
+    call neko_log%message("Approximation : " // trim(type))
+    write (log_buf, '(A,"[",2(ES12.6,","),ES12.6,"]")') "Value         : ", &
+         uinf(1), uinf(2), uinf(3)
+    call neko_log%message(log_buf)
+
+    select case (trim(type))
+    case ('linear')
        bla => blasius_linear
-    case('quadratic')
+    case ('quadratic')
        bla => blasius_quadratic
-    case('cubic')
+    case ('cubic')
        bla => blasius_cubic
-    case('quartic')
+    case ('quartic')
        bla => blasius_quartic
-    case('sin')
+    case ('sin')
        bla => blasius_sin
     case default
        call neko_error('Invalid Blasius approximation')
@@ -240,11 +305,20 @@ subroutine set_flow_ic_point_zone(u, v, w, base_value, zone_name, zone_value)
     real(kind=rp), intent(in), dimension(3) :: base_value
     character(len=*), intent(in) :: zone_name
     real(kind=rp), intent(in) :: zone_value(:)
+    character(len=LOG_SIZE) :: log_buf
 
     ! Internal variables
     class(point_zone_t), pointer :: zone
     integer :: size
 
+    call neko_log%message("Type       : point_zone")
+    write (log_buf, '(A,ES12.6)') "Base value : ", base_value
+    call neko_log%message(log_buf)
+    call neko_log%message("Zone name : " // trim(zone_name))
+    write (log_buf, '(A,"[",2(ES12.6,","),ES12.6," ]")') "Value      : ", &
+         zone_value(1), zone_value(2), zone_value(3)
+    call neko_log%message(log_buf)
+
     call set_flow_ic_uniform(u, v, w, base_value)
     size = u%dof%size()
 
@@ -256,4 +330,174 @@ subroutine set_flow_ic_point_zone(u, v, w, base_value, zone_name, zone_value)
 
   end subroutine set_flow_ic_point_zone
 
+  !> Set the initial condition of the flow based on a field.
+  !! @detail The fields are read from an `fld` file. If enabled, interpolation
+  !! is also possible. In that case, the mesh coordinates can be read from
+  !! another field in the `fld` field series.
+  !! @param u The x-component of the velocity field.
+  !! @param v The y-component of the velocity field.
+  !! @param w The z-component of the velocity field.
+  !! @param p The pressure field.
+  !! @param file_name The name of the "fld" file series.
+  !! @param sample_idx index of the field file .f000* to read, default is
+  !! -1.
+  !! @param interpolate Flag to indicate wether or not to interpolate the
+  !! values onto the current mesh.
+  !! @param tolerance If interpolation is enabled, tolerance for finding the
+  !! points in the mesh.
+  !! @param sample_mesh_idx If interpolation is enabled, index of the field
+  !! file where the mesh coordinates are located.
+  subroutine set_flow_ic_fld(u, v, w, p, file_name, &
+       interpolate, tolerance, mesh_file_name)
+    type(field_t), intent(inout) :: u
+    type(field_t), intent(inout) :: v
+    type(field_t), intent(inout) :: w
+    type(field_t), intent(inout) :: p
+    character(len=*), intent(in) :: file_name
+    logical, intent(in) :: interpolate
+    real(kind=rp), intent(in) :: tolerance
+    character(len=*), intent(inout) :: mesh_file_name
+
+    character(len=LOG_SIZE) :: log_buf
+    integer :: sample_idx, sample_mesh_idx
+    integer :: last_index
+    type(fld_file_data_t) :: fld_data
+    type(file_t) :: f
+    logical :: mesh_mismatch
+
+    ! ---- For the mesh to mesh interpolation
+    type(global_interpolation_t) :: global_interp
+    ! -----
+
+    ! ---- For space to space interpolation
+    type(space_t) :: prev_Xh
+    type(interpolator_t) :: space_interp
+    ! ----
+
+    call neko_log%message("Type          : field")
+    call neko_log%message("File name     : " // trim(file_name))
+    write (log_buf, '(A,L1)') "Interpolation : ", interpolate
+    call neko_log%message(log_buf)
+    if (interpolate) then
+    end if
+
+    ! Extract sample index from the file name
+    sample_idx = extract_fld_file_index(file_name, -1)
+
+    if (sample_idx .eq. -1) &
+         call neko_error("Invalid file name for the initial condition. The&
+         & file format must be e.g. 'mean0.f00001'")
+
+    ! Change from "field0.f000*" to "field0.fld" for the fld reader
+    call filename_chsuffix(file_name, file_name, 'fld')
+
+    call fld_data%init
+    f = file_t(trim(file_name))
+
+    if (interpolate) then
+
+       ! If no mesh file is specified, use the default file name
+       if (mesh_file_name .eq. "none") then
+          mesh_file_name = trim(file_name)
+          sample_mesh_idx = sample_idx
+       else
+
+          ! Extract sample index from the mesh file name
+          sample_mesh_idx = extract_fld_file_index(mesh_file_name, -1)
+
+          if (sample_mesh_idx .eq. -1) &
+               call neko_error("Invalid file name for the initial condition. &
+&The file format must be e.g. 'mean0.f00001'")
+
+          write (log_buf, '(A,ES12.6)') "Tolerance     : ", tolerance
+          call neko_log%message(log_buf)
+          write (log_buf, '(A,A)')     "Mesh file     : ", &
+               trim(mesh_file_name)
+          call neko_log%message(log_buf)
+
+       end if ! if mesh_file_name .eq. none
+
+       ! Read the mesh coordinates if they are not in our fld file
+       if (sample_mesh_idx .ne. sample_idx) then
+          call f%set_counter(sample_mesh_idx)
+          call f%read(fld_data)
+       end if
+
+    end if
+
+    ! Read the field file containing (u,v,w,p)
+    call f%set_counter(sample_idx)
+    call f%read(fld_data)
+
+    !
+    ! Check that the data in the fld file matches the current case.
+    ! Note that this is a safeguard and there are corner cases where
+    ! two different meshes have the same dimension and same # of elements
+    ! but this should be enough to cover obvious cases.
+    !
+    mesh_mismatch = (fld_data%glb_nelv .ne. u%msh%glb_nelv .or. &
+         fld_data%gdim .ne. u%msh%gdim)
+
+    if (mesh_mismatch .and. .not. interpolate) then
+       call neko_error("The fld file must match the current mesh! &
+&Use 'interpolate': 'true' to enable interpolation.")
+    else if (.not. mesh_mismatch .and. interpolate) then
+       call neko_log%warning("You have activated interpolation but you might &
+&still be using the same mesh.")
+    end if
+
+    ! Mesh interpolation if specified
+    if (interpolate) then
+
+       ! Issue a warning if the mesh is in single precision
+       select type (ft => f%file_type)
+       type is (fld_file_t)
+          if (.not. ft%dp_precision) then
+             call neko_warning("The coordinates read from the field file are &
+&in single precision.")
+             call neko_log%message("It is recommended to use a mesh in double &
+&precision for better interpolation results.")
+             call neko_log%message("If the interpolation does not work, you&
+&can try to increase the tolerance.")
+          end if
+       class default
+       end select
+
+       ! Generates an interpolator object and performs the point search
+       global_interp = fld_data%generate_interpolator(u%dof, u%msh, &
+            tolerance)
+
+       ! Evaluate velocities and pressure
+       call global_interp%evaluate(u%x, fld_data%u%x)
+       call global_interp%evaluate(v%x, fld_data%v%x)
+       call global_interp%evaluate(w%x, fld_data%w%x)
+       call global_interp%evaluate(p%x, fld_data%p%x)
+
+       call global_interp%free
+
+    else ! No interpolation, but potentially just from different spaces
+
+       ! Build a space_t object from the data in the fld file
+       call prev_Xh%init(GLL, fld_data%lx, fld_data%ly, fld_data%lz)
+       call space_interp%init(u%Xh, prev_Xh)
+
+       ! Do the space-to-space interpolation
+       call space_interp%map_host(u%x, fld_data%u%x, fld_data%nelv, u%Xh)
+       call space_interp%map_host(v%x, fld_data%v%x, fld_data%nelv, u%Xh)
+       call space_interp%map_host(w%x, fld_data%w%x, fld_data%nelv, u%Xh)
+       call space_interp%map_host(p%x, fld_data%p%x, fld_data%nelv, u%Xh)
+
+       call space_interp%free
+
+    end if
+
+    ! NOTE: we do not copy (u,v,w) to the GPU since `set_flow_ic_common` does
+    ! the copy for us
+    if (NEKO_BCKND_DEVICE .eq. 1) call device_memcpy(p%x, p%x_d, p%dof%size(), &
+         HOST_TO_DEVICE, sync = .false.)
+
+    call fld_data%free
+
+  end subroutine set_flow_ic_fld
+
 end module flow_ic
diff --git a/src/io/fld_file_data.f90 b/src/io/fld_file_data.f90
index 4721de827a6..16ddd87db5e 100644
--- a/src/io/fld_file_data.f90
+++ b/src/io/fld_file_data.f90
@@ -2,12 +2,19 @@
 !! Provides an interface to the different fields sotred in a fld file
 !! Also provides simple functions to scale and add different fld files.
 !! An example of using this module is shown in contrib/average_fields.f90
-!! The fld_file_data_t should dynamically update each time one reads a new fld file
+!! The fld_file_data_t should dynamically update each time one reads a new fld
+!! file.
 !! Martin Karp 1/2-2023
 module fld_file_data
   use num_types, only : rp
-  use math
+  use math, only: cmult, add2
   use vector, only : vector_t, vector_ptr_t
+  use field, only: field_t
+  use dofmap, only: dofmap_t
+  use space, only: space_t, GLL
+  use global_interpolation, only: global_interpolation_t
+  use utils, only: neko_error
+  use mesh, only: mesh_t
   implicit none
   private
 
@@ -33,18 +40,23 @@ module fld_file_data
      integer :: lz = 0
      integer :: t_counter = 0 !< counter of samples
      ! meta file information (if any)
-     integer :: meta_nsamples = 0 !< number of samples specified in .nek5000 file
+     !> number of samples specified in .nek5000 file
+     integer :: meta_nsamples = 0
      integer :: meta_start_counter = 0 !< number of first field
-     character(len=1024) :: fld_series_fname !< name of fld series as specified in .nek5000 (meta) file
+     !> name of fld series as specified in .nek5000 (meta) file
+     character(len=1024) :: fld_series_fname
    contains
      procedure, pass(this) :: init => fld_file_data_init
      procedure, pass(this) :: free => fld_file_data_free
      procedure, pass(this) :: scale => fld_file_data_scale
      procedure, pass(this) :: add => fld_file_data_add
      procedure, pass(this) :: size => fld_file_data_size
-     procedure, pass(this) :: get_list => fld_file_get_list
-     procedure, pass(this) :: init_same => fld_file_init_same
-     procedure, pass(this) :: init_n_fields => fld_file_init_n_fields
+     procedure, pass(this) :: get_list => fld_file_data_get_list
+     procedure, pass(this) :: init_same => fld_file_data_init_same
+     procedure, pass(this) :: init_n_fields => fld_file_data_init_n_fields
+     !> Generates a global_interpolation object to interpolate the fld data.
+     procedure, pass(this) :: generate_interpolator => &
+          fld_file_data_generate_interpolator
   end type fld_file_data_t
 
 contains
@@ -53,6 +65,7 @@ module fld_file_data
   subroutine fld_file_data_init(this, nelv, offset_el)
     class(fld_file_data_t), intent(inout) :: this
     integer, intent(in), optional :: nelv, offset_el
+
     call this%free()
     if (present(nelv)) this%nelv = nelv
     if (present(offset_el)) this%offset_el = offset_el
@@ -63,17 +76,17 @@ function fld_file_data_size(this) result(i)
     class(fld_file_data_t) :: this
     integer :: i
     i = 0
-    if(this%u%n .gt. 0) i = i + 1
-    if(this%v%n .gt. 0) i = i + 1
-    if(this%w%n .gt. 0) i = i + 1
-    if(this%p%n .gt. 0) i = i + 1
-    if(this%t%n .gt. 0) i = i + 1
+    if (this%u%n .gt. 0) i = i + 1
+    if (this%v%n .gt. 0) i = i + 1
+    if (this%w%n .gt. 0) i = i + 1
+    if (this%p%n .gt. 0) i = i + 1
+    if (this%t%n .gt. 0) i = i + 1
     i = i + this%n_scalars
 
   end function fld_file_data_size
 
   !> Genereate same fields as in another fld_file
-  subroutine fld_file_init_same(this, fld_file, n)
+  subroutine fld_file_data_init_same(this, fld_file, n)
     class(fld_file_data_t), target, intent(inout) :: this
     class(fld_file_data_t), target, intent(in) :: fld_file
     integer, intent(in) :: n
@@ -100,10 +113,10 @@ subroutine fld_file_init_same(this, fld_file, n)
        call this%s(j)%init(n)
     end do
 
-  end subroutine fld_file_init_same
+  end subroutine fld_file_data_init_same
 
   !> Genereate same fields as in another fld_file
-  subroutine fld_file_init_n_fields(this, n_fields, n)
+  subroutine fld_file_data_init_n_fields(this, n_fields, n)
     class(fld_file_data_t), target, intent(inout) :: this
     integer, intent(in) :: n, n_fields
     integer :: i, j
@@ -132,35 +145,32 @@ subroutine fld_file_init_n_fields(this, n_fields, n)
        end do
     end if
 
-  end subroutine fld_file_init_n_fields
-
-
-
+  end subroutine fld_file_data_init_n_fields
 
   !> Get a list with pointers to the fields in the fld file
-  subroutine fld_file_get_list(this, ptr_list, n)
+  subroutine fld_file_data_get_list(this, ptr_list, n)
     class(fld_file_data_t), target, intent(in) :: this
     integer, intent(in) :: n
     integer :: i, j
     type(vector_ptr_t), intent(inout) :: ptr_list(n)
     i = 1
-    if(this%u%n .gt. 0) then
+    if (this%u%n .gt. 0) then
        ptr_list(i)%ptr => this%u
        i = i + 1
     end if
-    if(this%v%n .gt. 0) then
+    if (this%v%n .gt. 0) then
        ptr_list(i)%ptr => this%v
        i = i + 1
     end if
-    if(this%w%n .gt. 0) then
+    if (this%w%n .gt. 0) then
        ptr_list(i)%ptr => this%w
        i = i + 1
     end if
-    if(this%p%n .gt. 0) then
+    if (this%p%n .gt. 0) then
        ptr_list(i)%ptr => this%p
        i = i + 1
     end if
-    if(this%t%n .gt. 0) then
+    if (this%t%n .gt. 0) then
        ptr_list(i)%ptr => this%t
        i = i + 1
     end if
@@ -169,7 +179,7 @@ subroutine fld_file_get_list(this, ptr_list, n)
        i = i +1
     end do
 
-  end subroutine fld_file_get_list
+  end subroutine fld_file_data_get_list
 
 
 
@@ -179,14 +189,14 @@ subroutine fld_file_data_scale(this, c)
     real(kind=rp), intent(in) :: c
     integer :: i
 
-    if(this%u%n .gt. 0) call cmult(this%u%x,c,this%u%n)
-    if(this%v%n .gt. 0) call cmult(this%v%x,c,this%v%n)
-    if(this%w%n .gt. 0) call cmult(this%w%x,c,this%w%n)
-    if(this%p%n .gt. 0) call cmult(this%p%x,c,this%p%n)
-    if(this%t%n .gt. 0) call cmult(this%t%x,c,this%t%n)
+    if (this%u%n .gt. 0) call cmult(this%u%x, c, this%u%n)
+    if (this%v%n .gt. 0) call cmult(this%v%x, c, this%v%n)
+    if (this%w%n .gt. 0) call cmult(this%w%x, c, this%w%n)
+    if (this%p%n .gt. 0) call cmult(this%p%x, c, this%p%n)
+    if (this%t%n .gt. 0) call cmult(this%t%x, c, this%t%n)
 
     do i = 1, this%n_scalars
-       if(this%s(i)%n .gt. 0) call cmult(this%s(i)%x,c,this%s(i)%n)
+       if (this%s(i)%n .gt. 0) call cmult(this%s(i)%x, c, this%s(i)%n)
     end do
 
   end subroutine fld_file_data_scale
@@ -197,14 +207,15 @@ subroutine fld_file_data_add(this, fld_data_add)
     class(fld_file_data_t), intent(in) :: fld_data_add
     integer :: i
 
-    if(this%u%n .gt. 0) call add2(this%u%x,fld_data_add%u%x,this%u%n)
-    if(this%v%n .gt. 0) call add2(this%v%x,fld_data_add%v%x,this%v%n)
-    if(this%w%n .gt. 0) call add2(this%w%x,fld_data_add%w%x,this%w%n)
-    if(this%p%n .gt. 0) call add2(this%p%x,fld_data_add%p%x,this%p%n)
-    if(this%t%n .gt. 0) call add2(this%t%x,fld_data_add%t%x,this%t%n)
+    if (this%u%n .gt. 0) call add2(this%u%x, fld_data_add%u%x, this%u%n)
+    if (this%v%n .gt. 0) call add2(this%v%x, fld_data_add%v%x, this%v%n)
+    if (this%w%n .gt. 0) call add2(this%w%x, fld_data_add%w%x, this%w%n)
+    if (this%p%n .gt. 0) call add2(this%p%x, fld_data_add%p%x, this%p%n)
+    if (this%t%n .gt. 0) call add2(this%t%x, fld_data_add%t%x, this%t%n)
 
     do i = 1, this%n_scalars
-       if(this%s(i)%n .gt. 0) call add2(this%s(i)%x,fld_data_add%s(i)%x,this%s(i)%n)
+       if (this%s(i)%n .gt. 0) call add2(this%s(i)%x, fld_data_add%s(i)%x, &
+            this%s(i)%n)
     end do
   end subroutine fld_file_data_add
 
@@ -240,4 +251,82 @@ subroutine fld_file_data_free(this)
     if(allocated(this%idx)) deallocate(this%idx)
   end subroutine fld_file_data_free
 
+  !> Generates a global_interpolation object to interpolate the fld data.
+  !! @param to_dof Dofmap on which to interpolate.
+  !! @param to_msh Mesh on which to interpolate.
+  !! @param tolerance Tolerance for the newton iterations.
+  function fld_file_data_generate_interpolator(this, to_dof, &
+       to_msh, tolerance) result(global_interp)
+    class(fld_file_data_t), intent(in) :: this
+    type(dofmap_t), intent(in), target :: to_dof
+    type(mesh_t), intent(in), target :: to_msh
+    real(kind=rp), intent(in) :: tolerance
+
+    type(global_interpolation_t) :: global_interp
+
+    ! --- variables for interpolation
+    type(space_t) :: fld_Xh
+    real(kind=rp), allocatable :: x_coords(:,:,:,:), y_coords(:,:,:,:), &
+         z_coords(:,:,:,:)
+    real(kind=rp) :: center_x,  center_y, center_z
+    integer :: e, i
+    ! ---
+
+    type(space_t), pointer :: to_Xh
+    to_Xh => to_dof%Xh
+
+    ! Safeguard in case we didn't read mesh information
+    if (.not. allocated(this%x%x) .or. &
+         .not. allocated(this%y%x) .or. &
+         .not. allocated(this%z%x)) call neko_error("Unable to retrieve &
+&mesh information from fld data.")
+
+    ! Create a space based on the fld data
+    call fld_Xh%init(GLL, this%lx,  this%ly, this%lz)
+
+    ! These are the coordinates of our current dofmap
+    ! that we use for the interpolation
+    allocate(x_coords(to_Xh%lx,  to_Xh%ly, to_Xh%lz, to_msh%nelv))
+    allocate(y_coords(to_Xh%lx,  to_Xh%ly, to_Xh%lz, to_msh%nelv))
+    allocate(z_coords(to_Xh%lx,  to_Xh%ly, to_Xh%lz, to_msh%nelv))
+
+    !> To ensure that each point is within an element
+    !! Remedies issue with points on the boundary
+    !! Technically gives each point a slightly different value
+    !! but still within the specified tolerance
+    do e = 1, to_msh%nelv
+       center_x = 0d0
+       center_y = 0d0
+       center_z = 0d0
+       do i = 1, to_Xh%lxyz
+          center_x = center_x + to_dof%x(i, 1, 1, e)
+          center_y = center_y + to_dof%y(i, 1, 1, e)
+          center_z = center_z + to_dof%z(i, 1, 1, e)
+       end do
+       center_x = center_x / to_Xh%lxyz
+       center_y = center_y / to_Xh%lxyz
+       center_z = center_z / to_Xh%lxyz
+       do i = 1, to_Xh%lxyz
+          x_coords(i, 1, 1, e) = to_dof%x(i, 1, 1, e) - &
+               tolerance * (to_dof%x(i, 1, 1, e) - center_x)
+          y_coords(i, 1, 1, e) = to_dof%y(i, 1, 1, e) - &
+               tolerance * (to_dof%y(i, 1, 1, e) - center_y)
+          z_coords(i, 1, 1, e) = to_dof%z(i, 1, 1, e) - &
+               tolerance * (to_dof%z(i, 1, 1, e) - center_z)
+       end do
+    end do
+
+    ! The initialization is done based on the variables created from
+    ! fld data
+    call global_interp%init(this%x%x, this%y%x, this%z%x, this%gdim, &
+         this%nelv, fld_Xh, tol = tolerance)
+    call global_interp%find_points(x_coords, y_coords, z_coords, &
+         to_dof%size())
+
+    deallocate(x_coords)
+    deallocate(y_coords)
+    deallocate(z_coords)
+
+  end function fld_file_data_generate_interpolator
+
 end module fld_file_data
diff --git a/src/scalar/scalar_ic.f90 b/src/scalar/scalar_ic.f90
index 644bbf4f973..1ad26417943 100644
--- a/src/scalar/scalar_ic.f90
+++ b/src/scalar/scalar_ic.f90
@@ -38,14 +38,24 @@ module scalar_ic
   use device_math, only : device_col2
   use device, only : device_memcpy, HOST_TO_DEVICE
   use field, only : field_t
-  use utils, only : neko_error
+  use utils, only : neko_error, filename_chsuffix, filename_suffix, &
+       neko_warning, NEKO_FNAME_LEN, extract_fld_file_index
   use coefs, only : coef_t
   use math, only : col2, cfill, cfill_mask
   use user_intf, only : useric_scalar
   use json_module, only : json_file
-  use json_utils, only: json_get
+  use json_utils, only: json_get, json_get_or_default
   use point_zone, only: point_zone_t
   use point_zone_registry, only: neko_point_zone_registry
+  use field_registry, only: neko_field_registry
+  use logger, only: neko_log, LOG_SIZE
+  use fld_file_data, only: fld_file_data_t
+  use fld_file, only: fld_file_t
+  use checkpoint, only: chkp_t
+  use file, only: file_t
+  use global_interpolation, only: global_interpolation_t
+  use interpolation, only: interpolator_t
+  use space, only: space_t, GLL
   implicit none
   private
 
@@ -60,7 +70,9 @@ module scalar_ic
   !> Set scalar initial condition (builtin)
   !! @details Set scalar initial condition using one of the builtin types
   !! currently supported:
-  !! - uniform.
+  !! - uniform
+  !! - point zone
+  !! - field
   !! @param s Scalar field.
   !! @param coef Coefficient.
   !! @param gs Gather-Scatter object.
@@ -75,20 +87,44 @@ subroutine set_scalar_ic_int(s, coef, gs, type, params)
 
     ! Variables for retrieving JSON parameters
     real(kind=rp) :: ic_value
-    character(len=:), allocatable :: zone_name
-    real(kind=rp) :: zone_value
+    character(len=:), allocatable :: read_str
+    character(len=NEKO_FNAME_LEN) :: fname, mesh_fname
+    real(kind=rp) :: zone_value, tol
+    logical :: interpolate
 
     if (trim(type) .eq. 'uniform') then
+
        call json_get(params, 'case.scalar.initial_condition.value', ic_value)
        call set_scalar_ic_uniform(s, ic_value)
+
     else if (trim(type) .eq. 'point_zone') then
+
        call json_get(params, 'case.scalar.initial_condition.base_value', &
-                     ic_value)
+            ic_value)
        call json_get(params, 'case.scalar.initial_condition.zone_name', &
-                     zone_name)
+            read_str)
        call json_get(params, 'case.scalar.initial_condition.zone_value', &
-                     zone_value)
-       call set_scalar_ic_point_zone(s, ic_value, zone_name, zone_value)
+            zone_value)
+
+       call set_scalar_ic_point_zone(s, ic_value, read_str, zone_value)
+
+    else if (trim(type) .eq. 'field') then
+
+       call json_get(params, 'case.scalar.initial_condition.file_name', &
+            read_str)
+       fname = trim(read_str)
+       call json_get_or_default(params, &
+            'case.scalar.initial_condition.interpolate', interpolate, &
+            .false.)
+       call json_get_or_default(params, &
+            'case.scalar.initial_condition.tolerance', tol, 0.000001_rp)
+       call json_get_or_default(params, &
+            'case.scalar.initial_condition.mesh_file_name', read_str, &
+            "none")
+       mesh_fname = trim(read_str)
+
+       call set_scalar_ic_fld(s, fname, interpolate, tol, mesh_fname)
+
     else
        call neko_error('Invalid initial condition')
     end if
@@ -111,6 +147,7 @@ subroutine set_scalar_ic_usr(s, coef, gs, usr_ic, params)
     procedure(useric_scalar) :: usr_ic
     type(json_file), intent(inout) :: params
 
+    call neko_log%message("Type: user")
     call usr_ic(s, params)
 
     call set_scalar_ic_common(s, coef, gs)
@@ -132,7 +169,7 @@ subroutine set_scalar_ic_common(s, coef, gs)
     n = s%dof%size()
     if (NEKO_BCKND_DEVICE .eq. 1) then
        call device_memcpy(s%x, s%x_d, n, &
-                          HOST_TO_DEVICE, sync=.false.)
+                          HOST_TO_DEVICE, sync = .false.)
     end if
 
     ! Ensure continuity across elements for initial conditions
@@ -154,6 +191,12 @@ subroutine set_scalar_ic_uniform(s, ic_value)
     type(field_t), intent(inout) :: s
     real(kind=rp), intent(in) :: ic_value
     integer :: n
+    character(len=LOG_SIZE) :: log_buf
+
+    call neko_log%message("Type : uniform")
+    write (log_buf, '(A,ES12.6)') "Value: ", ic_value
+    call neko_log%message(log_buf)
+
     s = ic_value
     n = s%dof%size()
     if (NEKO_BCKND_DEVICE .eq. 1) then
@@ -171,14 +214,22 @@ end subroutine set_scalar_ic_uniform
   !! @param zone_value Desired value of the scalar field in the point zone.
   subroutine set_scalar_ic_point_zone(s, base_value, zone_name, zone_value)
     type(field_t), intent(inout) :: s
-    real(kind=rp), intent(in):: base_value
+    real(kind=rp), intent(in) :: base_value
     character(len=*), intent(in) :: zone_name
     real(kind=rp), intent(in) :: zone_value
 
     ! Internal variables
+    character(len=LOG_SIZE) :: log_buf
     class(point_zone_t), pointer :: zone
     integer :: size
 
+    call neko_log%message("Type       : point_zone")
+    write (log_buf, '(A,ES12.6)') "Base value: ", base_value
+    call neko_log%message(log_buf)
+    call neko_log%message("Zone name : " // trim(zone_name))
+    write (log_buf, '(A,ES12.6)') "Zone value: ", zone_value
+    call neko_log%message(log_buf)
+
     size = s%dof%size()
     zone => neko_point_zone_registry%get_point_zone(trim(zone_name))
 
@@ -187,4 +238,155 @@ subroutine set_scalar_ic_point_zone(s, base_value, zone_name, zone_value)
 
   end subroutine set_scalar_ic_point_zone
 
+  !> Set the initial condition of the scalar based on a field.
+  !! @detail The field is read from an `fld` file. If enabled, interpolation
+  !! is also possible. In that case, the mesh coordinates can be read from
+  !! another file in the `fld` field series.
+  !! @param s The scalar field.
+  !! @param file_name The name of the "fld" file series.
+  !! @param sample_idx index of the field file .f000* to read, default is
+  !! -1.
+  !! @param interpolate Flag to indicate wether or not to interpolate the
+  !! values onto the current mesh.
+  !! @param tolerance If interpolation is enabled, tolerance for finding the
+  !! points in the mesh.
+  !! @param sample_mesh_idx If interpolation is enabled, index of the field
+  !! file where the mesh coordinates are located.
+  subroutine set_scalar_ic_fld(s, file_name, &
+       interpolate, tolerance, mesh_file_name)
+    type(field_t), intent(inout) :: s
+    character(len=*), intent(in) :: file_name
+    logical, intent(in) :: interpolate
+    real(kind=rp), intent(in) :: tolerance
+    character(len=*), intent(inout) :: mesh_file_name
+
+    character(len=LOG_SIZE) :: log_buf
+    integer :: sample_idx, sample_mesh_idx
+    integer :: last_index
+    type(fld_file_data_t) :: fld_data
+    type(file_t) :: f
+    logical :: mesh_mismatch
+
+    ! ---- For the mesh to mesh interpolation
+    type(global_interpolation_t) :: global_interp
+    ! -----
+
+    ! ---- For space to space interpolation
+    type(space_t) :: prev_Xh
+    type(interpolator_t) :: space_interp
+    ! ----
+
+    call neko_log%message("Type          : field")
+    call neko_log%message("File name     : " // trim(file_name))
+    write (log_buf, '(A,L1)') "Interpolation : ", interpolate
+    call neko_log%message(log_buf)
+    if (interpolate) then
+    end if
+
+    ! Extract sample index from the file name
+    sample_idx = extract_fld_file_index(file_name, -1)
+
+    if (sample_idx .eq. -1) &
+         call neko_error("Invalid file name for the initial condition. The&
+         & file format must be e.g. 'mean0.f00001'")
+
+    ! Change from "field0.f000*" to "field0.fld" for the fld reader
+    call filename_chsuffix(file_name, file_name, 'fld')
+
+    call fld_data%init
+    f = file_t(trim(file_name))
+
+    if (interpolate) then
+
+       ! If no mesh file is specified, use the default file name
+       if (mesh_file_name .eq. "none") then
+          mesh_file_name = trim(file_name)
+          sample_mesh_idx = sample_idx
+       else
+
+          ! Extract sample index from the mesh file name
+          sample_mesh_idx = extract_fld_file_index(mesh_file_name, -1)
+
+          if (sample_mesh_idx .eq. -1) &
+               call neko_error("Invalid file name for the initial condition. &
+&The file format must be e.g. 'mean0.f00001'")
+
+          write (log_buf, '(A,ES12.6)') "Tolerance     : ", tolerance
+          call neko_log%message(log_buf)
+          write (log_buf, '(A,A)')     "Mesh file     : ", &
+               trim(mesh_file_name)
+          call neko_log%message(log_buf)
+
+       end if ! if mesh_file_name .eq. none
+
+       ! Read the mesh coordinates if they are not in our fld file
+       if (sample_mesh_idx .ne. sample_idx) then
+          call f%set_counter(sample_mesh_idx)
+          call f%read(fld_data)
+       end if
+
+    end if
+
+    ! Read the field file containing (u,v,w,p)
+    call f%set_counter(sample_idx)
+    call f%read(fld_data)
+
+    !
+    ! Check that the data in the fld file matches the current case.
+    ! Note that this is a safeguard and there are corner cases where
+    ! two different meshes have the same dimension and same # of elements
+    ! but this should be enough to cover obvious cases.
+    !
+    mesh_mismatch = (fld_data%glb_nelv .ne. s%msh%glb_nelv .or. &
+         fld_data%gdim .ne. s%msh%gdim)
+
+    if (mesh_mismatch .and. .not. interpolate) then
+       call neko_error("The fld file must match the current mesh! &
+&Use 'interpolate': 'true' to enable interpolation.")
+    else if (.not. mesh_mismatch .and. interpolate) then
+       call neko_log%warning("You have activated interpolation but you might &
+&still be using the same mesh.")
+    end if
+
+
+    ! Mesh interpolation if specified
+    if (interpolate) then
+       ! Issue a warning if the mesh is in single precision
+       select type (ft => f%file_type)
+       type is (fld_file_t)
+          if (.not. ft%dp_precision) then
+             call neko_warning("The coordinates read from the field file are &
+&in single precision.")
+             call neko_log%message("It is recommended to use a mesh in double &
+&precision for better interpolation results.")
+             call neko_log%message("If the interpolation does not work, you&
+&can try to increase the tolerance.")
+          end if
+       class default
+       end select
+
+       ! Generates an interpolator object and performs the point search
+       global_interp = fld_data%generate_interpolator(s%dof, s%msh, tolerance)
+
+       ! Evaluate scalar
+       call global_interp%evaluate(s%x, fld_data%t%x)
+       call global_interp%free
+
+    else ! No interpolation, just potentially from different spaces
+
+       ! Build a space_t object from the data in the fld file
+       call prev_Xh%init(GLL, fld_data%lx, fld_data%ly, fld_data%lz)
+       call space_interp%init(s%Xh, prev_Xh)
+
+       ! Do the space-to-space interpolation
+       call space_interp%map_host(s%x, fld_data%t%x, fld_data%nelv, s%Xh)
+
+       call space_interp%free
+
+    end if
+
+    call fld_data%free
+
+  end subroutine set_scalar_ic_fld
+
 end module scalar_ic
diff --git a/src/simulation_components/probes.F90 b/src/simulation_components/probes.F90
index 1755781c33b..aea070a99e0 100644
--- a/src/simulation_components/probes.F90
+++ b/src/simulation_components/probes.F90
@@ -179,7 +179,7 @@ subroutine probes_init_from_json(this, json, case)
 
          case ('file')
           call this%read_file(json_point)
-         case ('point')
+         case ('points')
           call this%read_point(json_point)
          case ('line')
           call this%read_line(json_point)
@@ -187,10 +187,8 @@ subroutine probes_init_from_json(this, json, case)
           call neko_error('Plane probes not implemented yet.')
          case ('circle')
           call this%read_circle(json_point)
-
          case ('point_zone')
           call this%read_point_zone(json_point, case%fluid%dm_Xh)
-
          case ('none')
           call json_point%print()
           call neko_error('No point type specified.')
@@ -241,17 +239,11 @@ subroutine read_point(this, json)
 
     real(kind=rp), dimension(:,:), allocatable :: point_list
     real(kind=rp), dimension(:), allocatable :: rp_list_reader
-    logical :: found
 
     ! Ensure only rank 0 reads the coordinates.
     if (pe_rank .ne. 0) return
     call json_get(json, 'coordinates', rp_list_reader)
 
-    ! Check if the coordinates were found and were valid
-    if (.not. found) then
-       call neko_error('No coordinates found.')
-    end if
-
     if (mod(size(rp_list_reader), 3) /= 0) then
        call neko_error('Invalid number of coordinates.')
     end if