program adsorption_breakthrough_bed
    implicit none
    integer :: iostat_val
    double precision :: Z_bed_m, D_col_m, rho_bulk_gL, Q_Lh, C0_mgL, q0_mgg, kTh_mLmgmin
    double precision :: kTh_Lmg_h, A_col_m2, V_bed_L, M_ads_g, u0_mh
    double precision :: tb_h, ts_h, V_treated_L, MTZ_m, Z0_min_m, t50_h
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) Z_bed_m        ! Bed Height [m] (e.g. 1.20)
    read(*,*,iostat=iostat_val) D_col_m        ! Column Diameter [m] (e.g. 0.30)
    read(*,*,iostat=iostat_val) rho_bulk_gL    ! Adsorbent Bulk Density [g/L] (e.g. 500.0)
    read(*,*,iostat=iostat_val) Q_Lh           ! Liquid Flow Rate [L/h] (e.g. 200.0)
    read(*,*,iostat=iostat_val) C0_mgL         ! Inlet Solute Concentration [mg/L] (e.g. 50.0)
    read(*,*,iostat=iostat_val) q0_mgg         ! Maximum Adsorption Capacity [mg/g] (e.g. 45.0)
    read(*,*,iostat=iostat_val) kTh_mLmgmin    ! Thomas Kinetic Constant [mL/(mg.min)] (e.g. 0.15)

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

    if (Z_bed_m <= 0.0d0 .or. D_col_m <= 0.0d0 .or. Q_Lh <= 0.0d0 .or. C0_mgL <= 0.0d0) then
        write(*,*) 'ERROR: Bed geometry, flow rate, and concentration must be positive.'
        stop
    end if

    ! Unit conversion
    kTh_Lmg_h = kTh_mLmgmin * 1.0d-3 * 60.0d0 ! L/(mg.h)

    A_col_m2 = (PI / 4.0d0) * (D_col_m**2)
    V_bed_L = A_col_m2 * Z_bed_m * 1000.0d0
    M_ads_g = V_bed_L * rho_bulk_gL
    u0_mh = (Q_Lh * 1.0d-3) / A_col_m2

    ! 50% Breakthrough midpoint time t50 [hours]
    t50_h = (q0_mgg * M_ads_g) / (C0_mgL * Q_Lh)

    ! Breakthrough time tb at C/C0 = 0.05 (5% leakage)
    tb_h = t50_h - (log(0.95d0 / 0.05d0) / (kTh_Lmg_h * C0_mgL))
    if (tb_h < 0.0d0) tb_h = 0.0d0

    ! Exhaustion / Saturation time ts at C/C0 = 0.95 (95% saturation)
    ts_h = t50_h + (log(0.95d0 / 0.05d0) / (kTh_Lmg_h * C0_mgL))

    ! Total treated effluent volume before breakthrough
    V_treated_L = Q_Lh * tb_h

    ! Mass Transfer Zone length MTZ [m]
    MTZ_m = Z_bed_m * (ts_h - tb_h) / max(0.01d0, ts_h)

    ! Critical minimum bed depth Z0 [m] (Bohart-Adams BDST)
    Z0_min_m = (u0_mh / (kTh_Lmg_h * (q0_mgg * rho_bulk_gL))) * log(1.0d0 / 0.05d0 - 1.0d0)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — FIXED-BED ADSORPTION BREAKTHROUGH'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Bed Height / Diameter     = ', Z_bed_m, ' m / ', D_col_m, ' m'
    write(*,'(A,F10.1,A,F10.1,A)') 'Adsorbent Mass in Bed     = ', M_ads_g / 1000.0d0, ' kg (', V_bed_L, ' L)'
    write(*,'(A,F10.1,A,F10.1,A)') 'Liquid Feed Flow / C0     = ', Q_Lh, ' L/h / ', C0_mgL, ' mg/L'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  '5% Breakthrough Time (tb) = ', tb_h, ' hours'
    write(*,'(A,F10.2,A)')  '95% Exhaustion Time (ts)  = ', ts_h, ' hours'
    write(*,'(A,F10.1,A)')  'Treated Effluent Volume   = ', V_treated_L, ' Liters'
    write(*,'(A,F10.3,A)')  'Mass Transfer Zone (MTZ)  = ', MTZ_m, ' meters'
    write(*,'(A,F10.3,A)')  'Critical Bed Depth (Z0)   = ', Z0_min_m, ' meters (BDST)'
    write(*,'(A,F10.2,A)')  'Bed Service Lifetime      = ', tb_h / 24.0d0, ' days continuous'
    write(*,'(A)') '============================================================'

end program adsorption_breakthrough_bed