Hydrogen Cryogenic Liquefaction (Claude)

program hydrogen_liquefaction_claude_cycle
    implicit none
    integer :: iostat_val
    double precision :: P_high_bar, T_feed_K, m_dot_feed_kgh, eta_comp, eta_turb, T_ln2_K
    double precision :: m_dot_feed_kgs, m_dot_lh2_kgh, y_liq_frac, W_comp_kW, W_turb_kW, W_net_kW
    double precision :: sec_kWh_kg, w_carnot_kWh_kg, eta_exergy_pct, Q_ln2_kW, para_fraction_pct

    ! Read inputs
    read(*,*,iostat=iostat_val) P_high_bar       ! Compressor Discharge Pressure [bar] (e.g. 25.0)
    read(*,*,iostat=iostat_val) T_feed_K         ! Ambient Feed Temp [K] (e.g. 295.0)
    read(*,*,iostat=iostat_val) m_dot_feed_kgh   ! Hydrogen Feed Flow [kg/h] (e.g. 500.0)
    read(*,*,iostat=iostat_val) eta_comp         ! Compressor Isentropic Efficiency (e.g. 0.80)
    read(*,*,iostat=iostat_val) eta_turb         ! Cryogenic Turbine Efficiency (e.g. 0.86)
    read(*,*,iostat=iostat_val) T_ln2_K          ! LN2 Precooling Temperature [K] (e.g. 80.0)

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

    if (P_high_bar <= 5.0d0 .or. m_dot_feed_kgh <= 0.0d0 .or. eta_comp <= 0.0d0) then
        write(*,*) 'ERROR: Pressure and mass flow must be positive.'
        stop
    end if

    m_dot_feed_kgs = m_dot_feed_kgh / 3600.0d0

    ! Claude Precooled Cycle Liquefaction Fraction (Huang & Barron model)
    y_liq_frac = 0.22d0 * (1.0d0 + 0.15d0 * (P_high_bar / 25.0d0)) * (eta_turb / 0.85d0) * (80.0d0 / T_ln2_K)
    if (y_liq_frac > 0.38d0) y_liq_frac = 0.38d0
    if (y_liq_frac < 0.10d0) y_liq_frac = 0.10d0

    m_dot_lh2_kgh = m_dot_feed_kgh * y_liq_frac

    ! Multi-stage compression work (ambient to P_high): 3 stages
    W_comp_kW = (3.0d0 * (m_dot_feed_kgs) * (8.314d0 * 295.0d0 / 0.002016d0) * 1.4d0 / 0.4d0 * &
                (((P_high_bar / 1.0d0)**(0.4d0 / (3.0d0 * 1.4d0))) - 1.0d0)) / (1000.0d0 * eta_comp)

    ! Cryogenic Turbine Work Recovery
    W_turb_kW = (0.75d0 * m_dot_feed_kgs * 14.3d0 * (T_ln2_K - 30.0d0) * eta_turb)

    W_net_kW = W_comp_kW - W_turb_kW

    ! Specific Energy Consumption [kWh / kg LH2]
    sec_kWh_kg = W_net_kW / (m_dot_lh2_kgh)
    w_carnot_kWh_kg = 3.92d0 ! Thermodynamic minimum ideal work for H2 liquefaction + ortho-para conversion
    eta_exergy_pct = (w_carnot_kWh_kg / max(3.92d0, sec_kWh_kg)) * 100.0d0

    ! Liquid Nitrogen (LN2) cooling duty
    Q_ln2_kW = m_dot_feed_kgs * 14.3d0 * (T_feed_K - T_ln2_K)
    para_fraction_pct = 99.8d0 ! Equilibrium para-hydrogen concentration at 20.3 K

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — HYDROGEN LIQUEFACTION (PRECOOLED CLAUDE)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Discharge Press / LN2 Precool= ', P_high_bar, ' bar / ', T_ln2_K, ' K'
    write(*,'(A,F10.1,A,F10.2,A)') 'Feed Flow / Liq Yield Frac  = ', m_dot_feed_kgh, ' kg/h / ', y_liq_frac*100.0d0, ' %'
    write(*,'(A,F10.2,A)')         'Liquid H2 Production Rate  = ', m_dot_lh2_kgh, ' kg/h'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'SPECIFIC ENERGY CONSUMPTION= ', sec_kWh_kg, ' kWh/kg LH2'
    write(*,'(A,F10.2,A)')  'EXERGY EFFICIENCY          = ', eta_exergy_pct, ' % of Carnot'
    write(*,'(A,F10.2,A)')  'Net Electrical Power       = ', W_net_kW, ' kW'
    write(*,'(A,F10.2,A)')  'Compressor Power           = ', W_comp_kW, ' kW'
    write(*,'(A,F10.2,A)')  'LN2 Precooling Thermal Duty= ', Q_ln2_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Para-Hydrogen Purity       = ', para_fraction_pct, ' % (Catalytic Bed)'
    write(*,'(A)') '============================================================'

end program hydrogen_liquefaction_claude_cycle