program merkel_cooling_tower_transfer
    implicit none
    integer :: iostat_val, i
    double precision :: Twin_C, Twout_C, Tdb_C, Twb_C, L_flow_m3h, LG_ratio
    double precision :: range_C, approach_C, tower_eff_pct, KaV_L, evap_m3h
    double precision :: ha_in_kJkg, T_points(4), h_sat(4), h_air(4), dh(4)
    double precision, parameter :: c_w = 4.186d0 ! kJ/(kg.K)

    ! Read inputs
    read(*,*,iostat=iostat_val) Twin_C         ! Hot Water Inlet Temp [deg C] (e.g. 37.0)
    read(*,*,iostat=iostat_val) Twout_C        ! Cold Water Outlet Temp [deg C] (e.g. 29.0)
    read(*,*,iostat=iostat_val) Tdb_C          ! Air Dry Bulb Temp [deg C] (e.g. 32.0)
    read(*,*,iostat=iostat_val) Twb_C          ! Air Wet Bulb Temp [deg C] (e.g. 24.0)
    read(*,*,iostat=iostat_val) L_flow_m3h     ! Water Flow Rate [m3/h] (e.g. 500.0)
    read(*,*,iostat=iostat_val) LG_ratio       ! Liquid-to-Gas Ratio L/G (e.g. 1.25)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for cooling tower Merkel calculation.'
        stop
    end if

    if (Twin_C <= Twout_C .or. Twout_C <= Twb_C .or. LG_ratio <= 0.0d0) then
        write(*,*) 'ERROR: Must satisfy Twin > Twout > Twb and LG_ratio > 0.'
        stop
    end if

    range_C = Twin_C - Twout_C
    approach_C = Twout_C - Twb_C
    tower_eff_pct = (range_C / (Twin_C - Twb_C)) * 100.0d0

    ! Inlet air enthalpy at Twb (approx saturated air enthalpy at Twb)
    ha_in_kJkg = 1.006d0 * Twb_C + (0.62198d0 * (0.61078d0 * exp(17.27d0 * Twb_C / (Twb_C + 237.3d0))) / &
                 (101.325d0 - 0.61078d0 * exp(17.27d0 * Twb_C / (Twb_C + 237.3d0)))) * (2501.0d0 + 1.86d0 * Twb_C)

    ! Chebyshev 4-point quadrature temperatures (CTI Bulletin P-130)
    T_points(1) = Twout_C + 0.1d0 * range_C
    T_points(2) = Twout_C + 0.4d0 * range_C
    T_points(3) = Twout_C + 0.6d0 * range_C
    T_points(4) = Twout_C + 0.9d0 * range_C

    KaV_L = 0.0d0
    do i = 1, 4
        ! Saturated air enthalpy at water temp T_points(i)
        h_sat(i) = 1.006d0 * T_points(i) + (0.62198d0 * (0.61078d0 * exp(17.27d0 * T_points(i) / (T_points(i) + 237.3d0))) / &
                   (101.325d0 - 0.61078d0 * exp(17.27d0 * T_points(i) / (T_points(i) + 237.3d0)))) * (2501.0d0 + 1.86d0 * T_points(i))
        ! Air enthalpy in tower
        h_air(i) = ha_in_kJkg + LG_ratio * c_w * (T_points(i) - Twout_C)
        dh(i) = max(0.5d0, h_sat(i) - h_air(i))
        KaV_L = KaV_L + (1.0d0 / dh(i))
    end do
    KaV_L = (c_w * range_C / 4.0d0) * KaV_L

    ! Evaporation loss rate [m3/h]
    evap_m3h = 0.0018d0 * L_flow_m3h * range_C

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — COOLING TOWER MERKEL INTEGRAL ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Hot Water In / Basin Out  = ', Twin_C, ' C / ', Twout_C, ' C'
    write(*,'(A,F10.1,A,F10.1,A)') 'Dry Bulb / Wet Bulb Temp  = ', Tdb_C, ' C / ', Twb_C, ' C'
    write(*,'(A,F10.2,A,F10.2,A)') 'Cooling Range / Approach  = ', range_C, ' C / ', approach_C, ' C'
    write(*,'(A,F10.2,A,F10.1,A)') 'Liquid-to-Gas Ratio L/G   = ', LG_ratio, ' (Water: ', L_flow_m3h, ' m3/h)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.3)')    'MERKEL NUMBER (KaV / L)   = ', KaV_L
    write(*,'(A,F10.2,A)')  'Thermal Effectiveness     = ', tower_eff_pct, ' %'
    write(*,'(A,F10.2,A)')  'Evaporation Water Loss    = ', evap_m3h, ' m3/h'
    write(*,'(A,F10.2,A)')  'Air Flow Rate Required    = ', (L_flow_m3h * 998.0d0 / LG_ratio) / 1.18d0, ' m3/h'
    write(*,'(A)') '============================================================'

end program merkel_cooling_tower_transfer