Falling Film Evaporator & Desalination Sizer

program falling_film_evaporator_desalination
    implicit none
    integer :: iostat_val, feed_type, num_tubes
    double precision :: m_dot_feed_kgh, T_feed_C, T_steam_C, T_sat_evap_C, D_tube_mm, L_tube_m
    double precision :: rho_L, cp_L, mu_L, k_L, hfg_kJkg, sigma_Nm, h_steam, t_wall_mm, k_wall
    double precision :: m_dot_feed_kgs, D_i_m, Gamma_wetting, Re_film, Pr_L, h_film, U_overall
    double precision :: A_total_m2, delta_T_eff, Q_duty_kW, m_dot_evap_kgh, steam_economy
    character(len=32) :: wetting_status, film_regime
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: G_ACC = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) feed_type       ! 1=Seawater Desalination (MED), 2=Milk / Dairy Concentration, 3=Black Liquor, 4=Sugar Juice
    read(*,*,iostat=iostat_val) m_dot_feed_kgh  ! Feed Mass Flow Rate [kg/h] (e.g. 5000.0)
    read(*,*,iostat=iostat_val) T_feed_C        ! Feed Temperature [deg C] (e.g. 68.0)
    read(*,*,iostat=iostat_val) T_steam_C       ! Heating Steam Sat Temp [deg C] (e.g. 75.0)
    read(*,*,iostat=iostat_val) T_sat_evap_C    ! Evaporation Chamber Sat Temp [deg C] (e.g. 65.0)
    read(*,*,iostat=iostat_val) num_tubes       ! Number of Vertical Tubes (e.g. 120)
    read(*,*,iostat=iostat_val) D_tube_mm       ! Tube Inner Diameter [mm] (e.g. 38.0)
    read(*,*,iostat=iostat_val) L_tube_m        ! Tube Length [m] (e.g. 4.50)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for falling film evaporator calculation.'
        stop
    end if

    if (m_dot_feed_kgh <= 0.0d0 .or. num_tubes <= 0 .or. D_tube_mm <= 0.0d0 .or. L_tube_m <= 0.0d0) then
        write(*,*) 'ERROR: Feed rate, tubes count, and dimensions must be positive.'
        stop
    end if

    ! Feed liquid properties
    if (feed_type == 1) then ! Seawater MED Desalination
        rho_L = 1015.0d0; cp_L = 4050.0d0; mu_L = 0.00045d0; k_L = 0.64d0; hfg_kJkg = 2345.0d0
    else if (feed_type == 2) then ! Milk / Dairy
        rho_L = 1030.0d0; cp_L = 3900.0d0; mu_L = 0.00120d0; k_L = 0.58d0; hfg_kJkg = 2345.0d0
    else if (feed_type == 3) then ! Black liquor
        rho_L = 1180.0d0; cp_L = 3200.0d0; mu_L = 0.00850d0; k_L = 0.42d0; hfg_kJkg = 2345.0d0
    else ! Sugar juice
        rho_L = 1060.0d0; cp_L = 3700.0d0; mu_L = 0.00180d0; k_L = 0.55d0; hfg_kJkg = 2345.0d0
    end if

    m_dot_feed_kgs = m_dot_feed_kgh / 3600.0d0
    D_i_m = D_tube_mm / 1000.0d0
    t_wall_mm = 1.65d0    ! Titanium / 316L SS tube
    k_wall = 22.0d0
    h_steam = 9500.0d0    ! Condensing steam film coefficient

    A_total_m2 = dble(num_tubes) * PI * D_i_m * L_tube_m

    ! Peripheral Wetting Rate Gamma [kg/(m.s)]
    Gamma_wetting = m_dot_feed_kgs / (dble(num_tubes) * PI * D_i_m)
    Re_film = 4.0d0 * Gamma_wetting / mu_L
    Pr_L = (cp_L * mu_L) / k_L

    if (Gamma_wetting < 0.04d0) then
        wetting_status = 'WARNING: RISK OF DRYOUT'
    else
        wetting_status = 'SAFE (FULL WETTING FILM)'
    end if

    ! Chun & Seban Falling Film Evaporation Correlation
    if (Re_film < 1600.0d0) then
        film_regime = 'LAMINAR-WAVY FILM'
        h_film = 0.822d0 * (((k_L**3) * (rho_L**2) * G_ACC / (mu_L**2))**(1.0d0/3.0d0)) * (Re_film**(-0.22d0))
    else
        film_regime = 'TURBULENT FILM'
        h_film = 0.0038d0 * (((k_L**3) * (rho_L**2) * G_ACC / (mu_L**2))**(1.0d0/3.0d0)) * (Re_film**0.40d0) * (Pr_L**0.65d0)
    end if

    ! Overall Heat Transfer Coefficient U [W/(m2.K)]
    U_overall = 1.0d0 / ((1.0d0 / h_film) + ((t_wall_mm/1000.0d0) / k_wall) + (1.0d0 / h_steam))

    delta_T_eff = max(0.5d0, T_steam_C - T_sat_evap_C)
    Q_duty_kW = (U_overall * A_total_m2 * delta_T_eff) / 1000.0d0

    m_dot_evap_kgh = (Q_duty_kW / hfg_kJkg) * 3600.0d0
    steam_economy = m_dot_evap_kgh / max(1.0d0, (Q_duty_kW / 2250.0d0) * 3600.0d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — FALLING FILM EVAPORATOR SIZING'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.3,A,A)')       'Peripheral Wetting Rate (Γ)  = ', Gamma_wetting, ' kg/(m.s) | ', trim(wetting_status)
    write(*,'(A,F10.0,A,A)')       'Film Reynolds Number (Re_f)  = ', Re_film, ' | ', trim(film_regime)
    write(*,'(A,F10.2,A,I6,A)')    'Total Evaporator Area        = ', A_total_m2, ' m2 (', num_tubes, ' tubes)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'THERMAL HEAT DUTY (Q)        = ', Q_duty_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Falling Film Coeff (h_film)  = ', h_film, ' W/(m2.K)'
    write(*,'(A,F10.1,A)')  'Overall Heat Transfer Coeff U= ', U_overall, ' W/(m2.K)'
    write(*,'(A,F10.1,A)')  'EVAPORATED VAPOR PRODUCTION  = ', m_dot_evap_kgh, ' kg/h'
    write(*,'(A,F10.2)')    'Single-Effect Steam Economy  = ', steam_economy
    write(*,'(A)') '============================================================'

end program falling_film_evaporator_desalination