program absorption_packed_bed
    implicit none
    integer :: i, n_points, iostat_val
    double precision :: y_in, y_out_t, x_in, m, G, L, KGa, area, packingFactor
    double precision :: LG, Aabs, x_out, y_star_bottom, y_star_top, removal
    double precision :: NTU, HTU, Z, Nstage, Gflux, y, x, ystar, drive, ntu_cum, dy, y_hi, y_lo
    double precision :: minL, pinch_top, pinch_bottom

    read(*,*,iostat=iostat_val) y_in
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inlet gas composition input.'
        stop
    end if
    read(*,*,iostat=iostat_val) y_out_t
    read(*,*,iostat=iostat_val) x_in
    read(*,*,iostat=iostat_val) m
    read(*,*,iostat=iostat_val) G
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) KGa
    read(*,*,iostat=iostat_val) area
    read(*,*,iostat=iostat_val) packingFactor
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all absorber inputs.'
        stop
    end if

    if (y_in <= 0.0d0 .or. y_in >= 1.0d0 .or. y_out_t < 0.0d0 .or. y_out_t >= y_in) then
        write(*,*) 'ERROR: Require 0 <= y_out < y_in < 1.'
        stop
    end if
    if (x_in < 0.0d0 .or. x_in >= 1.0d0 .or. m <= 0.0d0) then
        write(*,*) 'ERROR: Liquid inlet composition and Henry slope are invalid.'
        stop
    end if
    if (G <= 0.0d0 .or. L <= 0.0d0 .or. KGa <= 0.0d0 .or. area <= 0.0d0) then
        write(*,*) 'ERROR: G, L, KGa, and column area must be positive.'
        stop
    end if
    if (packingFactor <= 0.0d0) packingFactor = 1.0d0

    LG = L/G
    Aabs = L/(m*G)
    x_out = x_in + G/L * (y_in - y_out_t)
    y_star_top = m*x_in
    y_star_bottom = m*x_out
    pinch_top = y_out_t - y_star_top
    pinch_bottom = y_in - y_star_bottom
    if (pinch_top <= 0.0d0 .or. pinch_bottom <= 0.0d0) then
        write(*,*) 'ERROR: Pinch detected. Increase L/G or reduce target removal.'
        stop
    end if

    removal = (y_in - y_out_t)/y_in * 100.0d0
    NTU = integrate_ntu(y_in, y_out_t, x_in, m, G, L)
    Gflux = G/area
    HTU = Gflux / KGa * packingFactor
    Z = NTU * HTU

    if (abs(Aabs - 1.0d0) < 1.0d-8) then
        Nstage = (y_in - y_out_t) / max(1.0d-12, (y_out_t - y_star_top))
    else
        ! Practical Kremser-style estimate using end-point driving force ratio.
        Nstage = log((y_in - y_star_bottom)/(y_out_t - y_star_top)) / log(Aabs)
        if (Nstage < 0.0d0) Nstage = abs(Nstage)
    end if
    if (Nstage < 0.0d0 .or. Nstage /= Nstage) Nstage = NTU

    minL = m*G*(y_in - y_out_t)/max(1.0d-12, y_in - m*x_in)

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ABSORPTION COLUMN DESIGN - PACKED BED ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4)')   '  Gas Inlet y_in           = ', y_in
    write(*,'(A,ES12.4)')   '  Outlet Gas y_out         = ', y_out_t
    write(*,'(A,ES12.4)')   '  Lean Liquid x_in         = ', x_in
    write(*,'(A,ES12.4)')   '  Henry Slope (m)          = ', m
    write(*,'(A,ES12.4,A)') '  Gas Molar Flow (G)       = ', G, ' mol/s'
    write(*,'(A,ES12.4,A)') '  Liquid Molar Flow (L)    = ', L, ' mol/s'
    write(*,'(A,ES12.4,A)') '  Volumetric K_Ga          = ', KGa, ' mol/m3.s'
    write(*,'(A,ES12.4,A)') '  Column Area              = ', area, ' m2'
    write(*,*)
    write(*,'(A)') '--- EQUILIBRIUM / OPERATING LINE ---------------------------'
    write(*,'(A,ES12.4)')   '  L/G Ratio                = ', LG
    write(*,'(A,ES12.4)')   '  Absorption Factor        = ', Aabs
    write(*,'(A,ES12.4)')   '  Outlet Liquid x_out      = ', x_out
    write(*,'(A,ES12.4)')   '  y*_top                   = ', y_star_top
    write(*,'(A,ES12.4)')   '  y*_bottom                = ', y_star_bottom
    write(*,'(A,ES12.4,A)') '  Minimum Liquid Flow      = ', minL, ' mol/s'
    write(*,*)
    write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
    write(*,'(A,ES12.4)')   '  Theoretical Stages       = ', Nstage
    write(*,'(A,ES12.4)')   '  Overall Gas NTU          = ', NTU
    write(*,'(A,ES12.4,A)') '  Overall Gas HTU          = ', HTU, ' m'
    write(*,'(A,ES12.4,A)') '  Column Height            = ', Z, ' m'
    write(*,'(A,ES12.4,A)') '  Removal Efficiency       = ', removal, ' percent'
    write(*,*)

    write(*,'(A)') '--- COLUMN PROFILE ------------------------------------------'
    write(*,'(A)') '  z/H          y_gas         x_liq         y_star        driving       NTU_cum'
    write(*,'(A)') '  ---------------------------------------------------------------------------'
    n_points = 60
    ntu_cum = 0.0d0
    do i = 0, n_points-1
        y = y_out_t + (y_in - y_out_t)*dble(i)/dble(n_points-1)
        x = x_in + G/L*(y - y_out_t)
        ystar = m*x
        drive = y - ystar
        if (i > 0) then
            y_hi = y_out_t + (y_in - y_out_t)*dble(i)/dble(n_points-1)
            y_lo = y_out_t + (y_in - y_out_t)*dble(i-1)/dble(n_points-1)
            dy = y_hi - y_lo
            ntu_cum = ntu_cum + dy / max(1.0d-12, drive)
        end if
        write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            dble(i)/dble(n_points-1), y, x, ystar, drive, ntu_cum
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Henry law: y* = m x.'
    write(*,'(A)') '  Operating line: x = x_in + (G/L)(y - y_out).'
    write(*,'(A)') '  NTU_OG = integral dy/(y - y*) from y_out to y_in.'
    write(*,'(A)') '  HTU_OG = G_molar_flux / K_Ga; Z = NTU_OG * HTU_OG.'
    write(*,'(A)') '  Kremser-style stage estimate based on absorption factor A = L/(mG).'

contains
    double precision function integrate_ntu(yin, yout, xin, mh, Gf, Lf)
        implicit none
        double precision, intent(in) :: yin, yout, xin, mh, Gf, Lf
        integer :: j, N
        double precision :: y, x, yst, f, sumv, dy, w
        N = 2000
        dy = (yin - yout)/dble(N)
        sumv = 0.0d0
        do j = 0, N
            y = yout + dy*dble(j)
            x = xin + Gf/Lf*(y - yout)
            yst = mh*x
            f = 1.0d0/max(1.0d-12, y - yst)
            if (j == 0 .or. j == N) then
                w = 0.5d0
            else
                w = 1.0d0
            end if
            sumv = sumv + w*f
        end do
        integrate_ntu = sumv*dy
    end function integrate_ntu
end program absorption_packed_bed
