program underfloor_radiant_hvac
    implicit none
    integer :: mode_type, iostat_val
    double precision :: Afloor_m2, W_spacing_mm, R_cover_m2KW, Twin_C, Twout_C, Ti_C, RH_pct
    double precision :: dTw_C, dtheta_H, B_factor, q_flux_Wm2, TF_mean_C, T_dew_C
    double precision :: Pvs_kPa, Pv_kPa, Q_total_W, V_water_Lh
    character(len=32) :: limit_status

    ! Read inputs
    read(*,*,iostat=iostat_val) mode_type       ! 1=Heating, 2=Cooling
    read(*,*,iostat=iostat_val) Afloor_m2       ! Room Floor Area [m2] (e.g. 35.0)
    read(*,*,iostat=iostat_val) W_spacing_mm    ! Pipe Pitch / Spacing [mm] (e.g. 150.0)
    read(*,*,iostat=iostat_val) R_cover_m2KW    ! Floor Covering Resistance [m2.K/W] (e.g. 0.02)
    read(*,*,iostat=iostat_val) Twin_C          ! Water Supply Temp [deg C] (e.g. 38.0 for Heat, 16.0 for Cool)
    read(*,*,iostat=iostat_val) Twout_C         ! Water Return Temp [deg C] (e.g. 33.0 for Heat, 19.0 for Cool)
    read(*,*,iostat=iostat_val) Ti_C            ! Room Air Temp [deg C] (e.g. 21.0)
    read(*,*,iostat=iostat_val) RH_pct          ! Room Relative Humidity [%] (e.g. 50.0)

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

    ! Dew point calculation
    Pvs_kPa = 0.61078d0 * exp((17.27d0 * Ti_C) / (Ti_C + 237.3d0))
    Pv_kPa = (RH_pct / 100.0d0) * Pvs_kPa
    T_dew_C = (237.3d0 * log(Pv_kPa / 0.61078d0)) / (17.27d0 - log(Pv_kPa / 0.61078d0))

    dTw_C = abs(Twin_C - Twout_C)
    if (abs(Twin_C - Ti_C) > 0.5d0 .and. abs(Twout_C - Ti_C) > 0.5d0) then
        dtheta_H = abs(Twin_C - Twout_C) / log(abs(Twin_C - Ti_C) / abs(Twout_C - Ti_C))
    else
        dtheta_H = abs(0.5d0 * (Twin_C + Twout_C) - Ti_C)
    end if

    ! Characteristic emission factor B per EN 1264
    B_factor = 6.5d0 / (1.0d0 + (W_spacing_mm * 1.0d-3 * 0.8d0) + (R_cover_m2KW / 0.05d0))

    if (mode_type == 1) then
        ! Heating Mode
        q_flux_Wm2 = B_factor * (dtheta_H**1.0d0) * 1.85d0
        TF_mean_C = Ti_C + (q_flux_Wm2 / 10.8d0)
        if (TF_mean_C > 29.0d0) then
            limit_status = 'EXCEEDS COMFORT LIMIT (TF > 29 C)'
        else
            limit_status = 'COMPLIANT WITH EN 1264 (TF <= 29 C)'
        end if
    else
        ! Cooling Mode
        q_flux_Wm2 = B_factor * (dtheta_H**1.0d0) * 1.35d0
        TF_mean_C = Ti_C - (q_flux_Wm2 / 6.5d0)
        if (TF_mean_C < T_dew_C + 1.0d0) then
            limit_status = 'RISK OF CONDENSATION (TF < Tdew+1)'
        else
            limit_status = 'SAFE AGAINST CONDENSATION'
        end if
    end if

    Q_total_W = q_flux_Wm2 * Afloor_m2
    V_water_Lh = (Q_total_W * 3600.0d0) / (998.0d0 * 4186.0d0 * max(1.0d0, dTw_C)) * 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — EN 1264 RADIANT HEATING & COOLING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F6.1,A)') 'Floor Area / Spacing      = ', Afloor_m2, ' m2 / ', W_spacing_mm, ' mm'
    write(*,'(A,F10.1,A,F6.1,A)') 'Water Supply / Return     = ', Twin_C, ' C / ', Twout_C, ' C'
    write(*,'(A,F10.1,A,F6.1,A)') 'Room Temp / RH            = ', Ti_C, ' C / ', RH_pct, ' %'
    write(*,'(A,F10.2,A)')  'Room Dew Point T_dew      = ', T_dew_C, ' C'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.1,A)')  'Specific Heat Flux q      = ', q_flux_Wm2, ' W/m2'
    write(*,'(A,F10.2,A)')  'Mean Floor Surface Temp TF= ', TF_mean_C, ' C'
    write(*,'(A,A)')        'Safety / Limit Compliance = ', trim(limit_status)
    write(*,'(A,F10.2,A)')  'Total Output Capacity Q   = ', Q_total_W / 1000.0d0, ' kW'
    write(*,'(A,F10.1,A)')  'Hydronic Water Flow Rate  = ', V_water_Lh, ' L/h'
    write(*,'(A)') '============================================================'

end program underfloor_radiant_hvac