program gas_absorption_reactive
    implicit none
    integer :: model, i, n_points, iostat_val
    double precision :: kL, DA, k1, k2, CB, CAi, CAb, H, pA, Ei, area
    double precision :: kp, Ha, E, Nphys, Nreact, ratio, totalAbs, driving
    double precision :: Ha_i, E_i, Nphys_i, Nreact_i, ratio_i
    character(len=90) :: modelName, regime

    read(*,*,iostat=iostat_val) model
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid model input.'
        stop
    end if
    read(*,*,iostat=iostat_val) kL
    read(*,*,iostat=iostat_val) DA
    read(*,*,iostat=iostat_val) k1
    read(*,*,iostat=iostat_val) k2
    read(*,*,iostat=iostat_val) CB
    read(*,*,iostat=iostat_val) CAi
    read(*,*,iostat=iostat_val) CAb
    read(*,*,iostat=iostat_val) H
    read(*,*,iostat=iostat_val) pA
    read(*,*,iostat=iostat_val) Ei
    read(*,*,iostat=iostat_val) area
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all reactive absorption inputs.'
        stop
    end if
    if (kL <= 0.0d0 .or. DA <= 0.0d0 .or. Ei < 1.0d0) then
        write(*,*) 'ERROR: kL, DA must be positive and Ei must be at least 1.'
        stop
    end if
    if (CAi <= 0.0d0 .and. H > 0.0d0 .and. pA > 0.0d0) CAi = H*pA
    if (CAi < CAb) then
        write(*,*) 'ERROR: Interfacial concentration must be greater than or equal to bulk concentration.'
        stop
    end if
    if (area < 0.0d0) area = 0.0d0

    if (model == 1) then
        modelName = 'Second-order A+B with pseudo-first-order approximation'
        if (k2 < 0.0d0 .or. CB < 0.0d0) then
            write(*,*) 'ERROR: k2 and CB must be non-negative.'
            stop
        end if
        kp = k2*CB
    else if (model == 2) then
        modelName = 'Direct pseudo-first-order reaction in liquid film'
        if (k1 < 0.0d0) then
            write(*,*) 'ERROR: k1 must be non-negative.'
            stop
        end if
        kp = k1
    else
        write(*,*) 'ERROR: Model must be 1 second-order pseudo-first-order or 2 first-order.'
        stop
    end if

    driving = CAi - CAb
    Ha = sqrt(kp*DA)/kL
    E = enhancement(Ha, Ei)
    Nphys = kL*driving
    Nreact = E*Nphys
    if (Nphys > 0.0d0) then
        ratio = Nreact/Nphys
    else
        ratio = 0.0d0
    end if
    totalAbs = Nreact*area

    if (Ha < 0.3d0) then
        regime = 'Slow reaction; physical absorption nearly controls'
    else if (Ha < 3.0d0) then
        regime = 'Moderate reaction enhancement in liquid film'
    else if (E >= 0.95d0*Ei) then
        regime = 'Near instantaneous-reaction enhancement limit'
    else
        regime = 'Fast reaction; strong chemical enhancement'
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   GAS ABSORPTION - REACTIVE SYSTEMS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Model                    = ', trim(modelName)
    write(*,'(A,ES12.4,A)') '  Liquid Mass Transfer Coeff (kL) = ', kL, ' m/s'
    write(*,'(A,ES12.4,A)') '  Liquid Diffusivity D_A   = ', DA, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Interfacial C_Ai         = ', CAi, ' mol/m3'
    write(*,'(A,ES12.4,A)') '  Bulk C_Ab                = ', CAb, ' mol/m3'
    write(*,'(A,ES12.4)')   '  Pseudo-first-order kprime= ', kp
    write(*,'(A,ES12.4)')   '  Instantaneous Limit (Ei) = ', Ei
    write(*,*)
    write(*,'(A)') '--- REACTIVE ABSORPTION RESULTS -----------------------------'
    write(*,'(A,ES12.4)')   '  Hatta Number (Ha)        = ', Ha
    write(*,'(A,ES12.4)')   '  Enhancement Factor (E)   = ', E
    write(*,'(A,ES12.4,A)') '  Physical Absorption Flux = ', Nphys, ' mol/m2.s'
    write(*,'(A,ES12.4,A)') '  Reactive Absorption Flux = ', Nreact, ' mol/m2.s'
    write(*,'(A,ES12.4)')   '  Flux Ratio               = ', ratio
    write(*,'(A,ES12.4,A)') '  Total Absorption Rate    = ', totalAbs, ' mol/s'
    write(*,'(A,A)')        '  Regime Assessment        = ', trim(regime)
    write(*,*)
    write(*,'(A)') '--- ENHANCEMENT PROFILE -------------------------------------'
    write(*,'(A)') '  Ha           E             N_phys        N_react       ratio'
    write(*,'(A)') '  ----------------------------------------------------------------'
    n_points = 80
    do i=0,n_points-1
        Ha_i = 0.01d0*(1000.0d0)**(dble(i)/dble(n_points-1))
        E_i = enhancement(Ha_i, Ei)
        Nphys_i = Nphys
        Nreact_i = E_i*Nphys_i
        if (Nphys_i > 0.0d0) then
            ratio_i = Nreact_i/Nphys_i
        else
            ratio_i = 0.0d0
        end if
        write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') Ha_i, E_i, Nphys_i, Nreact_i, ratio_i
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Hatta number: Ha = sqrt(kprime*D_A)/kL.'
    write(*,'(A)') '  Pseudo-first-order approximation for A+B: kprime = k2*C_B.'
    write(*,'(A)') '  Enhancement factor approximation: E = Ha/tanh(Ha), capped by Ei.'
    write(*,'(A)') '  Reactive flux: N_A,react = E*kL*(C_Ai-C_Ab).'

contains
    double precision function enhancement(Ha_in, Ei_in)
        implicit none
        double precision, intent(in) :: Ha_in, Ei_in
        double precision :: val
        if (Ha_in < 1.0d-8) then
            val = 1.0d0
        else if (Ha_in > 50.0d0) then
            val = Ha_in
        else
            val = Ha_in/tanh(Ha_in)
        end if
        if (val > Ei_in) val = Ei_in
        if (val < 1.0d0) val = 1.0d0
        enhancement = val
    end function enhancement
end program gas_absorption_reactive
