fld_file.f90

Go to the documentation of this file.

! Copyright (c) 2020-2023, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module fld_file
  use num_types, only : rp, dp, sp, i8
  use generic_file, only : generic_file_t
  use field, only : field_t
  use field_list, only : field_list_t
  use dofmap, only : dofmap_t
  use space, only : space_t
  use structs, only : array_ptr_t
  use vector, only : vector_t
  use fld_file_data, only : fld_file_data_t
  use vector, only : vector_t
  use space, only : space_t
  use logger, only : neko_log, log_size
  use mesh, only : mesh_t
  use mask, only : mask_t
  use utils, only : filename_suffix_pos, filename_chsuffix, filename_name, &
       filename_path, neko_error, nonlinear_index
  use comm
  use datadist, only : linear_dist_t
  use math, only : vlmin, vlmax, sabscmp
  use neko_mpi_types, only : mpi_character_size, mpi_double_precision_size, &
       mpi_real_size, mpi_integer_size
  use device, only : device_sync, host_to_device
  use mpi_f08
  implicit none
  private
 
  real(kind=dp), private, allocatable :: tmp_dp(:)
  real(kind=sp), private, allocatable :: tmp_sp(:)
 
  type, public, extends(generic_file_t) :: fld_file_t
     logical :: dp_precision = .false. 
     ! These flags can be manipulated to skip writing certain fields
     ! that will then be put into scalars.
     logical :: skip_pressure = .true. 
     logical :: skip_velocity = .true. 
     logical :: skip_temperature = .true. 
     logical :: write_mesh = .false. 
     type(mask_t) :: mask
   contains
     procedure :: read => fld_file_read
     procedure :: write => fld_file_write_manager
     procedure, private :: write_all => fld_file_write
     procedure, private :: write_masked => fld_file_write_masked
     procedure :: set_precision => fld_file_set_precision
     procedure :: set_mask => fld_file_set_mask
     procedure :: get_fld_fname => fld_file_get_fld_fname
     procedure :: get_meta_fname => fld_file_get_meta_fname
  end type fld_file_t
 
 
