program evaporation_humidification
    implicit none
    integer :: geom_type, i, n_points, iostat_val
    double precision :: L, W, Dhyd, u, Tdb, Ts, P_kPa, RH, Dab, rho, mu, kair, cp, MWA, activity, flowArea
    double precision :: area, Lc, nu, alpha, Le, Re, Sc, Sh, hm, P, Rgas, Ctot
    double precision :: psat_s, psat_db, pA_s, pA_inf, y_s, y_inf, omega_in, omega_s, omega_out, RH_out
    double precision :: NA, mass_flux, evap_rate, mdot_air, T_wb, wb_dep, hfg, h_in, h_out_sat
    double precision :: RHp, pA_in_p, yinfp, omegap, NAp, ratep, Twbp, depp, RHoutp, omegaoutp
    character(len=80) :: geom_name, corr_name

    read(*,*,iostat=iostat_val) geom_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid geometry type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) W
    read(*,*,iostat=iostat_val) Dhyd
    read(*,*,iostat=iostat_val) u
    read(*,*,iostat=iostat_val) Tdb
    read(*,*,iostat=iostat_val) Ts
    read(*,*,iostat=iostat_val) P_kPa
    read(*,*,iostat=iostat_val) RH
    read(*,*,iostat=iostat_val) Dab
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) kair
    read(*,*,iostat=iostat_val) cp
    read(*,*,iostat=iostat_val) MWA
    read(*,*,iostat=iostat_val) activity
    read(*,*,iostat=iostat_val) flowArea
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all evaporation/humidification parameters.'
        stop
    end if

    if (L <= 0.0d0 .or. W <= 0.0d0 .or. Dhyd <= 0.0d0) then
        write(*,*) 'ERROR: Geometry dimensions must be positive.'
        stop
    end if
    if (P_kPa <= 0.0d0 .or. Dab <= 0.0d0 .or. rho <= 0.0d0 .or. mu <= 0.0d0) then
        write(*,*) 'ERROR: Pressure and fluid properties must be positive.'
        stop
    end if
    if (RH < 0.0d0 .or. RH > 100.0d0) then
        write(*,*) 'ERROR: Relative humidity must be between 0 and 100 percent.'
        stop
    end if
    if (activity < 0.0d0) then
        write(*,*) 'ERROR: Activity/Raoult factor cannot be negative.'
        stop
    end if

    if (geom_type == 1) then
        geom_name = 'Open Surface / Flat Plate'
        Lc = L
        area = L * W
    else
        geom_name = 'Humidifier Duct / Wetted Pad'
        Lc = Dhyd
        area = L * W
    end if

    P = P_kPa * 1000.0d0
    Rgas = 8.314462618d0
    Ctot = P / (Rgas * (Tdb + 273.15d0))
    nu = mu / rho
    alpha = kair / (rho * cp)
    Le = alpha / Dab
    Re = rho * u * Lc / mu
    Sc = mu / (rho * Dab)
    call calc_sh(geom_type, Re, Sc, Sh, corr_name)
    hm = Sh * Dab / Lc

    psat_s = saturation_pressure_water(Ts)
    psat_db = saturation_pressure_water(Tdb)
    pA_s = min(activity * psat_s, 0.995d0 * P)
    pA_inf = min((RH/100.0d0) * psat_db, 0.995d0 * P)
    y_s = pA_s / P
    y_inf = pA_inf / P
    omega_in = humidity_ratio(P, pA_inf)
    omega_s = humidity_ratio(P, pA_s)

    if (y_s > y_inf) then
        NA = hm * Ctot * log((1.0d0 - y_inf)/(1.0d0 - y_s))
    else
        NA = hm * Ctot * (y_s - y_inf)
    end if
    mass_flux = NA * MWA / 1000.0d0
    evap_rate = mass_flux * area

    if (flowArea <= 0.0d0) flowArea = area
    mdot_air = rho * u * flowArea
    if (mdot_air > 0.0d0) then
        omega_out = omega_in + evap_rate / mdot_air
    else
        omega_out = omega_in
    end if
    RH_out = relative_humidity_from_omega(P, Tdb, omega_out) * 100.0d0
    if (RH_out > 100.0d0) RH_out = 100.0d0

    T_wb = wet_bulb_temperature(Tdb, P, omega_in)
    wb_dep = Tdb - T_wb
    hfg = latent_heat_water(Ts)
    h_in = moist_air_enthalpy(Tdb, omega_in)
    h_out_sat = moist_air_enthalpy(T_wb, humidity_ratio(P, saturation_pressure_water(T_wb)))

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   EVAPORATION AND HUMIDIFICATION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Geometry Name             = ', trim(geom_name)
    write(*,'(A,A)')        '  Correlation               = ', trim(corr_name)
    write(*,'(A,ES12.4,A)') '  Transfer Area             = ', area, ' m2'
    write(*,'(A,ES12.4,A)') '  Characteristic Length     = ', Lc, ' m'
    write(*,*)
    write(*,'(A)') '--- AIR / VAPOR STATE ---------------------------------------'
    write(*,'(A,F12.4,A)')  '  Dry Bulb Temperature      = ', Tdb, ' deg-C'
    write(*,'(A,F12.4,A)')  '  Surface Temperature       = ', Ts, ' deg-C'
    write(*,'(A,F12.4,A)')  '  Total Pressure            = ', P_kPa, ' kPa'
    write(*,'(A,F12.4,A)')  '  Inlet Relative Humidity   = ', RH, ' percent'
    write(*,'(A,ES12.4,A)') '  Surface Vapor Pressure    = ', pA_s, ' Pa'
    write(*,'(A,ES12.4,A)') '  Bulk Vapor Pressure       = ', pA_inf, ' Pa'
    write(*,'(A,ES12.4)')   '  Surface Mole Fraction     = ', y_s
    write(*,'(A,ES12.4)')   '  Bulk Mole Fraction        = ', y_inf
    write(*,'(A,ES12.4)')   '  Inlet Humidity Ratio      = ', omega_in
    write(*,'(A,ES12.4)')   '  Outlet Humidity Ratio     = ', omega_out
    write(*,'(A,F12.4,A)')  '  Outlet Relative Humidity  = ', RH_out, ' percent'
    write(*,*)
    write(*,'(A)') '--- DIMENSIONLESS GROUPS ------------------------------------'
    write(*,'(A,ES12.4)')   '  Reynolds Number (Re)      = ', Re
    write(*,'(A,ES12.4)')   '  Schmidt Number (Sc)       = ', Sc
    write(*,'(A,ES12.4)')   '  Sherwood Number (Sh)      = ', Sh
    write(*,'(A,ES12.4)')   '  Lewis Number (Le)         = ', Le
    write(*,*)
    write(*,'(A)') '--- EVAPORATION / HUMIDIFICATION RESULTS --------------------'
    write(*,'(A,ES12.4,A)') '  Mass Transfer Coeff (h_m) = ', hm, ' m/s'
    write(*,'(A,ES12.4,A)') '  Molar Flux (N_A)          = ', NA, ' mol/m2.s'
    write(*,'(A,ES12.4,A)') '  Mass Flux                 = ', mass_flux, ' kg/m2.s'
    write(*,'(A,ES12.4,A)') '  Evaporation Rate          = ', evap_rate, ' kg/s'
    write(*,'(A,F12.4,A)')  '  Wet Bulb Temperature      = ', T_wb, ' deg-C'
    write(*,'(A,F12.4,A)')  '  Wet Bulb Depression       = ', wb_dep, ' K'
    write(*,'(A,ES12.4,A)') '  Latent Heat h_fg          = ', hfg, ' J/kg'
    write(*,'(A,ES12.4,A)') '  Psychrometric h_in        = ', h_in, ' J/kg dry-air'
    write(*,'(A,ES12.4,A)') '  Psychrometric h_sat_wb    = ', h_out_sat, ' J/kg dry-air'
    write(*,*)

    write(*,'(A)') '--- PARAMETRIC HUMIDIFICATION PROFILE -----------------------'
    write(*,'(A)') '  RH_in[%]    omega_in      evap_rate[kg/s]  T_wb[C]   depression[K]  RH_out[%]'
    write(*,'(A)') '  ---------------------------------------------------------------------------'
    n_points = 41
    do i = 0, n_points-1
        RHp = 5.0d0 + 90.0d0 * dble(i) / dble(n_points-1)
        pA_in_p = min((RHp/100.0d0) * psat_db, 0.995d0 * P)
        yinfp = pA_in_p / P
        omegap = humidity_ratio(P, pA_in_p)
        if (y_s > yinfp) then
            NAp = hm * Ctot * log((1.0d0 - yinfp)/(1.0d0 - y_s))
        else
            NAp = hm * Ctot * (y_s - yinfp)
        end if
        ratep = NAp * MWA / 1000.0d0 * area
        Twbp = wet_bulb_temperature(Tdb, P, omegap)
        depp = Tdb - Twbp
        if (mdot_air > 0.0d0) then
            omegaoutp = omegap + ratep / mdot_air
        else
            omegaoutp = omegap
        end if
        RHoutp = min(100.0d0, relative_humidity_from_omega(P, Tdb, omegaoutp) * 100.0d0)
        write(*,'(F9.3,2X,ES12.4,2X,ES14.4,2X,F9.3,2X,F12.4,2X,F10.3)') RHp, omegap, ratep, Twbp, depp, RHoutp
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Raoult law: p_A,s = activity * p_sat(T_s).'
    write(*,'(A)') '  Stefan diffusion: N_A = h_m C ln[(1-y_A,inf)/(1-y_A,s)].'
    write(*,'(A)') '  Lewis number: Le = alpha/D_AB = k/(rho cp D_AB).'
    write(*,'(A)') '  Wet-bulb estimated by moist-air enthalpy equality.'

