Agitator Power & Mixing Time

program agitator_power_mixing
    implicit none
    integer :: iostat_val, impeller_type, is_baffled
    double precision :: T_tank_m, H_liq_m, D_imp_m, N_rpm, rho_kgm3, mu_cP
    double precision :: N_rps, mu_Pas, Re_imp, Np_val, Nq_val, P_shaft_W, P_shaft_kW
    double precision :: V_tank_m3, P_per_vol_Wm3, Q_pump_m3h, theta95_s
    character(len=32) :: imp_name
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) T_tank_m        ! Tank Inside Diameter [m]
    read(*,*,iostat=iostat_val) H_liq_m         ! Liquid Height [m]
    read(*,*,iostat=iostat_val) D_imp_m         ! Impeller Diameter [m]
    read(*,*,iostat=iostat_val) N_rpm           ! Agitator Speed [RPM]
    read(*,*,iostat=iostat_val) rho_kgm3        ! Fluid Density [kg/m3]
    read(*,*,iostat=iostat_val) mu_cP           ! Dynamic Viscosity [cP]
    read(*,*,iostat=iostat_val) impeller_type   ! 1=Rushton, 2=Hydrofoil, 3=Pitched Blade, 4=Anchor
    read(*,*,iostat=iostat_val) is_baffled      ! 1=Standard 4 Baffles, 0=Unbaffled

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

    ! Conversions
    N_rps = N_rpm / 60.0d0
    mu_Pas = mu_cP * 1.0d-3
    V_tank_m3 = (PI / 4.0d0) * (T_tank_m**2) * H_liq_m

    ! Impeller Reynolds Number
    Re_imp = (rho_kgm3 * N_rps * (D_imp_m**2)) / max(1.0d-5, mu_Pas)

    ! Impeller Power & Pumping Numbers (Np, Nq)
    if (impeller_type == 1) then
        imp_name = 'Rushton 6-Blade Turbine (Radial)'
        if (Re_imp < 10.0d0) then
            Np_val = 70.0d0 / max(1.0d-2, Re_imp)
        else if (Re_imp < 1.0d4) then
            Np_val = 5.0d0 * (1.0d0 + 10.0d0 / sqrt(Re_imp))
        else
            Np_val = 5.0d0
        end if
        Nq_val = 0.75d0

    else if (impeller_type == 2) then
        imp_name = 'Hydrofoil / Marine Propeller (Axial)'
        if (Re_imp < 10.0d0) then
            Np_val = 40.0d0 / max(1.0d-2, Re_imp)
        else
            Np_val = 0.35d0
        end if
        Nq_val = 0.55d0

    else if (impeller_type == 3) then
        imp_name = '4-Pitched-Blade Turbine 45 deg (Mixed)'
        if (Re_imp < 10.0d0) then
            Np_val = 50.0d0 / max(1.0d-2, Re_imp)
        else
            Np_val = 1.27d0
        end if
        Nq_val = 0.79d0

    else
        imp_name = 'Anchor Impeller (Viscous Laminar)'
        Np_val = max(0.8d0, 350.0d0 / max(1.0d-1, Re_imp))
        Nq_val = 0.30d0
    end if

    ! Unbaffled tank power reduction
    if (is_baffled == 0 .and. Re_imp > 300.0d0) then
        Np_val = Np_val * 0.65d0
    end if

    ! Shaft Power Draw [Watts]
    P_shaft_W = Np_val * rho_kgm3 * (N_rps**3) * (D_imp_m**5)
    P_shaft_kW = P_shaft_W * 1.0d-3
    P_per_vol_Wm3 = P_shaft_W / max(1.0d-3, V_tank_m3)

    ! Pumping Rate [m3/h]
    Q_pump_m3h = Nq_val * N_rps * (D_imp_m**3) * 3600.0d0

    ! Grenville 95% Mixing Time [seconds]
    if (Re_imp > 5000.0d0) then
        theta95_s = (5.4d0 / (N_rps * (Np_val**(1.0d0/3.0d0)))) * ((T_tank_m / D_imp_m)**2)
    else
        theta95_s = (300.0d0 / (N_rps * sqrt(max(1.0d0, Re_imp)))) * ((T_tank_m / D_imp_m)**2)
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — INDUSTRIAL AGITATOR & MIXING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,A)')        'Impeller Configuration    = ', trim(imp_name)
    write(*,'(A,F10.2,A)')  'Tank Diameter T           = ', T_tank_m, ' m'
    write(*,'(A,F10.2,A)')  'Impeller Diameter D       = ', D_imp_m, ' m'
    write(*,'(A,F10.1,A)')  'Agitator Speed N          = ', N_rpm, ' RPM'
    write(*,'(A,ES12.4)')   'Impeller Reynolds Re      = ', Re_imp
    write(*,'(A,F10.2)')    'Power Number Np           = ', Np_val
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Shaft Mechanical Power    = ', P_shaft_kW, ' kW'
    write(*,'(A,F10.2,A)')  'Specific Power (P/V)      = ', P_per_vol_Wm3, ' W/m3'
    write(*,'(A,F10.1,A)')  'Liquid Pumping Rate Qp    = ', Q_pump_m3h, ' m3/h'
    write(*,'(A,F10.2,A)')  '95% Blending Mixing Time  = ', theta95_s, ' seconds'
    write(*,'(A,F10.2,A)')  'Recommended Motor (1.25x) = ', P_shaft_kW * 1.25d0, ' kW'
    write(*,'(A)') '============================================================'

end program agitator_power_mixing