program brazed_plate_supercritical_co2_bphe
    implicit none
    integer :: iostat_val, cycle_app, num_plates
    double precision :: P_sco2_bar, T_sco2_in_C, m_dot_sco2_kgs, T_cool_in_C, m_dot_cool_kgs
    double precision :: W_plate_m, L_plate_m, chevron_deg, corrugation_depth_mm
    double precision :: rho_b, cp_b, mu_b, k_b, T_pc_C, h_in_kJkg, h_out_kJkg
    double precision :: D_h_m, A_plate_m2, A_total_m2, V_sco2, Re_sco2, Pr_sco2, Nu_sco2
    double precision :: h_sco2, h_cool, U_overall, Q_duty_kW, T_sco2_out_C, T_cool_out_C
    double precision :: dp_sco2_kPa, dp_cool_kPa, compactness_m2m3
    character(len=32) :: pc_regime
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) cycle_app           ! 1=Transcritical CO2 Heat Pump Gas Cooler, 2=sCO2 Brayton Recuperator, 3=sCO2 Waste Heat Exchanger
    read(*,*,iostat=iostat_val) P_sco2_bar          ! sCO2 Operating Pressure [bar] (e.g. 90.0)
    read(*,*,iostat=iostat_val) T_sco2_in_C         ! sCO2 Inlet Temperature [deg C] (e.g. 110.0)
    read(*,*,iostat=iostat_val) m_dot_sco2_kgs      ! sCO2 Mass Flow Rate [kg/s] (e.g. 2.50)
    read(*,*,iostat=iostat_val) T_cool_in_C         ! Coolant Water Inlet Temp [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) m_dot_cool_kgs      ! Coolant Water Mass Flow [kg/s] (e.g. 4.00)
    read(*,*,iostat=iostat_val) num_plates          ! Number of Brazed Plates (e.g. 60)
    read(*,*,iostat=iostat_val) W_plate_m           ! Plate Width W [m] (e.g. 0.25)
    read(*,*,iostat=iostat_val) L_plate_m           ! Plate Length L [m] (e.g. 0.60)
    read(*,*,iostat=iostat_val) chevron_deg         ! Chevron Corrugation Angle [deg] (e.g. 60.0)

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

    if (P_sco2_bar < 73.8d0 .or. m_dot_sco2_kgs <= 0.0d0 .or. num_plates <= 4) then
        write(*,*) 'ERROR: Pressure must be supercritical (P > 73.8 bar) and flow rates positive.'
        stop
    end if

    ! Pseudo-critical temperature correlation: Tpc ~ 31.04 + 0.35 * (P - 73.8)
    T_pc_C = 31.04d0 + 0.35d0 * (P_sco2_bar - 73.8d0)

    if (T_sco2_in_C > T_pc_C .and. T_cool_in_C < T_pc_C) then
        pc_regime = 'TRANS-PSEUDOCRITICAL COOLING'
    else if (T_sco2_in_C > T_pc_C) then
        pc_regime = 'SUPERCRITICAL GAS-LIKE'
    else
        pc_regime = 'SUPERCRITICAL LIQUID-LIKE'
    end if

    ! Supercritical thermophysical properties near 80 - 120 bar
    rho_b = 320.0d0 + 4.5d0 * (P_sco2_bar - 73.8d0) - 2.8d0 * (T_sco2_in_C - 40.0d0)
    rho_b = max(180.0d0, min(850.0d0, rho_b))

    cp_b = 3200.0d0 + 15000.0d0 / (1.0d0 + 0.05d0 * (T_sco2_in_C - T_pc_C)**2) ! Cp peak near Tpc
    mu_b = 0.000035d0
    k_b = 0.065d0

    corrugation_depth_mm = 2.4d0
    D_h_m = 2.0d0 * (corrugation_depth_mm / 1000.0d0)
    A_plate_m2 = W_plate_m * L_plate_m * 1.18d0 ! Enlargement factor 1.18
    A_total_m2 = dble(num_plates - 2) * A_plate_m2

    V_sco2 = m_dot_sco2_kgs / (rho_b * (dble(num_plates/2) * W_plate_m * (corrugation_depth_mm/1000.0d0)))
    Re_sco2 = (rho_b * V_sco2 * D_h_m) / mu_b
    Pr_sco2 = (cp_b * mu_b) / k_b

    ! Chevron Supercritical Nusselt correlation
    Nu_sco2 = 0.30d0 * (Re_sco2**0.68d0) * (Pr_sco2**0.40d0) * (sin(chevron_deg * PI / 180.0d0)**0.3d0)
    h_sco2 = Nu_sco2 * k_b / D_h_m

    h_cool = 5500.0d0 ! Water side high chevron convection
    U_overall = 1.0d0 / ((1.0d0 / h_sco2) + (0.0004d0 / 16.0d0) + (1.0d0 / h_cool))

    ! Heat duty and temperatures
    Q_duty_kW = (m_dot_sco2_kgs * cp_b * (T_sco2_in_C - (T_cool_in_C + 8.0d0))) / 1000.0d0 * 0.85d0
    T_sco2_out_C = T_sco2_in_C - (Q_duty_kW * 1000.0d0 / (m_dot_sco2_kgs * cp_b))
    T_cool_out_C = T_cool_in_C + (Q_duty_kW * 1000.0d0 / (m_dot_cool_kgs * 4184.0d0))

    ! Pressure drops [kPa]
    dp_sco2_kPa = (0.65d0 * (L_plate_m / D_h_m) * (rho_b * V_sco2**2 / 2.0d0)) / 1000.0d0
    dp_cool_kPa = 35.0d0

    compactness_m2m3 = A_total_m2 / (W_plate_m * L_plate_m * (dble(num_plates) * 0.0028d0))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SUPERCRITICAL CO2 BRAZED PLATE (BPHE)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.2,A)') 'Operating Pressure / T_pc    = ', P_sco2_bar, ' bar / ', T_pc_C, ' deg C'
    write(*,'(A,A)')               'Supercritical State Regime   = ', trim(pc_regime)
    write(*,'(A,F10.1,A,I4,A)')    'Total Plate Area / Count     = ', A_total_m2, ' m2 (', num_plates, ' plates)'
    write(*,'(A,F10.0,A)')         'Area Compactness Density     = ', compactness_m2m3, ' m2/m3 (>1000 compact)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'GAS COOLER THERMAL DUTY (Q)  = ', Q_duty_kW, ' kW'
    write(*,'(A,F10.1,A)')  'sCO2 Heat Transfer Coeff h   = ', h_sco2, ' W/(m2.K)'
    write(*,'(A,F10.1,A)')  'Overall Heat Transfer Coeff U= ', U_overall, ' W/(m2.K)'
    write(*,'(A,F10.2,A)')  'sCO2 Outlet Temperature      = ', T_sco2_out_C, ' deg C'
    write(*,'(A,F10.2,A)')  'Coolant Water Outlet Temp    = ', T_cool_out_C, ' deg C'
    write(*,'(A,F10.2,A)')  'sCO2 Channel Pressure Drop   = ', dp_sco2_kPa, ' kPa'
    write(*,'(A)') '============================================================'

end program brazed_plate_supercritical_co2_bphe
