program agitated_vessel_mixing_power
    implicit none
    integer :: iostat_val, impeller_type, has_baffles
    double precision :: N_rpm, D_imp_m, T_tank_m, H_liq_m, rho_kgm3, mu_Pas, Q_gas_vvm
    double precision :: N_revs, Re_imp, Np_turb, KL_lam, N_q, Np_calc, Power_W, Power_kW, Power_hp
    double precision :: Q_pump_m3h, blend_time_sec, tip_speed_ms, gassed_ratio
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) impeller_type    ! 1=Rushton 6-Blade, 2=Pitched Blade 45, 3=Marine Propeller, 4=Anchor, 5=Helical Ribbon
    read(*,*,iostat=iostat_val) has_baffles      ! 1=Standard 4 Baffles (T/10), 2=Unbaffled
    read(*,*,iostat=iostat_val) N_rpm            ! Impeller Speed [rpm] (e.g. 150.0)
    read(*,*,iostat=iostat_val) D_imp_m          ! Impeller Diameter D [m] (e.g. 0.40)
    read(*,*,iostat=iostat_val) T_tank_m         ! Tank Inner Diameter T [m] (e.g. 1.20)
    read(*,*,iostat=iostat_val) H_liq_m          ! Liquid Depth H [m] (e.g. 1.20)
    read(*,*,iostat=iostat_val) rho_kgm3         ! Fluid Density [kg/m3] (e.g. 1050.0)
    read(*,*,iostat=iostat_val) mu_Pas           ! Dynamic Viscosity [Pa.s] (e.g. 0.050)
    read(*,*,iostat=iostat_val) Q_gas_vvm        ! Aeration Rate [vvm] (e.g. 0.0 for ungassed, 0.8 for bioreactor)

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

    if (N_rpm <= 0.0d0 .or. D_imp_m <= 0.0d0 .or. T_tank_m <= 0.0d0 .or. rho_kgm3 <= 0.0d0) then
        write(*,*) 'ERROR: Speed, dimensions and density must be positive.'
        stop
    end if

    N_revs = N_rpm / 60.0d0
    Re_imp = (rho_kgm3 * N_revs * (D_imp_m**2)) / max(1.0d-6, mu_Pas)
    tip_speed_ms = PI * D_imp_m * N_revs

    ! Impeller specific constants
    if (impeller_type == 1) then ! Rushton Disc Turbine
        Np_turb = 5.50d0; KL_lam = 70.0d0; N_q = 0.75d0
    else if (impeller_type == 2) then ! Pitched Blade Turbine 45 deg
        Np_turb = 1.27d0; KL_lam = 45.0d0; N_q = 0.79d0
    else if (impeller_type == 3) then ! Marine Propeller
        Np_turb = 0.35d0; KL_lam = 40.0d0; N_q = 0.50d0
    else if (impeller_type == 4) then ! Anchor Impeller
        Np_turb = 0.85d0; KL_lam = 300.0d0; N_q = 0.25d0
    else ! Helical Ribbon
        Np_turb = 1.50d0; KL_lam = 600.0d0; N_q = 0.35d0
    end if

    ! Power number interpolation across laminar to turbulent
    if (Re_imp < 10.0d0) then
        Np_calc = KL_lam / max(0.1d0, Re_imp)
    else if (Re_imp > 10000.0d0) then
        Np_calc = Np_turb
    else
        Np_calc = (KL_lam / Re_imp) + Np_turb * (1.0d0 - exp(-Re_imp / 2000.0d0))
    end if

    ! Unbaffled vortex reduction
    if (has_baffles == 2 .and. Re_imp > 300.0d0) then
        Np_calc = Np_calc * 0.45d0
    end if

    ! Ungassed Shaft Power [W]
    Power_W = Np_calc * rho_kgm3 * (N_revs**3) * (D_imp_m**5)

    ! Aeration Gassing factor reduction (Michel & Miller model)
    if (Q_gas_vvm > 0.01d0 .and. impeller_type == 1) then
        gassed_ratio = 1.0d0 - 0.22d0 * log10(1.0d0 + 10.0d0 * Q_gas_vvm)
        if (gassed_ratio < 0.40d0) gassed_ratio = 0.40d0
        Power_W = Power_W * gassed_ratio
    end if

    Power_kW = Power_W / 1000.0d0
    Power_hp = Power_W / 745.7d0

    ! Pumping flow rate [m3/h]
    Q_pump_m3h = N_q * N_revs * (D_imp_m**3) * 3600.0d0

    ! 95% Blend time [sec]
    blend_time_sec = (5.4d0 / (N_revs * (Np_calc**(1.0d0/3.0d0)))) * ((T_tank_m / D_imp_m)**2)
    if (blend_time_sec < 1.0d0) blend_time_sec = 1.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — AGITATED VESSEL MIXING & IMPELLER POWER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.2,A)') 'Impeller Speed / Diameter   = ', N_rpm, ' rpm / ', D_imp_m, ' m'
    write(*,'(A,F10.1,A,F10.2,A)') 'Impeller Reynolds Number    = ', Re_imp, ' / Tip Speed = ', tip_speed_ms, ' m/s'
    write(*,'(A,F10.3)')           'Power Number (Np)           = ', Np_calc
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'IMPELLER SHAFT POWER        = ', Power_kW, ' kW (', Power_hp, ' HP)'
    write(*,'(A,F10.1,A)')         'Primary Pumping Flow Rate   = ', Q_pump_m3h, ' m3/h'
    write(*,'(A,F10.1,A)')         '95% Homogenization Blend Time= ', blend_time_sec, ' seconds'
    write(*,'(A)') '============================================================'

end program agitated_vessel_mixing_power