contains
    double precision function saturation_pressure_water(Tc)
        implicit none
        double precision, intent(in) :: Tc
        ! Buck/Magnus-style equation over liquid water, Pa.
        saturation_pressure_water = 611.21d0 * exp((18.678d0 - Tc/234.5d0) * (Tc/(257.14d0 + Tc)))
    end function saturation_pressure_water

    double precision function humidity_ratio(Pa, pv)
        implicit none
        double precision, intent(in) :: Pa, pv
        humidity_ratio = 0.621945d0 * pv / max(1.0d-9, Pa - pv)
    end function humidity_ratio

    double precision function vapor_pressure_from_omega(Pa, omega)
        implicit none
        double precision, intent(in) :: Pa, omega
        vapor_pressure_from_omega = Pa * omega / (0.621945d0 + omega)
    end function vapor_pressure_from_omega

    double precision function relative_humidity_from_omega(Pa, Tc, omega)
        implicit none
        double precision, intent(in) :: Pa, Tc, omega
        double precision :: pv, ps
        pv = vapor_pressure_from_omega(Pa, omega)
        ps = saturation_pressure_water(Tc)
        relative_humidity_from_omega = pv / ps
    end function relative_humidity_from_omega

    double precision function moist_air_enthalpy(Tc, omega)
        implicit none
        double precision, intent(in) :: Tc, omega
        moist_air_enthalpy = 1006.0d0*Tc + omega*(2501000.0d0 + 1860.0d0*Tc)
    end function moist_air_enthalpy

    double precision function latent_heat_water(Tc)
        implicit none
        double precision, intent(in) :: Tc
        latent_heat_water = (2501.0d0 - 2.361d0*Tc) * 1000.0d0
    end function latent_heat_water

    double precision function wet_bulb_temperature(Tdry, Pa, omega)
        implicit none
        double precision, intent(in) :: Tdry, Pa, omega
        integer :: it
        double precision :: lo, hi, mid, htarget, hmid, ws
        htarget = moist_air_enthalpy(Tdry, omega)
        lo = -50.0d0
        hi = Tdry
        do it = 1, 80
            mid = 0.5d0*(lo+hi)
            ws = humidity_ratio(Pa, saturation_pressure_water(mid))
            hmid = moist_air_enthalpy(mid, ws)
            if (hmid > htarget) then
                hi = mid
            else
                lo = mid
            end if
        end do
        wet_bulb_temperature = 0.5d0*(lo+hi)
    end function wet_bulb_temperature

    subroutine calc_sh(g, Re_in, Sc_in, Sh_out, corr)
        implicit none
        integer, intent(in) :: g
        double precision, intent(in) :: Re_in, Sc_in
        double precision, intent(out) :: Sh_out
        character(len=80), intent(out) :: corr
        if (g == 1) then
            if (Re_in < 5.0d5) then
                corr = 'Flat plate laminar: Sh = 0.664 Re^0.5 Sc^(1/3)'
                Sh_out = 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0)
            else
                corr = 'Flat plate turbulent: Sh = (0.037 Re^0.8 - 871) Sc^(1/3)'
                Sh_out = (0.037d0*max(Re_in,0.0d0)**0.8d0 - 871.0d0) * Sc_in**(1.0d0/3.0d0)
            end if
        else
            if (Re_in < 2300.0d0) then
                corr = 'Duct laminar approximate: Sh = 3.66'
                Sh_out = 3.66d0
            else
                corr = 'Duct turbulent mass analogy: Sh = 0.023 Re^0.8 Sc^(1/3)'
                Sh_out = 0.023d0 * max(Re_in,0.0d0)**0.8d0 * Sc_in**(1.0d0/3.0d0)
            end if
        end if
        if (Sh_out < 0.0d0) Sh_out = 0.0d0
    end subroutine calc_sh
end program evaporation_humidification
