Cross-Flow Ultrafiltration & Gel Polarization

program crossflow_ultrafiltration
    implicit none
    integer :: iostat_val
    double precision :: Cb_gL, Cgel_gL, dTMP_bar, u_cross_ms, dh_mm
    double precision :: Rm_m1, mu_cP, D_diff_m2s, Q_perm_Lh
    double precision :: dTMP_Pa, mu_Pas, dh_m, Re_val, Sc_val, km_ms, km_LMH
    double precision :: J_press_LMH, J_lim_LMH, J_op_LMH, Am_m2, Cp_wall_gL
    character(len=32) :: regime_status

    ! Read inputs
    read(*,*,iostat=iostat_val) Cb_gL         ! Feed Bulk Concentration [g/L] (e.g. 15.0)
    read(*,*,iostat=iostat_val) Cgel_gL       ! Gel / Limiting Concentration [g/L] (e.g. 220.0)
    read(*,*,iostat=iostat_val) dTMP_bar      ! Transmembrane Pressure TMP [bar] (e.g. 2.5)
    read(*,*,iostat=iostat_val) u_cross_ms    ! Tangential Cross-Flow Velocity [m/s] (e.g. 2.0)
    read(*,*,iostat=iostat_val) dh_mm         ! Channel Hydraulic Diameter [mm] (e.g. 4.0)
    read(*,*,iostat=iostat_val) Rm_m1         ! Clean Membrane Resistance [1/m] (e.g. 2.0e12)
    read(*,*,iostat=iostat_val) mu_cP         ! Fluid Dynamic Viscosity [cP] (e.g. 1.2)
    read(*,*,iostat=iostat_val) Q_perm_Lh     ! Required Permeate Flow [L/h] (e.g. 500.0)

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

    if (Cb_gL <= 0.0d0 .or. Cgel_gL <= Cb_gL .or. dTMP_bar <= 0.0d0 .or. u_cross_ms <= 0.0d0) then
        write(*,*) 'ERROR: Concentrations must be Cgel > Cb > 0 and TMP > 0.'
        stop
    end if

    dTMP_Pa = dTMP_bar * 1.0d5
    mu_Pas = mu_cP * 1.0d-3
    dh_m = dh_mm * 1.0d-3
    D_diff_m2s = 5.0d-11 ! Typical macromolecule / protein diffusivity

    ! Hydrodynamics & Mass transfer coefficient km
    Re_val = (1000.0d0 * u_cross_ms * dh_m) / mu_Pas
    Sc_val = mu_Pas / (1000.0d0 * D_diff_m2s)

    if (Re_val > 2100.0d0) then
        ! Turbulent Dittus-Boelter / Chilton-Colburn
        km_ms = 0.023d0 * (D_diff_m2s / dh_m) * (Re_val**0.8d0) * (Sc_val**0.33d0)
    else
        ! Laminar Leveque correlation (L ~ 1.0m)
        km_ms = 1.62d0 * ((u_cross_ms * (D_diff_m2s**2) / (dh_m * 1.0d0))**(1.0d0/3.0d0))
    end if
    km_LMH = km_ms * 3600.0d0 * 1000.0d0 ! LMH = L/(m2.h)

    ! Pressure-driven pure flux
    J_press_LMH = (dTMP_Pa / (mu_Pas * Rm_m1)) * 3600.0d0 * 1000.0d0

    ! Film theory limiting flux
    J_lim_LMH = km_LMH * log(Cgel_gL / Cb_gL)

    if (J_press_LMH < J_lim_LMH) then
        J_op_LMH = J_press_LMH
        regime_status = 'PRESSURE-CONTROLLED REGIME'
        Cp_wall_gL = Cb_gL * exp(J_op_LMH / km_LMH)
    else
        J_op_LMH = J_lim_LMH
        regime_status = 'GEL POLARIZATION LIMIT REGIME'
        Cp_wall_gL = Cgel_gL
    end if

    Am_m2 = Q_perm_Lh / max(0.1d0, J_op_LMH)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — CROSS-FLOW ULTRAFILTRATION SIZING'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Feed / Gel Concentration  = ', Cb_gL, ' g/L / ', Cgel_gL, ' g/L'
    write(*,'(A,F10.2,A,F10.2,A)') 'TMP / Cross-Flow Velocity = ', dTMP_bar, ' bar / ', u_cross_ms, ' m/s'
    write(*,'(A,F10.0,A,F10.2,A)') 'Reynolds / Mass Coeff km  = ', Re_val, ' / ', km_LMH, ' LMH'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,A)')        'Filtration Control Regime = ', trim(regime_status)
    write(*,'(A,F10.2,A)')  'Operating Permeate Flux   = ', J_op_LMH, ' LMH [L/(m2.h)]'
    write(*,'(A,F10.2,A)')  'Limiting Gel Flux (J_lim) = ', J_lim_LMH, ' LMH'
    write(*,'(A,F10.2,A)')  'Wall Concentration (Cw)   = ', Cp_wall_gL, ' g/L'
    write(*,'(A,F10.2,A)')  'Required Membrane Area Am = ', Am_m2, ' m2'
    write(*,'(A)') '============================================================'

end program crossflow_ultrafiltration