contains
 
  subroutine fld_file_select_from_field_list(this, data, p, u, v, w, tem, &
       scalar_fields, n_scalar_fields, write_pressure, write_velocity, &
       write_temperature, dof)
    class(fld_file_t), intent(in) :: this
    type(field_list_t), intent(in) :: data
    type(array_ptr_t), intent(inout) :: p, u, v, w, tem
    type(array_ptr_t), allocatable, intent(inout) :: scalar_fields(:)
    integer, intent(out) :: n_scalar_fields
    logical, intent(out) :: write_pressure, write_velocity, write_temperature
    type(dofmap_t), pointer, intent(out) :: dof
    integer :: i, idx, n_fields
 
    n_scalar_fields = 0
    write_pressure = .false.
    write_velocity = .false.
    write_temperature = .false.
    n_fields = data%size()
    idx = 1
 
    if (n_fields .eq. 0) then
       call neko_error('Empty field_list_t cannot be written to an fld file')
    end if
 
    if (.not. this%skip_pressure) then
       if (idx .le. n_fields) then
          p%ptr => data%items(idx)%ptr%x(:,1,1,1)
          write_pressure = .true.
          idx = idx + 1
       end if
    end if
 
    if (.not. this%skip_velocity) then
       if (idx + 2 .le. n_fields) then
          u%ptr => data%items(idx+0)%ptr%x(:,1,1,1)
          v%ptr => data%items(idx+1)%ptr%x(:,1,1,1)
          w%ptr => data%items(idx+2)%ptr%x(:,1,1,1)
          write_velocity = .true.
          idx = idx + 3
       end if
    end if
 
    if (.not. this%skip_temperature) then
       if (idx .le. n_fields) then
          tem%ptr => data%items(idx)%ptr%x(:,1,1,1)
          write_temperature = .true.
          idx = idx + 1
       end if
    end if
 
    if (idx .le. n_fields) then
       n_scalar_fields = n_fields - idx + 1
       allocate(scalar_fields(n_scalar_fields))
       do i = 1, n_scalar_fields
          scalar_fields(i)%ptr => data%items(idx + i - 1)%ptr%x(:,1,1,1)
       end do
    end if
 
    dof => data%dof(1)
  end subroutine fld_file_select_from_field_list
 
  subroutine fld_file_write_manager(this, data, t)
    class(fld_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    real(kind=rp), intent(in), optional :: t
 
    if (this%mask%is_set()) then
       call this%write_masked(data, this%mask, t)
    else
       call this%write_all(data, t)
    end if
 
  end subroutine fld_file_write_manager
 
  subroutine fld_file_write(this, data, t)
    class(fld_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    real(kind=rp), intent(in), optional :: t
    type(array_ptr_t) :: x, y, z, u, v, w, p, tem
    real(kind=rp), allocatable, target :: tempo(:)
    type(mesh_t), pointer :: msh
    type(dofmap_t), pointer :: dof
    type(space_t), pointer :: Xh
    real(kind=dp) :: time
    character(len= 132) :: hdr
    character :: rdcode(10)
    character(len=6) :: id_str
    character(len= 1024) :: fname
    character(len= 1024) :: name
    integer :: file_unit
    integer :: i, ierr, n, suffix_pos, tslash_pos
    integer :: lx, ly, lz, lxyz, gdim, glb_nelv, nelv, offset_el
    integer, allocatable :: idx(:)
    type(mpi_status) :: status
    type(mpi_file) :: fh
    integer (kind=MPI_OFFSET_KIND) :: mpi_offset, byte_offset, temp_offset
    real(kind=sp), parameter :: test_pattern = 6.54321
    type(array_ptr_t), allocatable :: scalar_fields(:)
    logical :: write_mesh, write_velocity, write_pressure, write_temperature
    integer :: fld_data_size, n_scalar_fields
    type(field_list_t) :: dummy_list
    type(field_t), pointer :: ptr
 
    if (present(t)) then
       time = real(t, dp)
    else
       time = 0d0
    end if
 
    nullify(msh)
    nullify(dof)
    nullify(xh)
    n_scalar_fields = 0
    write_pressure = .false.
    write_velocity = .false.
    write_temperature = .false.
 
    select type (data)
    type is (fld_file_data_t)
       nelv = data%nelv
       lx = data%lx
       ly = data%ly
       lz = data%lz
       gdim = data%gdim
       glb_nelv = data%glb_nelv
       offset_el = data%offset_el
 
       if (data%x%size() .gt. 0) x%ptr => data%x%x
       if (data%y%size() .gt. 0) y%ptr => data%y%x
       if (data%z%size() .gt. 0) z%ptr => data%z%x
       if (gdim .eq. 2) z%ptr => data%y%x
       if (data%u%size() .gt. 0) then
          u%ptr => data%u%x
          ! In case only u is actually allocated, point the other comps to u
          ! so that we don't die on trying to write them
          if (data%v%size() .le. 0) v%ptr => data%u%x
          if (data%w%size() .le. 0) w%ptr => data%u%x
          write_velocity = .true.
       end if
       if (data%v%size() .gt. 0) v%ptr => data%v%x
       if (data%w%size() .gt. 0) w%ptr => data%w%x
       if (data%p%size() .gt. 0) then
          p%ptr => data%p%x
          write_pressure = .true.
       end if
       if (data%t%size() .gt. 0) then
          write_temperature = .true.
          tem%ptr => data%t%x
       end if
       ! If gdim = 2 and Z-velocity component exists,
       ! it is stored in last scalar field
       if (gdim .eq. 2 .and. data%w%size() .gt. 0) then
          n_scalar_fields = data%n_scalars + 1
          allocate(scalar_fields(n_scalar_fields))
          do i = 1, n_scalar_fields -1
             scalar_fields(i)%ptr => data%s(i)%x
          end do
          scalar_fields(n_scalar_fields)%ptr => data%w%x
       else
          n_scalar_fields = data%n_scalars
          allocate(scalar_fields(n_scalar_fields+1))
          do i = 1, n_scalar_fields
             scalar_fields(i)%ptr => data%s(i)%x
          end do
          scalar_fields(n_scalar_fields+1)%ptr => data%w%x
       end if
       ! This is very stupid...
       ! Some compilers cannot handle that these pointers dont point to anything
       ! (although they are not used) this fixes a segfault due to this.
       if (nelv .eq. 0) then
          allocate(tempo(1))
          x%ptr => tempo
          y%ptr => tempo
          z%ptr => tempo
          u%ptr => tempo
          v%ptr => tempo
          w%ptr => tempo
          p%ptr => tempo
          tem%ptr => tempo
       end if
 
       allocate(idx(nelv))
       do i = 1, nelv
          idx(i) = data%idx(i)
       end do
    type is (field_t)
       call dummy_list%init(1)
       ptr => data ! For the sake of intel
       call dummy_list%assign(1, ptr)
       call fld_file_select_from_field_list(this, dummy_list, p, u, v, w, tem, &
            scalar_fields, n_scalar_fields, write_pressure, write_velocity, &
            write_temperature, dof)
       nullify(dummy_list%items(1)%ptr)
       nullify(ptr)
       deallocate(dummy_list%items)
    type is (field_list_t)
       call fld_file_select_from_field_list(this, data, p, u, v, w, tem, &
            scalar_fields, n_scalar_fields, write_pressure, write_velocity, &
            write_temperature, dof)
    class default
       call neko_error('Invalid data')
    end select
    ! Fix things for pointers that do not exist in all data types...
    if (associated(dof)) then
       x%ptr => dof%x(:,1,1,1)
       y%ptr => dof%y(:,1,1,1)
       z%ptr => dof%z(:,1,1,1)
       msh => dof%msh
       xh => dof%Xh
    end if
 
    if (associated(msh)) then
       nelv = msh%nelv
       glb_nelv = msh%glb_nelv
       offset_el = msh%offset_el
       gdim = msh%gdim
       ! Store global idx of each element
       allocate(idx(msh%nelv))
       do i = 1, msh%nelv
          idx(i) = msh%elements(i)%e%id()
       end do
    end if
 
    if (associated(xh)) then
       lx = xh%lx
       ly = xh%ly
       lz = xh%lz
    end if
 
    lxyz = lx*ly*lz
    n = nelv*lxyz
 
    if (this%dp_precision) then
       fld_data_size = mpi_double_precision_size
    else
       fld_data_size = mpi_real_size
    end if
    if (this%dp_precision) then
       allocate(tmp_dp(gdim*n))
    else
       allocate(tmp_sp(gdim*n))
    end if
 
 
    !
    ! Create fld header for NEKTON's multifile output
    !
 
    call this%increment_counter()
    ! Check if I should write the mesh. Always override at the start counters
    if (.not. this%write_mesh) then
       write_mesh = (this%get_counter() .eq. this%get_start_counter())
    else
       write_mesh = this%write_mesh
    end if
    call mpi_allreduce(mpi_in_place, write_mesh, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, write_velocity, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, write_pressure, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, write_temperature, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, n_scalar_fields, 1, &
         mpi_integer, mpi_max, neko_comm)
 
    ! Build rdcode note that for field_t, we only support scalar
    ! fields at the moment
    rdcode = ' '
    i = 1
    if (write_mesh) then
       rdcode(i) = 'X'
       i = i + 1
    end if
    if (write_velocity) then
       rdcode(i) = 'U'
       i = i + 1
    end if
    if (write_pressure) then
       rdcode(i) = 'P'
       i = i + 1
    end if
    if (write_temperature) then
       rdcode(i) = 'T'
       i = i + 1
    end if
    if (n_scalar_fields .gt. 0 ) then
       rdcode(i) = 'S'
       i = i + 1
       write(rdcode(i), '(i1)') (n_scalar_fields)/10
       i = i + 1
       write(rdcode(i), '(i1)') (n_scalar_fields) - 10*((n_scalar_fields)/10)
       i = i + 1
    end if
 
    write(hdr, 1) fld_data_size, lx, ly, lz, glb_nelv, glb_nelv,&
         time, this%get_counter(), 1, 1, (rdcode(i), i = 1, 10)
1   format('#std', 1x, i1, 1x, i2, 1x, i2, 1x, i2, 1x, i10, 1x, i10, &
         1x, e20.13, 1x, i9, 1x, i6, 1x, i6, 1x, 10a)
 
    ! Change to NEKTON's fld file format
    fname = this%get_fld_fname()
 
    call mpi_file_open(neko_comm, trim(fname), &
         mpi_mode_wronly + mpi_mode_create, mpi_info_null, fh, &
         ierr)
 
    call mpi_file_write_all(fh, hdr, 132, mpi_character, status, ierr)
    mpi_offset = 132 * mpi_character_size
 
    call mpi_file_write_all(fh, test_pattern, 1, mpi_real, status, ierr)
    mpi_offset = mpi_offset + mpi_real_size
 
    byte_offset = mpi_offset + &
         int(offset_el, i8) * int(mpi_integer_size, i8)
    call mpi_file_write_at_all(fh, byte_offset, idx, nelv, &
         mpi_integer, status, ierr)
    mpi_offset = mpi_offset + int(glb_nelv, i8) * int(mpi_integer_size, i8)
    deallocate(idx)
    if (write_mesh) then
 
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int(gdim*lxyz, i8) * &
            int(fld_data_size, i8))
       call fld_file_write_vector_field(this, fh, byte_offset, &
            x%ptr, y%ptr, z%ptr, &
            n, gdim, lxyz, nelv)
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int(gdim *lxyz, i8) * &
            int(fld_data_size, i8))
    end if
    if (write_velocity) then
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int(gdim * (lxyz), i8) * int(fld_data_size, i8))
       call fld_file_write_vector_field(this, fh, byte_offset, &
            u%ptr, v%ptr, w%ptr, n, gdim, lxyz, nelv)
 
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int(gdim * (lxyz), i8) * &
            int(fld_data_size, i8))
 
    end if
 
    if (write_pressure) then
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int((lxyz), i8) * int(fld_data_size, i8))
       call fld_file_write_field(this, fh, byte_offset, p%ptr, n)
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int((lxyz), i8) * int(fld_data_size, i8))
    end if
 
    if (write_temperature) then
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int((lxyz), i8) * &
            int(fld_data_size, i8))
       call fld_file_write_field(this, fh, byte_offset, tem%ptr, n)
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int((lxyz), i8) * &
            int(fld_data_size, i8))
    end if
 
    temp_offset = mpi_offset
 
    do i = 1, n_scalar_fields
       ! Without this redundant if statement Cray optimizes this loop to
       ! Oblivion
       if (i .eq. 2) then
          mpi_offset = int(temp_offset, i8) + int(1_i8*glb_nelv, i8) * &
               (int(lxyz, i8) * int(fld_data_size, i8))
       end if
       byte_offset = int(mpi_offset, i8) + int(offset_el, i8) * &
            (int((lxyz), i8) * &
            int(fld_data_size, i8))
       call fld_file_write_field(this, fh, byte_offset, scalar_fields(i)%ptr, n)
       mpi_offset = int(mpi_offset, i8) + int(glb_nelv, i8) * &
            (int(lxyz, i8) * &
            int(fld_data_size, i8))
    end do
 
    if (gdim .eq. 3) then
 
       if (write_mesh) then
          !The offset is:
          ! mpioff + element_off * 2(min max value)
          !           * 4(single precision) * gdim(dimensions)
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
          call fld_file_write_metadata_vector(this, fh, byte_offset, &
               x%ptr, y%ptr, z%ptr, gdim, lxyz, nelv)
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
       end if
 
       if (write_velocity) then
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
          call fld_file_write_metadata_vector(this, fh, byte_offset, &
               u%ptr, v%ptr, w%ptr, gdim, lxyz, nelv)
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
 
       end if
 
       if (write_pressure) then
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
          call fld_file_write_metadata_scalar(this, fh, byte_offset, &
               p%ptr, lxyz, nelv)
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
 
       end if
 
       if (write_temperature) then
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
          call fld_file_write_metadata_scalar(this, fh, byte_offset, &
               tem%ptr, lxyz, nelv)
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
 
       end if
 
 
 
       temp_offset = mpi_offset
 
       do i = 1, n_scalar_fields
          ! Without this redundant if statement, Cray optimizes this loop to
          ! Oblivion
          if (i .eq. 2) then
             mpi_offset = int(temp_offset, i8) + &
                  int(1_i8*glb_nelv, i8) * &
                  int(2, i8) * &
                  int(mpi_real_size, i8)
          end if
 
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
          call fld_file_write_metadata_scalar(this, fh, byte_offset, &
               scalar_fields(i)%ptr, lxyz, nelv)
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
       end do
    end if
 
 
    call mpi_file_sync(fh, ierr)
    call mpi_file_close(fh, ierr)
    ! Write metadata file
    if (pe_rank .eq. 0) then
       call filename_name(this%get_base_fname(), name)
 
       open(newunit = file_unit, &
            file = this%get_meta_fname(), status = 'replace')
       ! The following string will specify that the files in the file series
       ! are defined by the filename followed by a 0.
       ! This 0 is necessary as it specifies the index of number of files
       ! the output file is split across.
       ! In the past, many .f files were generated for each write.
       ! To be consistent with this the trailing 0 is still necessary today.
       write(file_unit, fmt = '(A,A,A)') 'filetemplate:         ', &
            trim(name), '%01d.f%05d'
       write(file_unit, fmt = '(A,i5)') 'firsttimestep: ', &
            this%get_start_counter()
       write(file_unit, fmt = '(A,i5)') 'numtimesteps: ', &
            (this%get_counter() + 1) - this%get_start_counter()
       close(file_unit)
    end if
 
    if (allocated(tmp_dp)) deallocate(tmp_dp)
    if (allocated(tmp_sp)) deallocate(tmp_sp)
    if (allocated(tempo)) deallocate(tempo)
    if (allocated(scalar_fields)) deallocate(scalar_fields)
 
  end subroutine fld_file_write
 
  subroutine fld_file_write_masked(this, data, mask, t)
    class(fld_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    type(mask_t), intent(in) :: mask
    real(kind=rp), intent(in), optional :: t
    type(array_ptr_t) :: x, y, z, u, v, w, p, tem
    real(kind=rp), allocatable, target :: tempo(:)
    type(mesh_t), pointer :: msh
    type(dofmap_t), pointer :: dof
    type(space_t), pointer :: Xh
    real(kind=dp) :: time
    character(len= 132) :: hdr
    character :: rdcode(10)
    character(len=6) :: id_str
    character(len= 1024) :: fname
    character(len= 1024) :: name
    integer :: file_unit
    integer :: i, ierr, n, suffix_pos, tslash_pos
    integer :: lx, ly, lz, lxyz, gdim, glb_nelv, nelv, offset_el
    integer, allocatable :: idx(:)
    type(mpi_status) :: status
    type(mpi_file) :: fh
    integer (kind=MPI_OFFSET_KIND) :: mpi_offset, byte_offset, temp_offset
    real(kind=sp), parameter :: test_pattern = 6.54321
    type(array_ptr_t), allocatable :: scalar_fields(:)
    logical :: write_mesh, write_velocity, write_pressure, write_temperature
    integer :: fld_data_size, n_scalar_fields
    type(field_list_t) :: dummy_list
    type(field_t), pointer :: ptr
 
    if (present(t)) then
       time = real(t, dp)
    else
       time = 0d0
    end if
 
    nullify(msh)
    nullify(dof)
    nullify(xh)
    n_scalar_fields = 0
    write_pressure = .false.
    write_velocity = .false.
    write_temperature = .false.
 
    select type (data)
    type is (fld_file_data_t)
       nelv = data%nelv
       lx = data%lx
       ly = data%ly
       lz = data%lz
       gdim = data%gdim
       glb_nelv = data%glb_nelv
       offset_el = data%offset_el
 
       if (data%x%size() .gt. 0) x%ptr => data%x%x
       if (data%y%size() .gt. 0) y%ptr => data%y%x
       if (data%z%size() .gt. 0) z%ptr => data%z%x
       if (gdim .eq. 2) z%ptr => data%y%x
       if (data%u%size() .gt. 0) then
          u%ptr => data%u%x
          ! In case only u is actually allocated, point the other comps to u
          ! so that we don't die on trying to write them
          if (data%v%size() .le. 0) v%ptr => data%u%x
          if (data%w%size() .le. 0) w%ptr => data%u%x
          write_velocity = .true.
       end if
       if (data%v%size() .gt. 0) v%ptr => data%v%x
       if (data%w%size() .gt. 0) w%ptr => data%w%x
       if (data%p%size() .gt. 0) then
          p%ptr => data%p%x
          write_pressure = .true.
       end if
       if (data%t%size() .gt. 0) then
          write_temperature = .true.
          tem%ptr => data%t%x
       end if
       ! If gdim = 2 and Z-velocity component exists,
       ! it is stored in last scalar field
       if (gdim .eq. 2 .and. data%w%size() .gt. 0) then
          n_scalar_fields = data%n_scalars + 1
          allocate(scalar_fields(n_scalar_fields))
          do i = 1, n_scalar_fields -1
             scalar_fields(i)%ptr => data%s(i)%x
          end do
          scalar_fields(n_scalar_fields)%ptr => data%w%x
       else
          n_scalar_fields = data%n_scalars
          allocate(scalar_fields(n_scalar_fields+1))
          do i = 1, n_scalar_fields
             scalar_fields(i)%ptr => data%s(i)%x
          end do
          scalar_fields(n_scalar_fields+1)%ptr => data%w%x
       end if
       ! This is very stupid...
       ! Some compilers cannot handle that these pointers dont point to anything
       ! (although they are not used) this fixes a segfault due to this.
       if (nelv .eq. 0) then
          allocate(tempo(1))
          x%ptr => tempo
          y%ptr => tempo
          z%ptr => tempo
          u%ptr => tempo
          v%ptr => tempo
          w%ptr => tempo
          p%ptr => tempo
          tem%ptr => tempo
       end if
 
       allocate(idx(nelv))
       do i = 1, nelv
          idx(i) = data%idx(i)
       end do
    type is (field_t)
       call dummy_list%init(1)
       ptr => data ! For the sake of intel
       call dummy_list%assign(1, ptr)
       call fld_file_select_from_field_list(this, dummy_list, p, u, v, w, tem, &
            scalar_fields, n_scalar_fields, write_pressure, write_velocity, &
            write_temperature, dof)
       nullify(dummy_list%items(1)%ptr)
       nullify(ptr)
       deallocate(dummy_list%items)
    type is (field_list_t)
       call fld_file_select_from_field_list(this, data, p, u, v, w, tem, &
            scalar_fields, n_scalar_fields, write_pressure, write_velocity, &
            write_temperature, dof)
    class default
       call neko_error('Invalid data')
    end select
    ! Fix things for pointers that do not exist in all data types...
    if (associated(dof)) then
       x%ptr => dof%x(:,1,1,1)
       y%ptr => dof%y(:,1,1,1)
       z%ptr => dof%z(:,1,1,1)
       msh => dof%msh
       xh => dof%Xh
    end if
 
    if (associated(msh)) then
       nelv = msh%nelv
       glb_nelv = msh%glb_nelv
       offset_el = msh%offset_el
       gdim = msh%gdim
       ! Store global idx of each element
       allocate(idx(msh%nelv))
       do i = 1, msh%nelv
          idx(i) = msh%elements(i)%e%id()
       end do
    end if
 
    if (associated(xh)) then
       lx = xh%lx
       ly = xh%ly
       lz = xh%lz
    end if
 
    ! Up to now, all data types have dealt with their stuff.
    ! Now overwrite with the masked info
    lxyz = lx*ly*lz
    nelv = mask%size() / lxyz
    if (mod(mask%size(), lxyz) /= 0) then
       call neko_error("Mask size must be a multiple of the number of elements in the mesh.")
    end if
    call mpi_allreduce(nelv, glb_nelv, 1, &
         mpi_integer, mpi_sum, neko_comm)
    call mpi_scan(nelv, offset_el, 1, &
         mpi_integer, mpi_sum, neko_comm, ierr)
    offset_el = offset_el - nelv
 
    if (allocated(idx)) then
       deallocate(idx)
    end if
 
    allocate(idx(nelv))
    do i = 1, nelv
       idx(i) = offset_el + i
    end do
    n = nelv*lxyz
 
    if (this%dp_precision) then
       fld_data_size = mpi_double_precision_size
    else
       fld_data_size = mpi_real_size
    end if
    if (this%dp_precision) then
       allocate(tmp_dp(gdim*n))
    else
       allocate(tmp_sp(gdim*n))
    end if
 
 
    !
    ! Create fld header for NEKTON's multifile output
    !
 
    call this%increment_counter()
    ! Check if I should write the mesh. Always override at the start counters
    if (.not. this%write_mesh) then
       write_mesh = (this%get_counter() .eq. this%get_start_counter())
    else
       write_mesh = this%write_mesh
    end if
    call mpi_allreduce(mpi_in_place, write_mesh, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, write_velocity, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, write_pressure, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, write_temperature, 1, &
         mpi_logical, mpi_lor, neko_comm)
    call mpi_allreduce(mpi_in_place, n_scalar_fields, 1, &
         mpi_integer, mpi_max, neko_comm)
 
    ! Build rdcode note that for field_t, we only support scalar
    ! fields at the moment
    rdcode = ' '
    i = 1
    if (write_mesh) then
       rdcode(i) = 'X'
       i = i + 1
    end if
    if (write_velocity) then
       rdcode(i) = 'U'
       i = i + 1
    end if
    if (write_pressure) then
       rdcode(i) = 'P'
       i = i + 1
    end if
    if (write_temperature) then
       rdcode(i) = 'T'
       i = i + 1
    end if
    if (n_scalar_fields .gt. 0 ) then
       rdcode(i) = 'S'
       i = i + 1
       write(rdcode(i), '(i1)') (n_scalar_fields)/10
       i = i + 1
       write(rdcode(i), '(i1)') (n_scalar_fields) - 10*((n_scalar_fields)/10)
       i = i + 1
    end if
 
    write(hdr, 1) fld_data_size, lx, ly, lz, glb_nelv, glb_nelv,&
         time, this%get_counter(), 1, 1, (rdcode(i), i = 1, 10)
1   format('#std', 1x, i1, 1x, i2, 1x, i2, 1x, i2, 1x, i10, 1x, i10, &
         1x, e20.13, 1x, i9, 1x, i6, 1x, i6, 1x, 10a)
 
    ! Change to NEKTON's fld file format
    fname = this%get_fld_fname()
 
    call mpi_file_open(neko_comm, trim(fname), &
         mpi_mode_wronly + mpi_mode_create, mpi_info_null, fh, &
         ierr)
 
    call mpi_file_write_all(fh, hdr, 132, mpi_character, status, ierr)
    mpi_offset = 132 * mpi_character_size
 
    call mpi_file_write_all(fh, test_pattern, 1, mpi_real, status, ierr)
    mpi_offset = mpi_offset + mpi_real_size
 
    byte_offset = mpi_offset + &
         int(offset_el, i8) * int(mpi_integer_size, i8)
    call mpi_file_write_at_all(fh, byte_offset, idx, nelv, &
         mpi_integer, status, ierr)
    mpi_offset = mpi_offset + int(glb_nelv, i8) * int(mpi_integer_size, i8)
    deallocate(idx)
    if (write_mesh) then
 
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int(gdim*lxyz, i8) * &
            int(fld_data_size, i8))
       call fld_file_write_vector_field_masked(this, fh, byte_offset, &
            x%ptr, y%ptr, z%ptr, &
            n, gdim, lxyz, nelv, lx, ly, lz, mask%get())
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int(gdim *lxyz, i8) * &
            int(fld_data_size, i8))
    end if
    if (write_velocity) then
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int(gdim * (lxyz), i8) * int(fld_data_size, i8))
       call fld_file_write_vector_field_masked(this, fh, byte_offset, &
            u%ptr, v%ptr, w%ptr, n, gdim, lxyz, nelv, lx, ly, lz, mask%get())
 
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int(gdim * (lxyz), i8) * &
            int(fld_data_size, i8))
 
    end if
 
    if (write_pressure) then
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int((lxyz), i8) * int(fld_data_size, i8))
       call fld_file_write_field_masked(this, fh, byte_offset, p%ptr, n, mask%get())
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int((lxyz), i8) * int(fld_data_size, i8))
    end if
 
    if (write_temperature) then
       byte_offset = mpi_offset + int(offset_el, i8) * &
            (int((lxyz), i8) * &
            int(fld_data_size, i8))
       call fld_file_write_field_masked(this, fh, byte_offset, tem%ptr, n, mask%get())
       mpi_offset = mpi_offset + int(glb_nelv, i8) * &
            (int((lxyz), i8) * &
            int(fld_data_size, i8))
    end if
 
    temp_offset = mpi_offset
 
    do i = 1, n_scalar_fields
       ! Without this redundant if statement Cray optimizes this loop to
       ! Oblivion
       if (i .eq. 2) then
          mpi_offset = int(temp_offset, i8) + int(1_i8*glb_nelv, i8) * &
               (int(lxyz, i8) * int(fld_data_size, i8))
       end if
       byte_offset = int(mpi_offset, i8) + int(offset_el, i8) * &
            (int((lxyz), i8) * &
            int(fld_data_size, i8))
       call fld_file_write_field_masked(this, fh, byte_offset, scalar_fields(i)%ptr, n, mask%get())
       mpi_offset = int(mpi_offset, i8) + int(glb_nelv, i8) * &
            (int(lxyz, i8) * &
            int(fld_data_size, i8))
    end do
 
    if (gdim .eq. 3) then
 
       if (write_mesh) then
          !The offset is:
          ! mpioff + element_off * 2(min max value)
          !           * 4(single precision) * gdim(dimensions)
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
          call fld_file_write_metadata_vector_masked(this, fh, byte_offset, &
               x%ptr, y%ptr, z%ptr, gdim, lxyz, nelv, lx, ly, lz, n, mask%get())
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
       end if
 
       if (write_velocity) then
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
          call fld_file_write_metadata_vector_masked(this, fh, byte_offset, &
               u%ptr, v%ptr, w%ptr, gdim, lxyz, nelv, lx, ly, lz, n, mask%get())
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8) * &
               int(gdim, i8)
 
       end if
 
       if (write_pressure) then
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
          call fld_file_write_metadata_scalar_masked(this, fh, byte_offset, &
               p%ptr, lxyz, nelv, lx, ly, lz, n, mask%get())
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
 
       end if
 
       if (write_temperature) then
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
          call fld_file_write_metadata_scalar_masked(this, fh, byte_offset, &
               tem%ptr, lxyz, nelv, lx, ly, lz, n, mask%get())
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
 
       end if
 
 
 
       temp_offset = mpi_offset
 
       do i = 1, n_scalar_fields
          ! Without this redundant if statement, Cray optimizes this loop to
          ! Oblivion
          if (i .eq. 2) then
             mpi_offset = int(temp_offset, i8) + &
                  int(1_i8*glb_nelv, i8) * &
                  int(2, i8) * &
                  int(mpi_real_size, i8)
          end if
 
          byte_offset = int(mpi_offset, i8) + &
               int(offset_el, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
          call fld_file_write_metadata_scalar_masked(this, fh, byte_offset, &
               scalar_fields(i)%ptr, lxyz, nelv, lx, ly, lz, n, mask%get())
          mpi_offset = int(mpi_offset, i8) + &
               int(glb_nelv, i8) * &
               int(2, i8) * &
               int(mpi_real_size, i8)
       end do
    end if
 
 
    call mpi_file_sync(fh, ierr)
    call mpi_file_close(fh, ierr)
    ! Write metadata file
    if (pe_rank .eq. 0) then
       call filename_name(this%get_base_fname(), name)
 
       open(newunit = file_unit, &
            file = this%get_meta_fname(), status = 'replace')
       ! The following string will specify that the files in the file series
       ! are defined by the filename followed by a 0.
       ! This 0 is necessary as it specifies the index of number of files
       ! the output file is split across.
       ! In the past, many .f files were generated for each write.
       ! To be consistent with this the trailing 0 is still necessary today.
       write(file_unit, fmt = '(A,A,A)') 'filetemplate:         ', &
            trim(name), '%01d.f%05d'
       write(file_unit, fmt = '(A,i5)') 'firsttimestep: ', &
            this%get_start_counter()
       write(file_unit, fmt = '(A,i5)') 'numtimesteps: ', &
            (this%get_counter() + 1) - this%get_start_counter()
       close(file_unit)
    end if
 
    if (allocated(tmp_dp)) deallocate(tmp_dp)
    if (allocated(tmp_sp)) deallocate(tmp_sp)
    if (allocated(tempo)) deallocate(tempo)
    if (allocated(scalar_fields)) deallocate(scalar_fields)
 
  end subroutine fld_file_write_masked
 
  subroutine fld_file_write_metadata_vector(this, fh, byte_offset, x, y, z, &
       gdim, lxyz, nelv)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(in) :: gdim, lxyz, nelv
    real(kind=rp), intent(in) :: x(lxyz, nelv), y(lxyz, nelv), z(lxyz, nelv)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer :: el, j, ierr, nout
    type(mpi_status) :: status
    real(kind=sp) :: buffer(2*gdim*nelv)
 
    j = 1
    do el = 1, nelv
       buffer(j+0) = real(vlmin(x(1, el), lxyz), sp)
       buffer(j+1) = real(vlmax(x(1, el), lxyz), sp)
       buffer(j+2) = real(vlmin(y(1, el), lxyz), sp)
       buffer(j+3) = real(vlmax(y(1, el), lxyz), sp)
       j = j + 4
       if (gdim .eq. 3) then
          buffer(j+0) = real(vlmin(z(1, el), lxyz), sp)
          buffer(j+1) = real(vlmax(z(1, el), lxyz), sp)
          j = j + 2
       end if
    end do
 
    ! write out data
    nout = 2*gdim*nelv
 
    call mpi_file_write_at_all(fh, byte_offset, buffer, nout, &
         mpi_real, status, ierr)
 
  end subroutine fld_file_write_metadata_vector
 
  subroutine fld_file_write_metadata_vector_masked(this, fh, byte_offset, x, y, z, &
       gdim, lxyz, nelv, lx, ly, lz, n, mask)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(in) :: gdim, lxyz, nelv, lx, ly, lz, n
    real(kind=rp), intent(in) :: x(lxyz, *), y(lxyz, *), z(lxyz, *)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer, intent(in) :: mask(n)
    integer :: el, j, ierr, nout, i_m, nidx(4), e_m
    type(mpi_status) :: status
    real(kind=sp) :: buffer(2*gdim*nelv)
 
    j = 1
    do el = 1, nelv
       i_m = 1 + lxyz * (el - 1) ! Offset in the mask array
       nidx = nonlinear_index(mask(i_m), lx, ly, lz)
       e_m = nidx(4) ! Actual element in the field array
       buffer(j+0) = real(vlmin(x(1, e_m), lxyz), sp)
       buffer(j+1) = real(vlmax(x(1, e_m), lxyz), sp)
       buffer(j+2) = real(vlmin(y(1, e_m), lxyz), sp)
       buffer(j+3) = real(vlmax(y(1, e_m), lxyz), sp)
       j = j + 4
       if (gdim .eq. 3) then
          buffer(j+0) = real(vlmin(z(1, e_m), lxyz), sp)
          buffer(j+1) = real(vlmax(z(1, e_m), lxyz), sp)
          j = j + 2
       end if
    end do
 
    ! write out data
    nout = 2*gdim*nelv
 
    call mpi_file_write_at_all(fh, byte_offset, buffer, nout, &
         mpi_real, status, ierr)
 
  end subroutine fld_file_write_metadata_vector_masked
 
  subroutine fld_file_write_metadata_scalar(this, fh, byte_offset, x, lxyz, &
       nelv)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(in) :: lxyz, nelv
    real(kind=rp), intent(in) :: x(lxyz, nelv)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer :: el, j, ierr, nout
    type(mpi_status) :: status
    real(kind=sp) :: buffer(2*nelv)
 
    j = 1
    do el = 1, nelv
       buffer(j+0) = real(vlmin(x(1, el), lxyz), sp)
       buffer(j+1) = real(vlmax(x(1, el), lxyz), sp)
       j = j + 2
    end do
 
    ! write out data
    nout = 2*nelv
 
    call mpi_file_write_at_all(fh, byte_offset, buffer, nout, &
         mpi_real, status, ierr)
 
  end subroutine fld_file_write_metadata_scalar
 
  subroutine fld_file_write_metadata_scalar_masked(this, fh, byte_offset, x, lxyz, &
       nelv, lx, ly, lz, n, mask)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(in) :: lxyz, nelv, lx, ly, lz, n
    real(kind=rp), intent(in) :: x(lxyz, *)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer, intent(in) :: mask(n)
    integer :: el, j, ierr, nout, i_m, nidx(4), e_m
    type(mpi_status) :: status
    real(kind=sp) :: buffer(2*nelv)
 
    j = 1
    do el = 1, nelv
       i_m = 1 + lxyz * (el - 1) ! Offset in the mask array
       nidx = nonlinear_index(mask(i_m), lx, ly, lz)
       e_m = nidx(4) ! Actual element in the field array
       buffer(j+0) = real(vlmin(x(1, e_m), lxyz), sp)
       buffer(j+1) = real(vlmax(x(1, e_m), lxyz), sp)
       j = j + 2
    end do
 
    ! write out data
    nout = 2*nelv
 
    call mpi_file_write_at_all(fh, byte_offset, buffer, nout, &
         mpi_real, status, ierr)
 
  end subroutine fld_file_write_metadata_scalar_masked
 
  subroutine fld_file_write_field(this, fh, byte_offset, p, n)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(inout) :: n
    real(kind=rp), intent(inout) :: p(n)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer :: i, ierr
    type(mpi_status) :: status
 
    if ( this%dp_precision) then
       do i = 1, n
          tmp_dp(i) = real(p(i), dp)
       end do
 
       call mpi_file_write_at_all(fh, byte_offset, tmp_dp, n, &
            mpi_double_precision, status, ierr)
    else
       do i = 1, n
          tmp_sp(i) = real(p(i), sp)
       end do
       call mpi_file_write_at_all(fh, byte_offset, tmp_sp, n, &
            mpi_real, status, ierr)
    end if
 
  end subroutine fld_file_write_field
 
  subroutine fld_file_write_field_masked(this, fh, byte_offset, p, n, mask)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(inout) :: n
    real(kind=rp), intent(inout) :: p(:)
    integer, intent(in) :: mask(n)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer :: i, ierr
    type(mpi_status) :: status
 
    if ( this%dp_precision) then
       do i = 1, n
          tmp_dp(i) = real(p(mask(i)), dp)
       end do
 
       call mpi_file_write_at_all(fh, byte_offset, tmp_dp, n, &
            mpi_double_precision, status, ierr)
    else
       do i = 1, n
          tmp_sp(i) = real(p(mask(i)), sp)
       end do
       call mpi_file_write_at_all(fh, byte_offset, tmp_sp, n, &
            mpi_real, status, ierr)
    end if
 
  end subroutine fld_file_write_field_masked
 
  subroutine fld_file_write_vector_field(this, fh, byte_offset, x, y, z, n, &
       gdim, lxyz, nelv)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(in) :: n, gdim, lxyz, nelv
    real(kind=rp), intent(in) :: x(lxyz, nelv), y(lxyz, nelv), z(lxyz, nelv)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer :: i, el, j, ierr
    type(mpi_status) :: status
 
    if (this%dp_precision) then
       i = 1
       do el = 1, nelv
          do j = 1, lxyz
             tmp_dp(i) = real(x(j, el), dp)
             i = i +1
          end do
          do j = 1, lxyz
             tmp_dp(i) = real(y(j, el), dp)
             i = i +1
          end do
          if (gdim .eq. 3) then
             do j = 1, lxyz
                tmp_dp(i) = real(z(j, el), dp)
                i = i +1
             end do
          end if
       end do
       call mpi_file_write_at_all(fh, byte_offset, tmp_dp, gdim*n, &
            mpi_double_precision, status, ierr)
    else
       i = 1
       do el = 1, nelv
          do j = 1, lxyz
             tmp_sp(i) = real(x(j, el), sp)
             i = i +1
          end do
          do j = 1, lxyz
             tmp_sp(i) = real(y(j, el), sp)
             i = i +1
          end do
          if (gdim .eq. 3) then
             do j = 1, lxyz
                tmp_sp(i) = real(z(j, el), sp)
                i = i +1
             end do
          end if
       end do
       call mpi_file_write_at_all(fh, byte_offset, tmp_sp, gdim*n, &
            mpi_real, status, ierr)
    end if
 
 
  end subroutine fld_file_write_vector_field
 
  subroutine fld_file_write_vector_field_masked(this, fh, byte_offset, x, y, z, n, &
       gdim, lxyz, nelv, lx, ly, lz, mask)
    class(fld_file_t), intent(inout) :: this
    type(mpi_file), intent(inout) :: fh
    integer, intent(in) :: n, gdim, lxyz, lx, ly, lz, nelv
    real(kind=rp), intent(in) :: x(lxyz, *), y(lxyz, *), z(lxyz, *)
    integer (kind=MPI_OFFSET_KIND), intent(in) :: byte_offset
    integer, intent(in) :: mask(n)
    integer :: i, el, j, ierr, i_m, e_m, nidx(4)
    type(mpi_status) :: status
 
    if (this%dp_precision) then
       i = 1
       do el = 1, nelv
          i_m = 1 + lxyz * (el - 1) ! Offset in the mask array
          nidx = nonlinear_index(mask(i_m), lx, ly, lz)
          e_m = nidx(4) ! Actual element in the field array
          do j = 1, lxyz
             tmp_dp(i) = real(x(j, e_m), dp)
             i = i +1
          end do
          do j = 1, lxyz
             tmp_dp(i) = real(y(j, e_m), dp)
             i = i +1
          end do
          if (gdim .eq. 3) then
             do j = 1, lxyz
                tmp_dp(i) = real(z(j, e_m), dp)
                i = i +1
             end do
          end if
       end do
       call mpi_file_write_at_all(fh, byte_offset, tmp_dp, gdim*n, &
            mpi_double_precision, status, ierr)
    else
       i = 1
       do el = 1, nelv
          i_m = 1 + lxyz * (el - 1) ! Offset in the mask array
          nidx = nonlinear_index(mask(i_m), lx, ly, lz)
          e_m = nidx(4) ! Actual element in the field array
          do j = 1, lxyz
             tmp_sp(i) = real(x(j, e_m), sp)
             i = i +1
          end do
          do j = 1, lxyz
             tmp_sp(i) = real(y(j, e_m), sp)
             i = i +1
          end do
          if (gdim .eq. 3) then
             do j = 1, lxyz
                tmp_sp(i) = real(z(j, e_m), sp)
                i = i +1
             end do
          end if
       end do
       call mpi_file_write_at_all(fh, byte_offset, tmp_sp, gdim*n, &
            mpi_real, status, ierr)
    end if
 
 
  end subroutine fld_file_write_vector_field_masked
 
  subroutine fld_file_read(this, data)
    class(fld_file_t) :: this
    class(*), target, intent(inout) :: data
    character(len= 132) :: hdr
    integer :: ierr, suffix_pos, i, j
    type(mpi_file) :: fh
    type(mpi_status) :: status
    character(len= 1024) :: fname, base_fname, meta_fname, string, path
    logical :: meta_file, read_mesh, read_velocity, read_pressure
    logical :: read_temp
    character(len=6) :: suffix
    integer (kind=MPI_OFFSET_KIND) :: mpi_offset, byte_offset
    integer :: lx, ly, lz, glb_nelv, counter, lxyz
    integer :: FLD_DATA_SIZE, n_scalars, n
    integer :: file_unit
    real(kind=rp) :: time
    real(kind=sp) :: temp
    type(linear_dist_t) :: dist
    real(kind=sp), parameter :: test_pattern = 6.54321
    character :: rdcode(10), temp_str(4)
    character(len=LOG_SIZE) :: log_buf
 
    select type (data)
    type is (fld_file_data_t)
       call filename_chsuffix(this%get_base_fname(), meta_fname, 'nek5000')
 
       inquire(file = trim(meta_fname), exist = meta_file)
       if (meta_file .and. data%meta_nsamples .eq. 0) then
          if (pe_rank .eq. 0) then
             open(newunit = file_unit, file = trim(meta_fname))
             read(file_unit, fmt = '(A)') string
             read(string(14:), fmt = '(A)') string
             string = trim(string)
 
             data%fld_series_fname = string(:scan(trim(string), '%')-1)
             data%fld_series_fname = adjustl(data%fld_series_fname)
             data%fld_series_fname = trim(data%fld_series_fname)//'0'
 
             read(file_unit, fmt = '(A)') string
             read(string(scan(string, ':')+1:), *) data%meta_start_counter
             read(file_unit, fmt = '(A)') string
             read(string(scan(string, ':')+1:), *) data%meta_nsamples
             close(file_unit)
 
             write(log_buf,*) 'Reading meta file for fld series'
             call neko_log%message(log_buf)
             write(log_buf,*) 'Name: ', trim(data%fld_series_fname)
             call neko_log%message(log_buf)
             write(log_buf,*) 'Start counter: ', data%meta_start_counter
             call neko_log%message(log_buf)
             write(log_buf,*) 'Nsamples: ', data%meta_nsamples
             call neko_log%message(log_buf)
 
          end if
          call mpi_bcast(data%fld_series_fname, 1024, mpi_character, 0, &
               neko_comm, ierr)
          call mpi_bcast(data%meta_start_counter, 1, mpi_integer, 0, &
               neko_comm, ierr)
          call mpi_bcast(data%meta_nsamples, 1, mpi_integer, 0, &
               neko_comm, ierr)
 
          if (this%get_counter() .eq. -1) then
             call this%set_start_counter(data%meta_start_counter)
             call this%set_counter(data%meta_start_counter)
          end if
       end if
 
       if (meta_file) then
          call filename_path(this%get_base_fname(), path)
          write(suffix, '(a,i5.5)') 'f', this%get_counter()
          fname = trim(path) // trim(data%fld_series_fname) // '.' // suffix
          if (this%get_counter() .ge. &
               data%meta_nsamples+data%meta_start_counter) then
             call neko_error('Trying to read more fld files than exist')
          end if
       else
          write(suffix, '(a,i5.5)') 'f', this%get_counter()
          call filename_chsuffix(trim(this%get_base_fname()), fname, suffix)
       end if
       call mpi_file_open(neko_comm, trim(fname), &
            mpi_mode_rdonly, mpi_info_null, fh, ierr)
 
       if (ierr .ne. 0) call neko_error("Could not read "//trim(fname))
 
       call neko_log%message('Reading fld file ' // trim(fname))
 
       call mpi_file_read_all(fh, hdr, 132, mpi_character, status, ierr)
       ! This read can prorbably be done wihtout the temp variables,
       !      temp_str, i, j
 
       read(hdr, 1) temp_str, fld_data_size, lx, ly, lz, glb_nelv, glb_nelv, &
            time, counter, i, j, (rdcode(i), i = 1, 10)
1      format(4a, 1x, i1, 1x, i2, 1x, i2, 1x, i2, 1x, i10, 1x, i10, &
            1x, e20.13, 1x, i9, 1x, i6, 1x, i6, 1x, 10a)
       if (data%nelv .eq. 0) then
          dist = linear_dist_t(glb_nelv, pe_rank, pe_size, neko_comm)
          data%nelv = dist%num_local()
          data%offset_el = dist%start_idx()
       end if
       data%lx = lx
       data%ly = ly
       data%lz = lz
       data%glb_nelv = glb_nelv
       data%t_counter = counter
       data%time = time
       lxyz = lx * ly * lz
       n = lxyz * data%nelv
 
       if (lz .eq. 1) then
          data%gdim = 2
       else
          data%gdim = 3
       end if
 
 
       if (fld_data_size .eq. mpi_double_precision_size) then
          this%dp_precision = .true.
       else
          this%dp_precision = .false.
       end if
       if (this%dp_precision) then
          allocate(tmp_dp(data%gdim*n))
       else
          allocate(tmp_sp(data%gdim*n))
       end if
 
 
       i = 1
       read_mesh = .false.
       read_velocity = .false.
       read_pressure = .false.
       read_temp = .false.
       if (rdcode(i) .eq. 'X') then
          read_mesh = .true.
          call data%x%init(n)
          call data%y%init(n)
          call data%z%init(n)
          i = i + 1
       end if
       if (rdcode(i) .eq. 'U') then
          read_velocity = .true.
          call data%u%init(n)
          call data%v%init(n)
          call data%w%init(n)
          i = i + 1
       end if
       if (rdcode(i) .eq. 'P') then
          read_pressure = .true.
          call data%p%init(n)
          i = i + 1
       end if
       if (rdcode(i) .eq. 'T') then
          read_temp = .true.
          call data%t%init(n)
          i = i + 1
       end if
       n_scalars = 0
       if (rdcode(i) .eq. 'S') then
          i = i + 1
          read(rdcode(i),*) n_scalars
          n_scalars = n_scalars*10
          i = i + 1
          read(rdcode(i),*) j
          n_scalars = n_scalars+j
          i = i + 1
          if (allocated(data%s)) then
             if (data%n_scalars .ne. n_scalars) then
                do j = 1, data%n_scalars
                   call data%s(j)%free()
                end do
                deallocate(data%s)
                data%n_scalars = n_scalars
                allocate(data%s(n_scalars))
                do j = 1, data%n_scalars
                   call data%s(j)%init(n)
                end do
             end if
          else
             data%n_scalars = n_scalars
             allocate(data%s(data%n_scalars))
             do j = 1, data%n_scalars
                call data%s(j)%init(n)
             end do
          end if
          i = i + 1
       end if
 
       mpi_offset = 132 * mpi_character_size
       call mpi_file_read_at_all(fh, mpi_offset, temp, 1, &
            mpi_real, status, ierr)
       if (.not. sabscmp(temp, test_pattern, epsilon(1.0_sp))) then
          call neko_error('Incorrect format for fld file, &
          &test pattern does not match.')
       end if
       mpi_offset = mpi_offset + mpi_real_size
 
 
       if (allocated(data%idx)) then
          if (size(data%idx) .ne. data%nelv) then
             deallocate(data%idx)
             allocate(data%idx(data%nelv))
          end if
       else
          allocate(data%idx(data%nelv))
       end if
 
       byte_offset = mpi_offset + &
            int(data%offset_el, i8) * int(mpi_integer_size, i8)
 
       call mpi_file_read_at_all(fh, byte_offset, data%idx, data%nelv, &
            mpi_integer, status, ierr)
 
       mpi_offset = mpi_offset + &
            int(data%glb_nelv, i8) * int(mpi_integer_size, i8)
 
       if (read_mesh) then
          byte_offset = mpi_offset + int(data%offset_el, i8) * &
               (int(data%gdim*lxyz, i8) * &
               int(fld_data_size, i8))
          call fld_file_read_vector_field(this, fh, byte_offset, &
               data%x, data%y, data%z, data)
          mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
               (int(data%gdim *lxyz, i8) * &
               int(fld_data_size, i8))
       end if
 
       if (read_velocity) then
          byte_offset = mpi_offset + int(data%offset_el, i8) * &
               (int(data%gdim*lxyz, i8) * &
               int(fld_data_size, i8))
          call fld_file_read_vector_field(this, fh, byte_offset, &
               data%u, data%v, data%w, data)
          mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
               (int(data%gdim *lxyz, i8) * &
               int(fld_data_size, i8))
       end if
 
       if (read_pressure) then
          byte_offset = mpi_offset + int(data%offset_el, i8) * &
               (int(lxyz, i8) * &
               int(fld_data_size, i8))
          call fld_file_read_field(this, fh, byte_offset, data%p, data)
          mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
               (int(lxyz, i8) * &
               int(fld_data_size, i8))
       end if
 
       if (read_temp) then
          byte_offset = mpi_offset + int(data%offset_el, i8) * &
               (int(lxyz, i8) * &
               int(fld_data_size, i8))
          call fld_file_read_field(this, fh, byte_offset, data%t, data)
          mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
               (int(lxyz, i8) * &
               int(fld_data_size, i8))
       end if
 
       do i = 1, n_scalars
          byte_offset = mpi_offset + int(data%offset_el, i8) * &
               (int(lxyz, i8) * &
               int(fld_data_size, i8))
          call fld_file_read_field(this, fh, byte_offset, data%s(i), data)
          mpi_offset = mpi_offset + int(data%glb_nelv, i8) * &
               (int(lxyz, i8) * &
               int(fld_data_size, i8))
       end do
 
       call device_sync()
       call this%increment_counter()
 
       if (allocated(tmp_dp)) deallocate(tmp_dp)
       if (allocated(tmp_sp)) deallocate(tmp_sp)
    class default
       call neko_error('Currently we only read into fld_file_data_t, &
       &please use that data structure instead.')
    end select
 
  end subroutine fld_file_read
 
  subroutine fld_file_read_field(this, fh, byte_offset, x, fld_data)
    class(fld_file_t), intent(inout) :: this
    type(vector_t), intent(inout) :: x
    type(fld_file_data_t) :: fld_data
    integer(kind=MPI_OFFSET_KIND) :: byte_offset
    type(mpi_file) :: fh
    type(mpi_status) :: status
    integer :: n, ierr, lxyz, i
 
    n = x%size()
    lxyz = fld_data%lx * fld_data%ly * fld_data%lz
 
    if (this%dp_precision) then
       call mpi_file_read_at_all(fh, byte_offset, tmp_dp, n, &
            mpi_double_precision, status, ierr)
    else
       call mpi_file_read_at_all(fh, byte_offset, tmp_sp, n, &
            mpi_real, status, ierr)
    end if
 
    if (this%dp_precision) then
       do i = 1, n
          x%x(i) = tmp_dp(i)
       end do
    else
       do i = 1, n
          x%x(i) = tmp_sp(i)
       end do
    end if
 
    call x%copy_from(host_to_device, sync = .false.)
 
  end subroutine fld_file_read_field
 
 
  subroutine fld_file_read_vector_field(this, fh, byte_offset, &
       x, y, z, fld_data)
    class(fld_file_t), intent(inout) :: this
    type(vector_t), intent(inout) :: x, y, z
    type(fld_file_data_t) :: fld_data
    integer(kind=MPI_OFFSET_KIND) :: byte_offset
    type(mpi_file) :: fh
    type(mpi_status) :: status
    integer :: n, ierr, lxyz, i, j, e, nd
 
    n = x%size()
    nd = n*fld_data%gdim
    lxyz = fld_data%lx*fld_data%ly*fld_data%lz
 
    if (this%dp_precision) then
       call mpi_file_read_at_all(fh, byte_offset, tmp_dp, nd, &
            mpi_double_precision, status, ierr)
    else
       call mpi_file_read_at_all(fh, byte_offset, tmp_sp, nd, &
            mpi_real, status, ierr)
    end if
 
 
    if (this%dp_precision) then
       i = 1
       do e = 1, fld_data%nelv
          do j = 1, lxyz
             x%x((e-1)*lxyz+j) = tmp_dp(i)
             i = i + 1
          end do
          do j = 1, lxyz
             y%x((e-1)*lxyz+j) = tmp_dp(i)
             i = i + 1
          end do
          if (fld_data%gdim .eq. 3) then
             do j = 1, lxyz
                z%x((e-1)*lxyz+j) = tmp_dp(i)
                i = i + 1
             end do
          end if
       end do
    else
       i = 1
       do e = 1, fld_data%nelv
          do j = 1, lxyz
             x%x((e-1)*lxyz+j) = tmp_sp(i)
             i = i +1
          end do
          do j = 1, lxyz
             y%x((e-1)*lxyz+j) = tmp_sp(i)
             i = i +1
          end do
          if (fld_data%gdim .eq. 3) then
             do j = 1, lxyz
                z%x((e-1)*lxyz+j) = tmp_sp(i)
                i = i +1
             end do
          end if
       end do
    end if
 
    call x%copy_from(host_to_device, sync = .false.)
    call y%copy_from(host_to_device, sync = .false.)
    if (fld_data%gdim .eq. 3) then
       call z%copy_from(host_to_device, sync = .false.)
    end if
 
  end subroutine fld_file_read_vector_field
 
  subroutine fld_file_set_precision(this, precision)
    class(fld_file_t) :: this
    integer, intent(in) :: precision
 
    if (precision .eq. dp) then
       this%dp_precision = .true.
    else if (precision .eq. sp) then
       this%dp_precision = .false.
    else
       call neko_error('Invalid precision')
    end if
 
  end subroutine fld_file_set_precision
 
  subroutine fld_file_set_mask(this, mask)
    class(fld_file_t) :: this
    type(mask_t), intent(inout), optional :: mask
 
    if (present(mask)) then
       call this%mask%init_from_mask(mask)
    else
       call this%mask%free()
    end if
 
  end subroutine fld_file_set_mask
 
  function fld_file_get_fld_fname(this) result(fname)
    class(fld_file_t), intent(in) :: this
    character(len=1024) :: fname
    character(len=1024) :: path, name, id_str
    integer :: suffix_pos
 
    call filename_path(this%get_base_fname(), path)
    call filename_name(this%get_base_fname(), name)
 
    write(fname, '(a,a,a,i5.5)') trim(path), trim(name), &
         '0.f', this%get_counter()
 
  end function fld_file_get_fld_fname
 
  function fld_file_get_meta_fname(this) result(fname)
    class(fld_file_t), intent(in) :: this
    character(len=1024) :: fname
    character(len=1024) :: path, name, id_str
 
    call filename_path(this%get_base_fname(), path)
    call filename_name(this%get_base_fname(), name)
 
    write(id_str, '(i5,a)') this%get_start_counter(), '.nek5000'
    write(fname, '(a,a,a)') trim(path), trim(name), trim(adjustl(id_str))
 
  end function fld_file_get_meta_fname
 
end module fld_file