libHDF5.f90 Source File

!! Module for writing sectional and integrated rotor results to HDF5
!! - fullresults.h5: cumulative file for entire simulation (open, append each write, close at end)
!! - results00000.h5, results00001.h5, etc.: one file per timestep (open, write, close)
!!     /rotor_<id>/
!!       /geom/                    geom.nml parameters for this rotor
!!       /integrated/              rotor-level integrated forces
!!         iter         (nt,)
!!         psi          (nt,)      reference azimuth (blade 1)
!!         CL           (nt,)      rotor-level
!!         CD           (nt,)
!!         CLu          (nt,)
!!         CDi          (nt,)
!!         CD0          (nt,)
!!         forceInertial (nt, 3)   rotor force vector [N]
!!         lift         (nt, 3)    rotor lift vector [N]
!!         drag         (nt, 3)    rotor drag vector [N]
!!         liftUnsteady (nt, 3)
!!         dragInduced  (nt, 3)
!!         dragProfile  (nt, 3)
!!         dragUnsteady (nt, 3)
!!       /blade_<id>/
!!         /sectional/
!!           r_R        (ns,)       radial stations, written once at init
!!           iter       (nt,)       timestep index, extendable
!!           psi        (nt,)       azimuth per timestep, extendable
!!           secCL      (nt, ns)    extendable
!!           secCD      (nt, ns)
!!           secCLu     (nt, ns)
!!           secAlpha   (nt, ns)    [deg]
!!           secTheta   (nt, ns)    [deg]
!!           secPhi     (nt, ns)    [deg]
!!           secVel     (nt, ns)
!!           secViz     (nt, ns)
!!           secVix     (nt, ns)
!!           secLift    (nt, ns)    magnitude per section
!!           secDrag    (nt, ns)
!!           secForceInertial (nt, ns)
!!           secLiftInPlane   (nt, ns)
!!           secLiftOutPlane  (nt, ns)
!!           secChord   (ns,)       geometry, written once
!!           secArea    (ns,)       geometry, written once
!!       /integrated/
!!         iter         (nt,)
!!         psi          (nt,)
!!         CL           (nt,)       blade-level
!!         CD           (nt,)       blade-level
!!         CLu          (nt,)       blade-level
!!         CDi          (nt,)       blade-level
!!         CD0          (nt,)       blade-level
!!         secLift      (nt,)       blade-integrated magnitude
!!         secDrag      (nt,)       blade-integrated magnitude
!!         secForceInertial (nt,)   blade-integrated magnitude
!!         totalArea    (nt,)       blade total planform area
!!         forceInertial (nt, 3)    blade force vector [N]
!!         lift         (nt, 3)     blade lift vector [N]
!!         drag         (nt, 3)     blade drag vector [N]
!!         liftUnsteady (nt, 3)     blade unsteady lift [N]
!!         dragInduced  (nt, 3)     blade induced drag [N]
!!         dragProfile  (nt, 3)     blade profile drag [N]
!!         dragUnsteady (nt, 3)     blade unsteady drag [N]
!!         flap         (nt,)       [deg]
!!         theta        (nt,)       [deg]
!!         dflap        (nt,)       [deg]

module libHDF5

  use hdf5
  use libMath, only: dp, radToDeg
  use classdef, only: rotor_class, switches_class
  use libErrorHandling

  implicit none

