program chilled_water_decoupler
    implicit none
    integer :: iostat_val
    double precision :: Q_prim_m3h, Q_sec_m3h, T_sup_C, T_ret_des_C, T_ret_act_C
    double precision :: v_bridge_max_ms, Q_bridge_m3h, D_bridge_min_mm, v_actual_bridge
    double precision :: dP_bridge_Pa, T_mixed_sup_C, T_mixed_ret_C, low_dt_penalty_pct
    character(len=32) :: bridge_flow_dir, low_dt_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Q_prim_m3h       ! Primary Chiller Total Flow [m3/h] (e.g. 150.0)
    read(*,*,iostat=iostat_val) Q_sec_m3h        ! Secondary Load Distribution Flow [m3/h] (e.g. 120.0)
    read(*,*,iostat=iostat_val) T_sup_C          ! Chiller Supply Temp [deg C] (e.g. 6.0)
    read(*,*,iostat=iostat_val) T_ret_des_C      ! Design Return Temp [deg C] (e.g. 12.0)
    read(*,*,iostat=iostat_val) T_ret_act_C      ! Actual Secondary Return Temp [deg C] (e.g. 9.5)
    read(*,*,iostat=iostat_val) v_bridge_max_ms  ! Max Allowable Bridge Velocity [m/s] (e.g. 0.50)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for chilled water decoupler sizing.'
        stop
    end if

    Q_bridge_m3h = Q_prim_m3h - Q_sec_m3h

    if (Q_bridge_m3h > 0.5d0) then
        bridge_flow_dir = 'FORWARD FLOW (Surplus to Return)'
        T_mixed_sup_C = T_sup_C
        T_mixed_ret_C = (Q_sec_m3h * T_ret_act_C + Q_bridge_m3h * T_sup_C) / Q_prim_m3h
    else if (Q_bridge_m3h < -0.5d0) then
        bridge_flow_dir = 'REVERSE FLOW (Deficit: Supply Degraded)'
        T_mixed_sup_C = (Q_prim_m3h * T_sup_C + abs(Q_bridge_m3h) * T_ret_act_C) / Q_sec_m3h
        T_mixed_ret_C = T_ret_act_C
    else
        bridge_flow_dir = 'BALANCED PERFECT EQUILIBRIUM'
        T_mixed_sup_C = T_sup_C
        T_mixed_ret_C = T_ret_act_C
    end if

    ! Decoupler Bridge Sizing (sized for single largest chiller or max mismatch flow)
    D_bridge_min_mm = sqrt((4.0d0 * (max(10.0d0, abs(Q_bridge_m3h)) / 3600.0d0)) / &
                           (PI * max(0.2d0, v_bridge_max_ms))) * 1000.0d0

    ! Low Delta T syndrome penalty
    if ((T_ret_act_C - T_sup_C) < (T_ret_des_C - T_sup_C) * 0.85d0) then
        low_dt_status = 'SEVERE LOW DELTA-T SYNDROME'
        low_dt_penalty_pct = (1.0d0 - (T_ret_act_C - T_sup_C) / max(0.5d0, T_ret_des_C - T_sup_C)) * 100.0d0
    else
        low_dt_status = 'ACCEPTABLE TEMPERATURE DELTA'
        low_dt_penalty_pct = 0.0d0
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CHILLED WATER DECOUPLER SIZING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Primary / Secondary Flow  = ', Q_prim_m3h, ' m3/h / ', Q_sec_m3h, ' m3/h'
    write(*,'(A,F10.1,A)')  'Decoupler Bypass Flow     = ', Q_bridge_m3h, ' m3/h'
    write(*,'(A,A)')        'Bridge Flow Direction     = ', trim(bridge_flow_dir)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.1,A)')  'Min Decoupler Bridge Diam = ', D_bridge_min_mm, ' mm (v <= 0.5 m/s)'
    write(*,'(A,F10.2,A)')  'Mixed Header Supply Temp  = ', T_mixed_sup_C, ' C'
    write(*,'(A,F10.2,A)')  'Mixed Header Return Temp  = ', T_mixed_ret_C, ' C'
    write(*,'(A,A)')        'Low Delta-T Status        = ', trim(low_dt_status)
    write(*,'(A,F10.1,A)')  'Plant Capacity Loss / Pen = ', low_dt_penalty_pct, ' %'
    write(*,'(A)') '============================================================'

end program chilled_water_decoupler
