program falling_film_absorber
    implicit none
    integer :: iostat_val
    double precision :: Dt_mm, Lt_m, Nt_tubes, ML_kgh, rhoL_kgm3, muL_cP, DAB_m2s
    double precision :: Dt_m, ML_kgs, Gamma_wetted, Re_film, delta_film_m, delta_film_mm
    double precision :: u_surf_ms, tc_s, kL_ms, kL_mh, Area_wetted_m2, Sh_film
    double precision, parameter :: PI = 3.141592653589793d0, g_acc = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Dt_mm         ! Tube Inside Diameter [mm] (e.g. 38.0)
    read(*,*,iostat=iostat_val) Lt_m          ! Tube Length [m] (e.g. 3.0)
    read(*,*,iostat=iostat_val) Nt_tubes      ! Number of Tubes (e.g. 12.0)
    read(*,*,iostat=iostat_val) ML_kgh        ! Total Liquid Flow Rate [kg/h] (e.g. 1500.0)
    read(*,*,iostat=iostat_val) rhoL_kgm3     ! Liquid Density [kg/m3] (e.g. 1000.0)
    read(*,*,iostat=iostat_val) muL_cP        ! Liquid Viscosity [cP] (e.g. 1.0)
    read(*,*,iostat=iostat_val) DAB_m2s       ! Solute Diffusivity in Liquid [m2/s] (e.g. 1.8e-9)

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

    if (Dt_mm <= 0.0d0 .or. Lt_m <= 0.0d0 .or. Nt_tubes <= 0.0d0 .or. ML_kgh <= 0.0d0) then
        write(*,*) 'ERROR: Tube geometry and liquid flow rate must be positive.'
        stop
    end if

    Dt_m = Dt_mm * 1.0d-3
    ML_kgs = (ML_kgh / 3600.0d0) / Nt_tubes ! per tube

    ! Wetted perimeter loading Gamma [kg/(m.s)]
    Gamma_wetted = ML_kgs / (PI * Dt_m)

    ! Film Reynolds number Re_f = 4 * Gamma / mu
    Re_film = (4.0d0 * Gamma_wetted) / (muL_cP * 1.0d-3)

    ! Nusselt film thickness delta [m]
    delta_film_m = ((3.0d0 * (muL_cP * 1.0d-3) * Gamma_wetted) / ((rhoL_kgm3**2) * g_acc))**(1.0d0/3.0d0)
    delta_film_mm = delta_film_m * 1000.0d0

    ! Surface velocity u_s [m/s]
    u_surf_ms = (rhoL_kgm3 * g_acc * (delta_film_m**2)) / (2.0d0 * (muL_cP * 1.0d-3))

    ! Liquid contact time t_c [s]
    tc_s = Lt_m / max(0.01d0, u_surf_ms)

    ! Higbie Penetration mass transfer coefficient k_L [m/s]
    kL_ms = 2.0d0 * sqrt(DAB_m2s / (PI * max(0.001d0, tc_s)))
    kL_mh = kL_ms * 3600.0d0

    Area_wetted_m2 = Nt_tubes * (PI * Dt_m * Lt_m)
    Sh_film = (kL_ms * delta_film_m) / DAB_m2s

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — FALLING FILM ABSORBER (NUSSELT & HIGBIE)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.2,A)') 'Tubes: (Nt x Dt x Lt)     = ', Nt_tubes, ' tubes x ', Dt_mm, ' mm x ', Lt_m, ' m'
    write(*,'(A,F10.2,A,F10.2,A)') 'Total Wetted Area / Liquid= ', Area_wetted_m2, ' m2 / ', ML_kgh, ' kg/h'
    write(*,'(A,F10.1,A,F10.4,A)') 'Film Reynolds / Gamma     = ', Re_film, ' / ', Gamma_wetted, ' kg/(m.s)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.3,A)')  'Nusselt Film Thickness    = ', delta_film_mm, ' mm'
    write(*,'(A,F10.3,A)')  'Surface Film Velocity (us)= ', u_surf_ms, ' m/s'
    write(*,'(A,F10.2,A)')  'Liquid Contact Time (tc)  = ', tc_s, ' seconds'
    write(*,'(A,F10.2,A,F10.2,A)') 'Mass Transfer Coeff (kL)  = ', kL_mh, ' m/h (', kL_ms*1.0d5, ' x 10^-5 m/s)'
    write(*,'(A,F10.2)')    'Film Sherwood Number (Sh) = ', Sh_film
    write(*,'(A)') '============================================================'

end program falling_film_absorber