contains

  ! -----------------------------------------------------------------------
  subroutine init_hdf5(rotorArray, nt, dt, density, velSound, switches, filename, file_id)
    !! Create HDF5 file with extendable structure. Call before results2hdf5, then close_hdf5.
    type(rotor_class), intent(in), dimension(:) :: rotorArray
    integer, intent(in) :: nt
    real(dp), intent(in) :: dt, density, velSound
    type(switches_class), intent(in) :: switches
    character(len=*), intent(in) :: filename
    integer(HID_T), intent(out) :: file_id

    integer :: hdferr, ir, ib, nr
    integer(HID_T) :: grp_rotor, grp_rotor_int, grp_blade, grp_sec, grp_int
    character(len=64) :: grpname
    integer(HSIZE_T), parameter :: CHUNK_T = 64

    nr = size(rotorArray)

    call h5open_f(hdferr)
    call h5fcreate_f(trim(filename), H5F_ACC_TRUNC_F, file_id, hdferr)
    if (hdferr /= 0) then
      call raise_error(ERR_MISC, 'Failed to create HDF5 file: '//trim(filename)// &
        & ' (ensure Results/ directory exists)')
    endif

    call write_params(file_id, nr, nt, dt, density, velSound, switches, rotorArray)

    do ir = 1, nr
      grpname = 'rotor_'//trim(rotorArray(ir)%id)
      call h5gcreate_f(file_id, trim(grpname), grp_rotor, hdferr)
      call write_geom_params(grp_rotor, rotorArray(ir))

      grpname = 'rotor_'//trim(rotorArray(ir)%id)//'/integrated'
      call h5gcreate_f(file_id, trim(grpname), grp_rotor_int, hdferr)
      call create_extendable_1d(grp_rotor_int, 'iter', CHUNK_T)
      call create_extendable_1d(grp_rotor_int, 'psi', CHUNK_T)
      call create_extendable_1d(grp_rotor_int, 'CL', CHUNK_T)
      call create_extendable_1d(grp_rotor_int, 'CD', CHUNK_T)
      call create_extendable_1d(grp_rotor_int, 'CLu', CHUNK_T)
      call create_extendable_1d(grp_rotor_int, 'CDi', CHUNK_T)
      call create_extendable_1d(grp_rotor_int, 'CD0', CHUNK_T)
      call create_extendable_2d(grp_rotor_int, 'forceInertial', CHUNK_T, 3_HSIZE_T)
      call create_extendable_2d(grp_rotor_int, 'lift', CHUNK_T, 3_HSIZE_T)
      call create_extendable_2d(grp_rotor_int, 'drag', CHUNK_T, 3_HSIZE_T)
      call create_extendable_2d(grp_rotor_int, 'liftUnsteady', CHUNK_T, 3_HSIZE_T)
      call create_extendable_2d(grp_rotor_int, 'dragInduced', CHUNK_T, 3_HSIZE_T)
      call create_extendable_2d(grp_rotor_int, 'dragProfile', CHUNK_T, 3_HSIZE_T)
      call create_extendable_2d(grp_rotor_int, 'dragUnsteady', CHUNK_T, 3_HSIZE_T)
      call h5gclose_f(grp_rotor_int, hdferr)

      do ib = 1, rotorArray(ir)%nbConvect
        grpname = 'rotor_'//trim(rotorArray(ir)%id)// &
          & '/blade_'//trim(rotorArray(ir)%blade(ib)%id)
        call h5gcreate_f(file_id, trim(grpname), grp_blade, hdferr)

        grpname = 'rotor_'//trim(rotorArray(ir)%id)// &
          & '/blade_'//trim(rotorArray(ir)%blade(ib)%id)//'/sectional'
        call h5gcreate_f(file_id, trim(grpname), grp_sec, hdferr)
        call write_fixed_1d(grp_sec, 'r_R', &
          rotorArray(ir)%blade(ib)%secChord / rotorArray(ir)%blade(ib)%secChord(1), &
          int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_1d(grp_sec, 'iter', CHUNK_T)
        call create_extendable_1d(grp_sec, 'psi', CHUNK_T)
        call create_extendable_2d(grp_sec, 'secCL', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secCD', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secCLu', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secAlpha', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secTheta', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secPhi', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secVel', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secViz', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secVix', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secLift', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secDrag', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secForceInertial', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secLiftInPlane', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call create_extendable_2d(grp_sec, 'secLiftOutPlane', CHUNK_T, int(rotorArray(ir)%ns, HSIZE_T))
        call write_fixed_1d(grp_sec, 'secChord', rotorArray(ir)%blade(ib)%secChord, &
          int(rotorArray(ir)%ns, HSIZE_T))
        call write_fixed_1d(grp_sec, 'secArea', rotorArray(ir)%blade(ib)%secArea, &
          int(rotorArray(ir)%ns, HSIZE_T))
        call h5gclose_f(grp_sec, hdferr)

        grpname = 'rotor_'//trim(rotorArray(ir)%id)// &
          & '/blade_'//trim(rotorArray(ir)%blade(ib)%id)//'/integrated'
        call h5gcreate_f(file_id, trim(grpname), grp_int, hdferr)
        call create_extendable_1d(grp_int, 'iter', CHUNK_T)
        call create_extendable_1d(grp_int, 'psi', CHUNK_T)
        call create_extendable_1d(grp_int, 'CL', CHUNK_T)
        call create_extendable_1d(grp_int, 'CD', CHUNK_T)
        call create_extendable_1d(grp_int, 'CLu', CHUNK_T)
        call create_extendable_1d(grp_int, 'CDi', CHUNK_T)
        call create_extendable_1d(grp_int, 'CD0', CHUNK_T)
        call create_extendable_1d(grp_int, 'secLift', CHUNK_T)
        call create_extendable_1d(grp_int, 'secDrag', CHUNK_T)
        call create_extendable_1d(grp_int, 'secForceInertial', CHUNK_T)
        call create_extendable_1d(grp_int, 'totalArea', CHUNK_T)
        call create_extendable_2d(grp_int, 'forceInertial', CHUNK_T, 3_HSIZE_T)
        call create_extendable_2d(grp_int, 'lift', CHUNK_T, 3_HSIZE_T)
        call create_extendable_2d(grp_int, 'drag', CHUNK_T, 3_HSIZE_T)
        call create_extendable_2d(grp_int, 'liftUnsteady', CHUNK_T, 3_HSIZE_T)
        call create_extendable_2d(grp_int, 'dragInduced', CHUNK_T, 3_HSIZE_T)
        call create_extendable_2d(grp_int, 'dragProfile', CHUNK_T, 3_HSIZE_T)
        call create_extendable_2d(grp_int, 'dragUnsteady', CHUNK_T, 3_HSIZE_T)
        call create_extendable_1d(grp_int, 'flap', CHUNK_T)
        call create_extendable_1d(grp_int, 'theta', CHUNK_T)
        call create_extendable_1d(grp_int, 'dflap', CHUNK_T)
        call h5gclose_f(grp_int, hdferr)
        call h5gclose_f(grp_blade, hdferr)
      enddo
      call h5gclose_f(grp_rotor, hdferr)
    enddo
  end subroutine init_hdf5

  ! -----------------------------------------------------------------------
  subroutine results2hdf5(rotor, file_id, iter)
    !! Append one timestep of data to open HDF5 file.
    type(rotor_class), intent(in) :: rotor
    integer(HID_T), intent(in) :: file_id
    integer, intent(in) :: iter

    integer :: ib, hdferr, i
    integer(HID_T) :: grp_sec, grp_int, grp_rotor_int
    character(len=64) :: grpname
    real(dp) :: signLiftBlade, bladeDenom, signLiftRotor

    signLiftRotor = sign(1._dp, dot_product(rotor%lift, rotor%zAxisBody))
    grpname = 'rotor_'//trim(rotor%id)//'/integrated'
    call h5gopen_f(file_id, trim(grpname), grp_rotor_int, hdferr)
    call append_scalar_1d(grp_rotor_int, 'iter', real(iter, dp))
    call append_scalar_1d(grp_rotor_int, 'psi', rotor%blade(1)%psi)
    call append_scalar_1d(grp_rotor_int, 'CL', &
      signLiftRotor * norm2(rotor%lift) / rotor%nonDimforceDenominator)
    call append_scalar_1d(grp_rotor_int, 'CD', norm2(rotor%drag) / rotor%nonDimforceDenominator)
    call append_scalar_1d(grp_rotor_int, 'CLu', &
      norm2(rotor%liftUnsteady) / rotor%nonDimforceDenominator)
    call append_scalar_1d(grp_rotor_int, 'CDi', &
      norm2(rotor%dragInduced) / rotor%nonDimforceDenominator)
    call append_scalar_1d(grp_rotor_int, 'CD0', &
      norm2(rotor%dragProfile) / rotor%nonDimforceDenominator)
    call append_row_2d(grp_rotor_int, 'forceInertial', rotor%forceInertial)
    call append_row_2d(grp_rotor_int, 'lift', rotor%lift)
    call append_row_2d(grp_rotor_int, 'drag', rotor%drag)
    call append_row_2d(grp_rotor_int, 'liftUnsteady', rotor%liftUnsteady)
    call append_row_2d(grp_rotor_int, 'dragInduced', rotor%dragInduced)
    call append_row_2d(grp_rotor_int, 'dragProfile', rotor%dragProfile)
    call append_row_2d(grp_rotor_int, 'dragUnsteady', rotor%dragUnsteady)
    call h5gclose_f(grp_rotor_int, hdferr)

    do ib = 1, rotor%nbConvect
      signLiftBlade = sign(1._dp, dot_product(rotor%blade(ib)%lift, rotor%zAxisBody))
      bladeDenom = rotor%nonDimforceDenominator / real(rotor%nbConvect, dp)

      grpname = 'rotor_'//trim(rotor%id)//'/blade_'//trim(rotor%blade(ib)%id)//'/sectional'
      call h5gopen_f(file_id, trim(grpname), grp_sec, hdferr)
      call append_scalar_1d(grp_sec, 'iter', real(iter, dp))
      call append_scalar_1d(grp_sec, 'psi', rotor%blade(ib)%psi)
      call append_row_2d(grp_sec, 'secCL', rotor%blade(ib)%secCL)
      call append_row_2d(grp_sec, 'secCD', rotor%blade(ib)%secCD)
      call append_row_2d(grp_sec, 'secCLu', rotor%blade(ib)%secCLu)
      call append_row_2d(grp_sec, 'secAlpha', rotor%blade(ib)%secAlpha * radToDeg)
      call append_row_2d(grp_sec, 'secTheta', rotor%blade(ib)%secTheta * radToDeg)
      call append_row_2d(grp_sec, 'secPhi', rotor%blade(ib)%secPhi * radToDeg)
      call append_row_2d(grp_sec, 'secVel', &
        [(norm2(rotor%blade(ib)%secChordwiseResVel(:, i)), i=1, rotor%ns)])
      call append_row_2d(grp_sec, 'secViz', rotor%blade(ib)%secViz)
      call append_row_2d(grp_sec, 'secVix', rotor%blade(ib)%secVix)
      call append_row_2d(grp_sec, 'secLift', &
        [(norm2(rotor%blade(ib)%secLift(:, i)), i=1, rotor%ns)])
      call append_row_2d(grp_sec, 'secDrag', &
        [(norm2(rotor%blade(ib)%secDrag(:, i)), i=1, rotor%ns)])
      call append_row_2d(grp_sec, 'secForceInertial', &
        [(norm2(rotor%blade(ib)%secForceInertial(:, i)), i=1, rotor%ns)])
      call append_row_2d(grp_sec, 'secLiftInPlane', &
        [(norm2(rotor%blade(ib)%secLiftInPlane(:, i)), i=1, rotor%ns)])
      call append_row_2d(grp_sec, 'secLiftOutPlane', &
        [(norm2(rotor%blade(ib)%secLiftOutPlane(:, i)), i=1, rotor%ns)])
      call h5gclose_f(grp_sec, hdferr)

      grpname = 'rotor_'//trim(rotor%id)//'/blade_'//trim(rotor%blade(ib)%id)//'/integrated'
      call h5gopen_f(file_id, trim(grpname), grp_int, hdferr)
      call append_scalar_1d(grp_int, 'iter', real(iter, dp))
      call append_scalar_1d(grp_int, 'psi', rotor%blade(ib)%psi)
      call append_scalar_1d(grp_int, 'CL', &
        signLiftBlade * norm2(rotor%blade(ib)%lift) / bladeDenom)
      call append_scalar_1d(grp_int, 'CD', norm2(rotor%blade(ib)%drag) / bladeDenom)
      call append_scalar_1d(grp_int, 'CLu', &
        norm2(rotor%blade(ib)%liftUnsteady) / bladeDenom)
      call append_scalar_1d(grp_int, 'CDi', &
        norm2(rotor%blade(ib)%dragInduced) / bladeDenom)
      call append_scalar_1d(grp_int, 'CD0', &
        norm2(rotor%blade(ib)%dragProfile) / bladeDenom)
      call append_scalar_1d(grp_int, 'secLift', norm2(rotor%blade(ib)%lift))
      call append_scalar_1d(grp_int, 'secDrag', norm2(rotor%blade(ib)%drag))
      call append_scalar_1d(grp_int, 'secForceInertial', norm2(rotor%blade(ib)%forceInertial))
      call append_scalar_1d(grp_int, 'totalArea', sum(rotor%blade(ib)%secArea))
      call append_row_2d(grp_int, 'forceInertial', rotor%blade(ib)%forceInertial)
      call append_row_2d(grp_int, 'lift', rotor%blade(ib)%lift)
      call append_row_2d(grp_int, 'drag', rotor%blade(ib)%drag)
      call append_row_2d(grp_int, 'liftUnsteady', rotor%blade(ib)%liftUnsteady)
      call append_row_2d(grp_int, 'dragInduced', rotor%blade(ib)%dragInduced)
      call append_row_2d(grp_int, 'dragProfile', rotor%blade(ib)%dragProfile)
      call append_row_2d(grp_int, 'dragUnsteady', rotor%blade(ib)%dragUnsteady)
      call append_scalar_1d(grp_int, 'flap', rotor%blade(ib)%flap * radToDeg)
      call append_scalar_1d(grp_int, 'theta', rotor%blade(ib)%theta * radToDeg)
      call append_scalar_1d(grp_int, 'dflap', rotor%blade(ib)%dflap * radToDeg)
      call h5gclose_f(grp_int, hdferr)
    enddo
  end subroutine results2hdf5

  ! -----------------------------------------------------------------------
  subroutine close_hdf5(file_id)
    !! Flush and close HDF5 file. Call after results2hdf5.
    integer(HID_T), intent(in) :: file_id
    integer :: hdferr
    call h5fflush_f(file_id, H5F_SCOPE_LOCAL_F, hdferr)
    call h5fclose_f(file_id, hdferr)
  end subroutine close_hdf5

  ! -----------------------------------------------------------------------
  subroutine deinit_hdf5(file_id)
    !! Deinitialize HDF5 library (h5close_f). file_id unused (kept for API).
    integer(HID_T), intent(in) :: file_id
    integer :: hdferr
    call h5close_f(hdferr)
  end subroutine deinit_hdf5

  ! =======================================================================
  ! Private helper routines
  ! =======================================================================

  subroutine write_geom_params(grp_rotor_id, rotor)
    !! Write geom.nml parameters to geom sub-group under the rotor
    integer(HID_T), intent(in) :: grp_rotor_id
    type(rotor_class), intent(in) :: rotor

    integer :: hdferr
    integer(HID_T) :: grp_params

    call h5gcreate_f(grp_rotor_id, 'geom', grp_params, hdferr)

    ! SURFACE
    call write_scalar_int(grp_params, 'surfaceType', rotor%surfaceType)
    call write_scalar_int(grp_params, 'imagePlane', rotor%imagePlane)
    call write_scalar_int(grp_params, 'imageRotorNum', rotor%imageRotorNum)

    ! PANELS
    call write_scalar_int(grp_params, 'nb', rotor%nb)
    call write_scalar_int(grp_params, 'propConvention', rotor%propConvention)
    call write_scalar_int(grp_params, 'spanSpacing', rotor%spanSpacing)
    call write_scalar_int(grp_params, 'chordSpacing', rotor%chordSpacing)
    call write_scalar_int(grp_params, 'nCamberFiles', rotor%nCamberFiles)
    call write_scalar_int(grp_params, 'nc', rotor%nc)
    call write_scalar_int(grp_params, 'ns', rotor%ns)
    call write_scalar_int(grp_params, 'nNwake', rotor%nNwake)

    ! ORIENT
    call write_vector_3d(grp_params, 'hubCoords', rotor%hubCoords)
    call write_vector_3d(grp_params, 'cgCoords', rotor%cgCoords)
    call write_vector_3d(grp_params, 'fromCoords', rotor%fromCoords)
    call write_vector_3d(grp_params, 'phiThetaPsi', rotor%pts)

    ! GEOMPARAMS
    call write_scalar_double(grp_params, 'radius', rotor%radius)
    call write_scalar_double(grp_params, 'root_cut', rotor%root_cut)
    call write_scalar_double(grp_params, 'chord', rotor%chord)
    call write_scalar_double(grp_params, 'preconeAngle', rotor%preconeAngle*radToDeg)
    call write_scalar_double(grp_params, 'Omega', rotor%Omega)
    call write_vector_3d(grp_params, 'shaftAxis', rotor%shaftAxis)
    call write_scalar_double(grp_params, 'theta0', rotor%controlPitch(1)*radToDeg)
    call write_scalar_double(grp_params, 'thetaC', rotor%controlPitch(2)*radToDeg)
    call write_scalar_double(grp_params, 'thetaS', rotor%controlPitch(3)*radToDeg)
    call write_scalar_double(grp_params, 'thetaTwist', rotor%thetaTwist*radToDeg)
    call write_scalar_int(grp_params, 'ductSwitch', rotor%ductSwitch)
    call write_scalar_int(grp_params, 'axisymmetrySwitch', rotor%axisymmetrySwitch)
    call write_scalar_double(grp_params, 'pivotLE', rotor%pivotLE)
    call write_scalar_double(grp_params, 'flapHinge', rotor%flapHinge)
    call write_scalar_int(grp_params, 'spanwiseLiftSwitch', rotor%spanwiseLiftSwitch)
    call write_scalar_int(grp_params, 'symmetricTau', rotor%symmetricTau)
    call write_scalar_int(grp_params, 'customTrajectorySwitch', rotor%customTrajectorySwitch)
    call write_vector_3d(grp_params, 'velBody', rotor%velBody)
    call write_vector_3d(grp_params, 'omegaBody', rotor%omegaBody)
    call write_scalar_int(grp_params, 'forceCalcSwitch', rotor%forceCalcSwitch)
    call write_scalar_int(grp_params, 'nAirfoils', rotor%nAirfoils)

    ! WAKEPARAMS
    call write_scalar_double(grp_params, 'apparentViscCoeff', rotor%apparentViscCoeff)
    call write_scalar_double(grp_params, 'decayCoeff', rotor%decayCoeff)
    call write_scalar_double(grp_params, 'wakeDispLimitFactor', rotor%wakeDispLimitFactor)
    call write_scalar_int(grp_params, 'wakeTruncateNt', rotor%wakeTruncateNt)
    call write_scalar_int(grp_params, 'prescWakeAfterTruncNt', rotor%prescWakeAfterTruncNt)
    call write_scalar_int(grp_params, 'prescWakeGenNt', rotor%prescWakeGenNt)
    call write_scalar_double(grp_params, 'spanwiseCore', rotor%spanwiseCore)
    call write_scalar_double(grp_params, 'rollupStartRadius', rotor%rollupStartRadius)
    call write_scalar_double(grp_params, 'rollupEndRadius', rotor%rollupEndRadius)
    call write_scalar_double(grp_params, 'initWakeVel', rotor%initWakeVel)
    call write_scalar_double(grp_params, 'psiStart', rotor%psiStart*radToDeg)
    call write_scalar_double(grp_params, 'skewLimit', rotor%skewLimit)

    ! DYNAMICS
    call write_scalar_int(grp_params, 'bladeDynamicsSwitch', rotor%bladeDynamicsSwitch)
    call write_scalar_double(grp_params, 'flapInitial', rotor%flapInitial)
    call write_scalar_double(grp_params, 'dflapInitial', rotor%dflapInitial)
    call write_scalar_double(grp_params, 'Iflap', rotor%Iflap)
    call write_scalar_double(grp_params, 'cflap', rotor%cflap)
    call write_scalar_double(grp_params, 'kflap', rotor%kflap)
    call write_scalar_double(grp_params, 'MflapConstant', rotor%MflapConstant)
    call write_scalar_double(grp_params, 'flap0', rotor%flap0*radToDeg)
    call write_scalar_double(grp_params, 'flapC', rotor%flapC*radToDeg)
    call write_scalar_double(grp_params, 'flapS', rotor%flapS*radToDeg)
    call write_scalar_int(grp_params, 'pitchDynamicsSwitch', rotor%pitchDynamicsSwitch)
    call write_scalar_double(grp_params, 'dpitch', rotor%dpitch)
    call write_scalar_int(grp_params, 'bodyDynamicsSwitch', rotor%bodyDynamicsSwitch)
    call write_scalar_int(grp_params, 'bodyDynamicsIOVars', rotor%bodyDynamicsIOVars)

    ! WINDFRAME
    call write_vector_3d(grp_params, 'dragUnitVec', rotor%dragUnitVec)
    call write_vector_3d(grp_params, 'sideUnitVec', rotor%sideUnitVec)
    call write_vector_3d(grp_params, 'liftUnitVec', rotor%liftUnitVec)

    ! Derived / computed
    call write_scalar_double(grp_params, 'nonDimForceDenom', rotor%nonDimforceDenominator)
    if (rotor%nAirfoils > 0) then
      call write_scalar_double(grp_params, 'alpha0', rotor%alpha0(1))
    else
      call write_scalar_double(grp_params, 'alpha0', 0._dp)
    endif

    call h5gclose_f(grp_params, hdferr)
  end subroutine write_geom_params

  subroutine write_params(file_id, nr, nt, dt, density, velSound, switches, rotorArray)
    integer(HID_T), intent(in) :: file_id
    integer, intent(in) :: nr, nt
    real(dp), intent(in) :: dt, density, velSound
    type(switches_class), intent(in) :: switches
    type(rotor_class), intent(in), dimension(:) :: rotorArray

    integer :: hdferr
    integer(HID_T) :: grp_params

    call h5gcreate_f(file_id, 'params', grp_params, hdferr)

    ! Config (matching params2file)
    call write_scalar_int(grp_params, 'nt', nt)
    call write_scalar_double(grp_params, 'dt', dt)
    call write_scalar_int(grp_params, 'nr', nr)
    call write_scalar_int(grp_params, 'restartWriteNt', switches%restartWriteNt)
    call write_scalar_int(grp_params, 'restartFromNt', switches%restartFromNt)
    call write_scalar_int(grp_params, 'hdf5WriteNt', switches%hdf5WriteNt)
    call write_scalar_int(grp_params, 'ntSub', switches%ntSub)
    call write_scalar_int(grp_params, 'ntSubInit', switches%ntSubInit)
    call write_scalar_double(grp_params, 'density', density)
    call write_scalar_double(grp_params, 'velSound', velSound)
    call write_scalar_int(grp_params, 'wakePlot', switches%wakePlot)
    call write_scalar_int(grp_params, 'wakeTipPlot', switches%wakeTipPlot)
    call write_scalar_int(grp_params, 'rotorForcePlot', switches%rotorForcePlot)
    call write_scalar_int(grp_params, 'gridPlot', switches%gridPlot)
    call write_scalar_int(grp_params, 'wakeDissipation', switches%wakeDissipation)
    call write_scalar_int(grp_params, 'wakeStrain', switches%wakeStrain)
    call write_scalar_int(grp_params, 'wakeBurst', switches%wakeBurst)
    call write_scalar_int(grp_params, 'slowStart', switches%slowStart)
    call write_scalar_int(grp_params, 'slowStartNt', switches%slowStartNt)
    call write_scalar_int(grp_params, 'fdScheme', switches%fdScheme)
    call write_scalar_int(grp_params, 'initWakeVelNt', switches%initWakeVelNt)
    call write_scalar_int(grp_params, 'probe', switches%probe)

    ! Per-rotor geom params are in /rotor_<id>/geom/

    call h5gclose_f(grp_params, hdferr)
  end subroutine write_params

  subroutine write_scalar_int(grp_id, dsetname, val)
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    integer, intent(in) :: val

    integer(HID_T) :: dspace_id, dset_id
    integer(HSIZE_T), dimension(1) :: dims
    integer :: hdferr

    dims(1) = 1
    call h5screate_simple_f(1, dims, dspace_id, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_INTEGER, dspace_id, dset_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, [val], dims, hdferr)
    call h5dclose_f(dset_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine write_scalar_int

  subroutine write_scalar_double(grp_id, dsetname, val)
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    real(dp), intent(in) :: val

    integer(HID_T) :: dspace_id, dset_id
    integer(HSIZE_T), dimension(1) :: dims
    integer :: hdferr

    dims(1) = 1
    call h5screate_simple_f(1, dims, dspace_id, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, [val], dims, hdferr)
    call h5dclose_f(dset_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine write_scalar_double

  subroutine write_vector_3d(grp_id, dsetname, vec)
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    real(dp), intent(in), dimension(3) :: vec

    integer(HID_T) :: dspace_id, dset_id
    integer(HSIZE_T), dimension(1) :: dims
    integer :: hdferr

    dims(1) = 3
    call h5screate_simple_f(1, dims, dspace_id, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, vec, dims, hdferr)
    call h5dclose_f(dset_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine write_vector_3d

  subroutine write_fixed_1d(grp_id, dsetname, data, n)
    !! Write a fixed-size 1D dataset (not extendable) -- used for r_R
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    real(dp), intent(in), dimension(:) :: data
    integer(HSIZE_T), intent(in) :: n

    integer(HID_T) :: dspace_id, dset_id
    integer(HSIZE_T), dimension(1) :: dims
    integer :: hdferr

    dims(1) = n
    call h5screate_simple_f(1, dims, dspace_id, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_DOUBLE, &
      dspace_id, dset_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, dims, hdferr)
    call h5dclose_f(dset_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine write_fixed_1d

  subroutine write_fixed_2d(grp_id, dsetname, data, nrows, ncols)
    !! Write a fixed-size 2D dataset (1 x ncols) for snapshot files
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    real(dp), intent(in), dimension(:) :: data
    integer(HSIZE_T), intent(in) :: nrows, ncols

    integer(HID_T) :: dspace_id, dset_id
    integer(HSIZE_T), dimension(2) :: dims
    integer :: hdferr

    dims = [ncols, nrows]  ! Fortran column-major order
    call h5screate_simple_f(2, dims, dspace_id, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_DOUBLE, &
      dspace_id, dset_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, data, dims, hdferr)
    call h5dclose_f(dset_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine write_fixed_2d

  ! -----------------------------------------------------------------------
  subroutine create_extendable_1d(grp_id, dsetname, chunk_size)
    !! Create an extendable 1D dataset of shape (0,) with unlimited max dim
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    integer(HSIZE_T), intent(in) :: chunk_size

    integer(HID_T) :: dspace_id, dset_id, dcpl_id
    integer(HSIZE_T), dimension(1) :: dims, maxdims, chunk_dims
    integer :: hdferr

    dims(1) = 0
    maxdims(1) = H5S_UNLIMITED_F
    chunk_dims(1) = chunk_size

    call h5screate_simple_f(1, dims, dspace_id, hdferr, maxdims)
    call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl_id, hdferr)
    call h5pset_chunk_f(dcpl_id, 1, chunk_dims, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_DOUBLE, &
      dspace_id, dset_id, hdferr, dcpl_id)
    call h5dclose_f(dset_id, hdferr)
    call h5pclose_f(dcpl_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine create_extendable_1d

  ! -----------------------------------------------------------------------
  subroutine create_extendable_2d(grp_id, dsetname, chunk_t, ns)
    !! Create an extendable 2D dataset of shape (0, ns) -- time x span
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    integer(HSIZE_T), intent(in) :: chunk_t, ns

    integer(HID_T) :: dspace_id, dset_id, dcpl_id
    integer(HSIZE_T), dimension(2) :: dims, maxdims, chunk_dims
    integer :: hdferr

    dims = [0_HSIZE_T, ns]
    maxdims = [H5S_UNLIMITED_F, ns]
    chunk_dims = [chunk_t, ns]

    call h5screate_simple_f(2, dims, dspace_id, hdferr, maxdims)
    call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl_id, hdferr)
    call h5pset_chunk_f(dcpl_id, 2, chunk_dims, hdferr)
    call h5dcreate_f(grp_id, dsetname, H5T_NATIVE_DOUBLE, &
      dspace_id, dset_id, hdferr, dcpl_id)
    call h5dclose_f(dset_id, hdferr)
    call h5pclose_f(dcpl_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
  end subroutine create_extendable_2d

  ! -----------------------------------------------------------------------
  subroutine append_scalar_1d(grp_id, dsetname, val)
    !! Append a single scalar to an extendable 1D dataset
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    real(dp), intent(in) :: val

    integer(HID_T) :: dset_id, dspace_id, memspace_id
    integer(HSIZE_T), dimension(1) :: cur_dims, new_dims, offset, count
    integer :: hdferr

    call h5dopen_f(grp_id, dsetname, dset_id, hdferr)
    call h5dget_space_f(dset_id, dspace_id, hdferr)
    call h5sget_simple_extent_dims_f(dspace_id, cur_dims, new_dims, hdferr)

    new_dims(1) = cur_dims(1) + 1
    call h5dset_extent_f(dset_id, new_dims, hdferr)
    call h5dget_space_f(dset_id, dspace_id, hdferr)

    offset(1) = cur_dims(1)
    count(1) = 1
    call h5sselect_hyperslab_f(dspace_id, H5S_SELECT_SET_F, &
      offset, count, hdferr)
    call h5screate_simple_f(1, count, memspace_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, [val], count, hdferr, &
      memspace_id, dspace_id)

    call h5sclose_f(memspace_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
    call h5dclose_f(dset_id, hdferr)
  end subroutine append_scalar_1d

  ! -----------------------------------------------------------------------
  subroutine append_row_2d(grp_id, dsetname, row)
    !! Append a 1D array as a new row in an extendable 2D dataset (nt x ns)
    integer(HID_T), intent(in) :: grp_id
    character(len=*), intent(in) :: dsetname
    real(dp), intent(in), dimension(:) :: row

    integer(HID_T) :: dset_id, dspace_id, memspace_id
    integer(HSIZE_T), dimension(2) :: cur_dims, max_dims, new_dims, offset, count
    integer :: hdferr
    integer(HSIZE_T) :: ns

    ns = int(size(row), HSIZE_T)

    call h5dopen_f(grp_id, dsetname, dset_id, hdferr)
    call h5dget_space_f(dset_id, dspace_id, hdferr)
    call h5sget_simple_extent_dims_f(dspace_id, cur_dims, max_dims, hdferr)

    new_dims = [cur_dims(1) + 1, ns]
    call h5dset_extent_f(dset_id, new_dims, hdferr)
    call h5dget_space_f(dset_id, dspace_id, hdferr)

    offset = [cur_dims(1), 0_HSIZE_T]
    count  = [1_HSIZE_T, ns]
    call h5sselect_hyperslab_f(dspace_id, H5S_SELECT_SET_F, &
      offset, count, hdferr)
    call h5screate_simple_f(2, count, memspace_id, hdferr)
    call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, row, count, hdferr, &
      memspace_id, dspace_id)

    call h5sclose_f(memspace_id, hdferr)
    call h5sclose_f(dspace_id, hdferr)
    call h5dclose_f(dset_id, hdferr)
  end subroutine append_row_2d

end module libHDF5