program convective_mass_transfer
    implicit none
    integer :: geom_type, flow_type, i, n_points, iostat_val
    double precision :: L, Dpipe, width, velocity, rho, mu, Dab, CA_bulk, CA_surface, MWA, lengthPipe
    double precision :: Lc, area, nu, Re, Sc, Sh, km, NA, mass_flux, total_mdot, delta_m, dC
    double precision :: x, Rex, Shx, kmx, NAx, filmx, localLc
    character(len=80) :: geom_name, regime_name, corr_name
    double precision, parameter :: PI = 3.1415926535897932384626433832795d0

    read(*,*,iostat=iostat_val) geom_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid geometry type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) flow_type
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) Dpipe
    read(*,*,iostat=iostat_val) width
    read(*,*,iostat=iostat_val) velocity
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) Dab
    read(*,*,iostat=iostat_val) CA_bulk
    read(*,*,iostat=iostat_val) CA_surface
    read(*,*,iostat=iostat_val) MWA
    read(*,*,iostat=iostat_val) lengthPipe
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all convective mass transfer parameters.'
        stop
    end if

    if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. Dab <= 0.0d0 .or. MWA <= 0.0d0) then
        write(*,*) 'ERROR: Fluid properties, diffusivity, and molecular weight must be positive.'
        stop
    end if
    if (velocity < 0.0d0) then
        write(*,*) 'ERROR: Velocity cannot be negative.'
        stop
    end if

    select case (geom_type)
    case (1)
        if (L <= 0.0d0) then
            write(*,*) 'ERROR: Plate length L must be positive.'
            stop
        end if
        if (width <= 0.0d0) width = 1.0d0
        geom_name = 'External Flat Plate'
        Lc = L
        area = L * width
    case (2)
        if (Dpipe <= 0.0d0) then
            write(*,*) 'ERROR: Pipe diameter must be positive.'
            stop
        end if
        if (lengthPipe <= 0.0d0) lengthPipe = 1.0d0
        geom_name = 'Circular Pipe Internal Flow'
        Lc = Dpipe
        area = PI * Dpipe * lengthPipe
    case default
        write(*,*) 'ERROR: Invalid geometry. Use 1 flat plate or 2 pipe.'
        stop
    end select

    nu = mu / rho
    Re = rho * velocity * Lc / mu
    Sc = mu / (rho * Dab)
    dC = CA_surface - CA_bulk

    call sherwood_average(geom_type, flow_type, Re, Sc, Sh, regime_name, corr_name)
    km = Sh * Dab / Lc
    NA = km * dC
    mass_flux = NA * MWA / 1000.0d0
    total_mdot = mass_flux * area
    if (km > 0.0d0) then
        delta_m = Dab / km
    else
        delta_m = 0.0d0
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   CONVECTIVE MASS TRANSFER ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Geometry Name           = ', trim(geom_name)
    write(*,'(A,I2)')       '  Geometry Type Code      = ', geom_type
    write(*,'(A,A)')        '  Flow Regime             = ', trim(regime_name)
    write(*,'(A,A)')        '  Correlation             = ', trim(corr_name)
    write(*,'(A,ES12.4,A)') '  Characteristic Length   = ', Lc, ' m'
    write(*,'(A,ES12.4,A)') '  Transfer Area           = ', area, ' m2'
    write(*,*)
    write(*,'(A)') '--- FLUID / DIFFUSION PROPERTIES ---------------------------'
    write(*,'(A,ES12.4,A)') '  Density (rho)           = ', rho, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Dynamic Viscosity (mu)  = ', mu, ' Pa.s'
    write(*,'(A,ES12.4,A)') '  Kinematic Viscosity     = ', nu, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Binary Diffusivity      = ', Dab, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Velocity                = ', velocity, ' m/s'
    write(*,'(A,ES12.4,A)') '  Species MW_A            = ', MWA, ' g/mol'
    write(*,*)
    write(*,'(A)') '--- DIMENSIONLESS GROUPS ------------------------------------'
    write(*,'(A,ES12.4)')   '  Reynolds Number (Re)    = ', Re
    write(*,'(A,ES12.4)')   '  Schmidt Number (Sc)     = ', Sc
    write(*,'(A,ES12.4)')   '  Sherwood Number (Sh)    = ', Sh
    write(*,*)
    write(*,'(A)') '--- MASS TRANSFER RESULTS -----------------------------------'
    write(*,'(A,ES12.4,A)') '  Mass Transfer Coeff (k_m) = ', km, ' m/s'
    write(*,'(A,ES12.4,A)') '  Molar Flux (N_A)          = ', NA, ' mol/m2.s'
    write(*,'(A,ES12.4,A)') '  Mass Flux                 = ', mass_flux, ' kg/m2.s'
    write(*,'(A,ES12.4,A)') '  Total Mass Rate           = ', total_mdot, ' kg/s'
    write(*,'(A,ES12.4,A)') '  Film Thickness            = ', delta_m, ' m'
    write(*,'(A,ES12.4,A)') '  Concentration Difference  = ', dC, ' mol/m3'
    write(*,*)

    write(*,'(A)') '--- LOCAL MASS TRANSFER PROFILE -----------------------------'
    write(*,'(A)') '  x [m]       Re_x          Sh_x          k_x [m/s]     N_A [mol/m2.s] film [m]'
    write(*,'(A)') '  ---------------------------------------------------------------------------'
    n_points = 50
    do i = 1, n_points
        if (geom_type == 1) then
            x = L * dble(i) / dble(n_points)
            localLc = max(x, 1.0d-12)
            Rex = rho * velocity * localLc / mu
            call sherwood_local_plate(flow_type, Rex, Sc, Shx)
            kmx = Shx * Dab / localLc
        else
            x = lengthPipe * dble(i) / dble(n_points)
            localLc = Dpipe
            Rex = Re
            call sherwood_average(geom_type, flow_type, Re, Sc, Shx, regime_name, corr_name)
            ! Simple entrance enhancement fading with x/D for visualization.
            if (Dpipe > 0.0d0) Shx = Shx * (1.0d0 + 0.25d0 * exp(-x/(10.0d0*Dpipe)))
            kmx = Shx * Dab / Dpipe
        end if
        NAx = kmx * dC
        if (kmx > 0.0d0) then
            filmx = Dab / kmx
        else
            filmx = 0.0d0
        end if
        write(*,'(F10.5,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') x, Rex, Shx, kmx, NAx, filmx
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Re = rho u Lc / mu; Sc = mu/(rho D_AB); k_m = Sh D_AB/Lc.'
    write(*,'(A)') '  Flat plate laminar average: Sh_L = 0.664 Re_L^0.5 Sc^(1/3).'
    write(*,'(A)') '  Flat plate turbulent average: Sh_L = (0.037 Re_L^0.8 - 871) Sc^(1/3).'
    write(*,'(A)') '  Pipe turbulent Dittus-Boelter mass analogy: Sh = 0.023 Re^0.8 Sc^n.'
    write(*,'(A)') '  Film thickness: delta_m = D_AB/k_m.'

