Thermal Comfort (PMV & PPD)

program thermal_comfort_pmv_ppd
    implicit none
    integer :: iostat_val, iter
    double precision :: Ta_C, Tr_C, vel_ms, RH_pct, met_rate, clo_val
    double precision :: M_Wm2, W_Wm2, Icl_m2KW, fcl_val, Pa_Pa
    double precision :: Tcl_C, Tcl_old, hc_val, hcf_val, hcn_val
    double precision :: HL1, HL2, HL3, HL4, HL_rad, HL_conv, L_thermal, PMV_val, PPD_val
    character(len=32) :: comfort_class

    ! Read inputs
    read(*,*,iostat=iostat_val) Ta_C       ! Air Temperature [deg C] (e.g. 23.0)
    read(*,*,iostat=iostat_val) Tr_C       ! Mean Radiant Temp [deg C] (e.g. 23.5)
    read(*,*,iostat=iostat_val) vel_ms     ! Relative Air Velocity [m/s] (e.g. 0.15)
    read(*,*,iostat=iostat_val) RH_pct     ! Relative Humidity [%] (e.g. 50.0)
    read(*,*,iostat=iostat_val) met_rate   ! Metabolic Rate [met] (1.0 = seated office)
    read(*,*,iostat=iostat_val) clo_val    ! Clothing Insulation [clo] (0.5 = summer, 1.0 = winter)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for PMV/PPD thermal comfort.'
        stop
    end if

    M_Wm2 = met_rate * 58.15d0
    W_Wm2 = 0.0d0
    Icl_m2KW = clo_val * 0.155d0

    ! Clothing area factor fcl
    if (Icl_m2KW <= 0.078d0) then
        fcl_val = 1.00d0 + 1.29d0 * Icl_m2KW
    else
        fcl_val = 1.05d0 + 0.645d0 * Icl_m2KW
    end if

    ! Partial water vapor pressure Pa [Pa]
    Pa_Pa = (RH_pct / 100.0d0) * 10.0d0 * exp(16.6536d0 - 4030.183d0 / (Ta_C + 235.0d0)) * 100.0d0

    ! Iterative solution for clothing temperature Tcl
    Tcl_C = Ta_C + (35.5d0 - Ta_C) / (3.5d0 * (6.45d0 * Icl_m2KW + 0.1d0))
    do iter = 1, 150
        Tcl_old = Tcl_C
        hcn_val = 2.38d0 * (abs(Tcl_C - Ta_C)**0.25d0)
        hcf_val = 12.1d0 * sqrt(max(0.01d0, vel_ms))
        hc_val = max(hcn_val, hcf_val)

        HL_rad = 3.96d-8 * fcl_val * (((Tcl_C + 273.15d0)**4) - ((Tr_C + 273.15d0)**4))
        HL_conv = fcl_val * hc_val * (Tcl_C - Ta_C)
        Tcl_C = 35.7d0 - 0.028d0 * (M_Wm2 - W_Wm2) - Icl_m2KW * (HL_rad + HL_conv)
        if (abs(Tcl_C - Tcl_old) < 1.0d-4) exit
    end do

    ! Heat loss components
    HL1 = 3.05d-3 * (5733.0d0 - 6.99d0 * (M_Wm2 - W_Wm2) - Pa_Pa)
    HL2 = max(0.0d0, 0.42d0 * ((M_Wm2 - W_Wm2) - 58.15d0))
    HL3 = 1.7d-5 * M_Wm2 * (5867.0d0 - Pa_Pa)
    HL4 = 0.0014d0 * M_Wm2 * (34.0d0 - Ta_C)

    L_thermal = (M_Wm2 - W_Wm2) - HL1 - HL2 - HL3 - HL4 - HL_rad - HL_conv
    PMV_val = (0.303d0 * exp(-0.036d0 * M_Wm2) + 0.028d0) * L_thermal
    PPD_val = 100.0d0 - 95.0d0 * exp(-0.03353d0 * (PMV_val**4) - 0.2179d0 * (PMV_val**2))

    if (abs(PMV_val) <= 0.20d0) then
        comfort_class = 'ISO 7730 Category A (PPD <= 6%)'
    else if (abs(PMV_val) <= 0.50d0) then
        comfort_class = 'ISO 7730 Category B (PPD <= 10%)'
    else if (abs(PMV_val) <= 0.70d0) then
        comfort_class = 'ISO 7730 Category C (PPD <= 15%)'
    else
        comfort_class = 'NON-COMPLIANT (Thermal Discomfort)'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — ISO 7730 / ASHRAE 55 PMV & PPD ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Air Temp Ta / Radiant Tr  = ', Ta_C, ' C / ', Tr_C, ' C'
    write(*,'(A,F10.2,A,F6.1,A)')  'Air Velocity / RH         = ', vel_ms, ' m/s / ', RH_pct, ' %'
    write(*,'(A,F10.2,A,F10.2,A)') 'Metabolic Rate / Clothing = ', met_rate, ' met / ', clo_val, ' clo'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Clothing Surface Temp Tcl = ', Tcl_C, ' C'
    write(*,'(A,F10.2,A)')  'Thermal Balance Load L    = ', L_thermal, ' W/m2'
    write(*,'(A,F10.2)')    'PREDICTED MEAN VOTE (PMV) = ', PMV_val
    write(*,'(A,F10.1,A)')  'PREDICTED DISSATISFIED PPD= ', PPD_val, ' %'
    write(*,'(A,A)')        'Comfort Classification    = ', trim(comfort_class)
    write(*,'(A)') '============================================================'

end program thermal_comfort_pmv_ppd