Indirect Evaporative Cooler Sizer (M-Cycle)

program indirect_evaporative_cooler_mcycle
    implicit none
    integer :: iostat_val, climate_type
    double precision :: T_db_in_C, RH_in_pct, V_dot_product_m3h, working_fraction_pct, eps_dp_target
    double precision :: P_atm_kPa, P_vsat_kPa, P_v_kPa, w_in_kgkg, T_wb_in_C, T_dp_in_C
    double precision :: T_db_out_C, eps_wb_pct, Q_cooling_kW, m_dot_product_kgs, m_dot_working_kgs
    double precision :: water_evap_Lph, fan_power_kW, COP_energy
    character(len=32) :: cooling_status
    double precision, parameter :: CP_AIR = 1006.0d0

    ! Read inputs
    read(*,*,iostat=iostat_val) climate_type          ! 1=Hot & Arid Desert (40C, 20% RH), 2=Mediterranean Summer (34C, 40% RH), 3=Warm Moderate (30C, 50% RH)
    read(*,*,iostat=iostat_val) T_db_in_C             ! Ambient Inlet Dry-Bulb Temp [deg C] (e.g. 38.0)
    read(*,*,iostat=iostat_val) RH_in_pct             ! Ambient Inlet Relative Humidity [%] (e.g. 25.0)
    read(*,*,iostat=iostat_val) V_dot_product_m3h     ! Supply Product Air Flow [m3/h] (e.g. 2500.0)
    read(*,*,iostat=iostat_val) working_fraction_pct  ! M-Cycle Working Air Fraction [%] (e.g. 35.0)
    read(*,*,iostat=iostat_val) eps_dp_target         ! Dew-Point Effectiveness [0.60 to 0.90] (e.g. 0.80)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for M-Cycle evaporative cooler calculation.'
        stop
    end if

    if (T_db_in_C <= 0.0d0 .or. RH_in_pct <= 0.0d0 .or. V_dot_product_m3h <= 0.0d0) then
        write(*,*) 'ERROR: Temperature, humidity, and airflow must be positive.'
        stop
    end if

    P_atm_kPa = 101.325d0

    ! Psychrometric calculations
    P_vsat_kPa = 0.61078d0 * exp((17.27d0 * T_db_in_C) / (T_db_in_C + 237.3d0))
    P_v_kPa = (RH_in_pct / 100.0d0) * P_vsat_kPa
    w_in_kgkg = 0.622d0 * P_v_kPa / (P_atm_kPa - P_v_kPa)

    ! Dew point temperature (Magnus formula)
    T_dp_in_C = (237.3d0 * log(max(0.001d0, P_v_kPa / 0.61078d0))) / (17.27d0 - log(max(0.001d0, P_v_kPa / 0.61078d0)))

    ! Wet bulb approximation (Stull formula)
    T_wb_in_C = T_db_in_C * atan(0.151977d0 * sqrt(RH_in_pct + 8.313659d0)) + &
                atan(T_db_in_C + RH_in_pct) - atan(RH_in_pct - 1.676331d0) + &
                0.00391838d0 * (RH_in_pct**1.5d0) * atan(0.023101d0 * RH_in_pct) - 4.686035d0

    ! Product Dry Air Outlet Temp via M-Cycle Dew-Point Effectiveness
    T_db_out_C = T_db_in_C - eps_dp_target * (T_db_in_C - T_dp_in_C)

    ! Wet-bulb effectiveness
    eps_wb_pct = ((T_db_in_C - T_db_out_C) / max(0.5d0, T_db_in_C - T_wb_in_C)) * 100.0d0

    ! Mass flows and sensible cooling capacity
    m_dot_product_kgs = (V_dot_product_m3h * 1.18d0) / 3600.0d0
    m_dot_working_kgs = m_dot_product_kgs * (working_fraction_pct / 100.0d0)

    Q_cooling_kW = (m_dot_product_kgs * CP_AIR * (T_db_in_C - T_db_out_C)) / 1000.0d0

    ! Water consumption rate [L/h]
    water_evap_Lph = m_dot_working_kgs * 0.015d0 * 3600.0d0

    ! Fan power and super-high COP
    fan_power_kW = (V_dot_product_m3h * 180.0d0) / (3600.0d0 * 0.60d0 * 1000.0d0) + 0.15d0
    COP_energy = Q_cooling_kW / max(0.1d0, fan_power_kW)

    if (T_db_out_C < T_wb_in_C) then
        cooling_status = 'DEW-POINT COOLING (SUB-WET-BULB SUCCESS)'
    else
        cooling_status = 'CONVENTIONAL EVAPORATIVE LIMIT'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — M-CYCLE DEW-POINT EVAPORATIVE COOLER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Ambient Dry-Bulb / RH        = ', T_db_in_C, ' deg C / ', RH_in_pct, ' %'
    write(*,'(A,F10.1,A,F10.1,A)') 'Wet-Bulb (Twb) / Dew-Pt (Tdp)= ', T_wb_in_C, ' deg C / ', T_dp_in_C, ' deg C'
    write(*,'(A,F10.1,A,A)')       'Wet-Bulb Effectiveness       = ', eps_wb_pct, ' % (>100% M-Cycle) | ', trim(cooling_status)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.1,A)')  'SUPPLY PRODUCT OUTLET TEMP   = ', T_db_out_C, ' deg C (Cooled WITHOUT Humidity Addition)'
    write(*,'(A,F10.2,A)')  'SENSIBLE COOLING CAPACITY    = ', Q_cooling_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Evaporative Water Consumption= ', water_evap_Lph, ' Liters/hour'
    write(*,'(A,F10.2,A)')  'Fan Electrical Power Demand  = ', fan_power_kW, ' kW'
    write(*,'(A,F10.1,A)')  'SYSTEM ENERGY COP (EER)      = ', COP_energy, ' (Ultra-Efficient)'
    write(*,'(A)') '============================================================'

end program indirect_evaporative_cooler_mcycle