program fluidized_bed_hydrodynamics_umf
    implicit none
    integer :: iostat_val, particle_type
    double precision :: dp_microns, rho_p_kgm3, T_gas_C, P_gas_bar, U0_gas_ms
    double precision :: H_bed_m, D_bed_m, sphericity_phi, voidage_emf
    double precision :: rho_g, mu_g, g, dp_m, Ar_num, Re_mf, Umf_ms, Ut_ms, Umb_ms
    double precision :: dp_bed_kPa, bed_expansion_ratio, fluidization_index
    character(len=16) :: geldart_group
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) particle_type   ! 1=Silica Sand, 2=FCC Catalyst, 3=Coal Ash / Biomass, 4=Glass Beads
    read(*,*,iostat=iostat_val) dp_microns      ! Mean Particle Diameter [microns] (e.g. 250.0)
    read(*,*,iostat=iostat_val) rho_p_kgm3      ! Particle Solid Density [kg/m3] (e.g. 2600.0)
    read(*,*,iostat=iostat_val) T_gas_C         ! Gas Temperature [deg C] (e.g. 200.0)
    read(*,*,iostat=iostat_val) P_gas_bar       ! Bed Operating Pressure [bar abs] (e.g. 1.2)
    read(*,*,iostat=iostat_val) U0_gas_ms       ! Superficial Gas Velocity U0 [m/s] (e.g. 0.35)
    read(*,*,iostat=iostat_val) H_bed_m         ! Static Bed Height [m] (e.g. 1.50)
    read(*,*,iostat=iostat_val) D_bed_m         ! Column Diameter [m] (e.g. 0.80)
    read(*,*,iostat=iostat_val) sphericity_phi  ! Particle Sphericity (e.g. 0.85)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for fluidized bed calculation.'
        stop
    end if

    if (dp_microns <= 0.0d0 .or. rho_p_kgm3 <= 0.0d0 .or. H_bed_m <= 0.0d0) then
        write(*,*) 'ERROR: Particle diameter, density and bed height must be positive.'
        stop
    end if

    g = 9.80665d0
    dp_m = dp_microns * 1.0d-6
    voidage_emf = 0.42d0 ! Standard packed bed voidage at Umf

    ! Air / Flue gas properties at T_gas and P_gas
    rho_g = (P_gas_bar * 1.0d5 * 0.02896d0) / (8.314d0 * (T_gas_C + 273.15d0))
    mu_g = 1.82d-5 * (((T_gas_C + 273.15d0) / 293.15d0)**0.7d0)

    ! Archimedes Number
    Ar_num = (dp_m**3 * rho_g * (rho_p_kgm3 - rho_g) * g) / (mu_g**2)

    ! Wen & Yu Correlation for Re_mf
    Re_mf = sqrt(33.7d0**2 + 0.0408d0 * Ar_num) - 33.7d0
    Umf_ms = (Re_mf * mu_g) / (rho_g * dp_m)

    ! Haider-Levenspiel Terminal Velocity Ut
    Ut_ms = ((mu_g / (rho_g * dp_m)) * ( (18.0d0 / (Ar_num**(2.0d0/3.0d0))) + &
            ((2.335d0 - 1.744d0 * sphericity_phi) / (Ar_num**(1.0d0/6.0d0))) )**(-1.0d0))

    ! Geldart Group Classification
    if (dp_microns < 30.0d0) then
        geldart_group = 'GROUP C (COHESIVE)'
        Umb_ms = Umf_ms
    else if (dp_microns < 100.0d0 .and. rho_p_kgm3 < 1500.0d0) then
        geldart_group = 'GROUP A (AERATABLE)'
        Umb_ms = 33.0d0 * dp_m * ((rho_g / mu_g)**0.1d0)
    else if (dp_microns < 800.0d0) then
        geldart_group = 'GROUP B (BUBBLING)'
        Umb_ms = Umf_ms
    else
        geldart_group = 'GROUP D (SPOUTABLE)'
        Umb_ms = Umf_ms
    end if

    ! Fluidization Quality Index (U0 / Umf)
    fluidization_index = U0_gas_ms / max(1.0d-5, Umf_ms)

    ! Total Bed Pressure Drop at Fluidization [kPa]
    dp_bed_kPa = ((1.0d0 - voidage_emf) * (rho_p_kgm3 - rho_g) * g * H_bed_m) / 1000.0d0

    ! Bed Expansion Ratio
    if (U0_gas_ms > Umf_ms) then
        bed_expansion_ratio = 1.0d0 + 0.20d0 * ((U0_gas_ms - Umf_ms) / max(0.1d0, Ut_ms - Umf_ms))**0.75d0
    else
        bed_expansion_ratio = 1.0d0
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — FLUIDIZED BED HYDRODYNAMICS & U_MF'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Particle Size / Density      = ', dp_microns, ' microns / ', rho_p_kgm3, ' kg/m3'
    write(*,'(A,F10.2,A,F10.1,A)') 'Archimedes Number / Gas Temp = ', Ar_num, ' / ', T_gas_C, ' deg C'
    write(*,'(A,A)')               'Geldart Particle Category    = ', trim(geldart_group)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.4,A)')  'MINIMUM FLUIDIZATION VELOCITY = ', Umf_ms, ' m/s'
    write(*,'(A,F10.3,A)')  'Terminal Particle Velocity Ut = ', Ut_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'TOTAL BED PRESSURE DROP (ΔP)  = ', dp_bed_kPa, ' kPa'
    write(*,'(A,F10.2)')    'Fluidization Ratio (U0 / Umf) = ', fluidization_index
    write(*,'(A,F10.2)')    'Bed Expansion Ratio (H / H0)  = ', bed_expansion_ratio
    write(*,'(A)') '============================================================'

end program fluidized_bed_hydrodynamics_umf