contains
    subroutine sherwood_average(g, mode, Re_in, Sc_in, Sh_out, regime, corr)
        implicit none
        integer, intent(in) :: g, mode
        double precision, intent(in) :: Re_in, Sc_in
        double precision, intent(out) :: Sh_out
        character(len=80), intent(out) :: regime, corr
        double precision :: n_exp
        n_exp = 1.0d0/3.0d0
        if (g == 1) then
            if (mode == 3 .or. (mode == 1 .and. Re_in < 5.0d5)) then
                regime = 'Laminar external boundary layer'
                corr = 'Flat plate average Sh = 0.664 Re^0.5 Sc^(1/3)'
                Sh_out = 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0)
            else
                regime = 'Turbulent / transition external boundary layer'
                corr = 'Flat plate average Sh = (0.037 Re^0.8 - 871) Sc^(1/3)'
                Sh_out = (0.037d0 * max(Re_in,0.0d0)**0.8d0 - 871.0d0) * Sc_in**(1.0d0/3.0d0)
                if (Sh_out < 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0)) &
                    Sh_out = 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0)
            end if
        else
            if (mode == 3 .or. (mode == 1 .and. Re_in < 2300.0d0)) then
                regime = 'Laminar internal flow'
                corr = 'Pipe fully developed laminar Sh = 3.66'
                Sh_out = 3.66d0
            else
                regime = 'Turbulent internal flow'
                corr = 'Dittus-Boelter mass analogy Sh = 0.023 Re^0.8 Sc^0.33'
                Sh_out = 0.023d0 * max(Re_in,0.0d0)**0.8d0 * Sc_in**n_exp
            end if
        end if
        if (Sh_out < 0.0d0) Sh_out = 0.0d0
    end subroutine sherwood_average

    subroutine sherwood_local_plate(mode, Rex, Sc_in, Shx)
        implicit none
        integer, intent(in) :: mode
        double precision, intent(in) :: Rex, Sc_in
        double precision, intent(out) :: Shx
        if (mode == 3 .or. (mode == 1 .and. Rex < 5.0d5)) then
            Shx = 0.332d0 * sqrt(max(Rex,0.0d0)) * Sc_in**(1.0d0/3.0d0)
        else
            Shx = 0.0296d0 * max(Rex,0.0d0)**0.8d0 * Sc_in**(1.0d0/3.0d0)
        end if
        if (Shx < 0.0d0) Shx = 0.0d0
    end subroutine sherwood_local_plate
end program convective_mass_transfer
