program scraped_surface_heat_exchanger_sshe
    implicit none
    integer :: iostat_val, product_type, num_blades
    double precision :: m_dot_kgs, Tin_C, Tout_C, T_jacket_C, N_rpm, D_inner_mm, L_cylinder_m
    double precision :: rho_p, cp_p, mu_p, k_p, latent_heat_kJkg, k_wall, t_wall_mm, h_jacket
    double precision :: D_i_m, A_heat_m2, V_ax, Re_rot, Pr_p, Nu_sshe, h_prod, U_overall
    double precision :: Q_duty_kW, LMTD_C, P_shaft_kW, shear_rate_s, tau_visc_Pa
    character(len=32) :: app_regime
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) product_type    ! 1=Ice Cream / Frozen Slush, 2=Tomato Paste / Puree, 3=Margarine / Palm Fat, 4=Gelatin / Starch
    read(*,*,iostat=iostat_val) m_dot_kgs       ! Product Mass Flow Rate [kg/s] (e.g. 0.80)
    read(*,*,iostat=iostat_val) Tin_C           ! Product Inlet Temp [deg C] (e.g. 4.0)
    read(*,*,iostat=iostat_val) T_jacket_C      ! Cooling/Heating Jacket Temp [deg C] (e.g. -20.0)
    read(*,*,iostat=iostat_val) N_rpm           ! Scraper Shaft Rotational Speed [rpm] (e.g. 350.0)
    read(*,*,iostat=iostat_val) num_blades      ! Number of Scraper Blades (e.g. 2, 3, or 4)
    read(*,*,iostat=iostat_val) D_inner_mm      ! Cylinder Inner Diameter [mm] (e.g. 150.0)
    read(*,*,iostat=iostat_val) L_cylinder_m    ! Cylinder Length [m] (e.g. 1.80)

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

    if (m_dot_kgs <= 0.0d0 .or. N_rpm <= 0.0d0 .or. D_inner_mm <= 0.0d0 .or. L_cylinder_m <= 0.0d0) then
        write(*,*) 'ERROR: Flow rate, rotational speed, and dimensions must be positive.'
        stop
    end if

    ! Product thermophysical properties
    if (product_type == 1) then ! Ice Cream / Frozen Slush
        rho_p = 1080.0d0; cp_p = 3300.0d0; mu_p = 1.20d0; k_p = 0.45d0
        latent_heat_kJkg = 45.0d0; app_regime = 'FREEZING CRYSTALLIZATION'
    else if (product_type == 2) then ! Tomato Paste
        rho_p = 1150.0d0; cp_p = 3600.0d0; mu_p = 2.50d0; k_p = 0.52d0
        latent_heat_kJkg = 0.0d0;  app_regime = 'HIGH VISCOSITY HEATING'
    else if (product_type == 3) then ! Margarine / Palm Fat
        rho_p = 920.0d0;  cp_p = 2200.0d0; mu_p = 0.85d0; k_p = 0.19d0
        latent_heat_kJkg = 35.0d0; app_regime = 'FAT CRYSTALLIZATION'
    else ! Gelatin
        rho_p = 1050.0d0; cp_p = 3800.0d0; mu_p = 0.65d0; k_p = 0.55d0
        latent_heat_kJkg = 0.0d0;  app_regime = 'VISCOUS COOLING'
    end if

    D_i_m = D_inner_mm / 1000.0d0
    t_wall_mm = 5.0d0
    k_wall = 16.0d0     ! 316L SS
    h_jacket = 3500.0d0 ! Evaporating ammonia or steam jacket
    A_heat_m2 = PI * D_i_m * L_cylinder_m

    ! Rotational & Scraping Fluid Dynamics
    Re_rot = (rho_p * (N_rpm / 60.0d0) * (D_i_m**2)) / mu_p
    Pr_p = (cp_p * mu_p) / k_p

    ! Trommelen & Beek Scraped Surface Nusselt Correlation
    Nu_sshe = 0.16d0 * (Re_rot**0.50d0) * (Pr_p**0.33d0) * (dble(num_blades)**0.33d0)
    h_prod = Nu_sshe * k_p / D_i_m

    ! Overall Heat Transfer Coefficient U [W/(m2.K)]
    U_overall = 1.0d0 / ((1.0d0 / h_prod) + ((t_wall_mm/1000.0d0) / k_wall) + (1.0d0 / h_jacket))

    ! Thermal sizing (NTU method with constant jacket temperature)
    ! NTU = U * A / (m_dot * cp)
    ! Tout = T_jacket + (Tin - T_jacket) * exp(-NTU)
    Tout_C = T_jacket_C + (Tin_C - T_jacket_C) * exp(-(U_overall * A_heat_m2) / (m_dot_kgs * cp_p))

    Q_duty_kW = (m_dot_kgs * cp_p * abs(Tin_C - Tout_C) / 1000.0d0) + (m_dot_kgs * latent_heat_kJkg)

    LMTD_C = abs(Tin_C - Tout_C) / log(max(0.01d0, abs((Tin_C - T_jacket_C)/(Tout_C - T_jacket_C))))

    ! Shaft mechanical power demand
    shear_rate_s = PI * D_i_m * (N_rpm / 60.0d0) / 0.005d0
    tau_visc_Pa = mu_p * shear_rate_s
    P_shaft_kW = (tau_visc_Pa * A_heat_m2 * (PI * D_i_m * N_rpm / 60.0d0)) / 1000.0d0 + 0.35d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SCRAPED SURFACE HEAT EXCHANGER (SSHE)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.1,A)') 'Scraper Speed / Blades Count = ', N_rpm, ' rpm / ', dble(num_blades), ' blades'
    write(*,'(A,F10.2,A,F10.1,A)') 'Heat Transfer Area           = ', A_heat_m2, ' m2 (Length = ', L_cylinder_m, ' m)'
    write(*,'(A,F10.1,A,A)')       'Rotational Reynolds (Re_r)   = ', Re_rot, ' | ', trim(app_regime)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'THERMAL HEAT DUTY (Q)        = ', Q_duty_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Product Scraped Film Coeff hp= ', h_prod, ' W/(m2.K)'
    write(*,'(A,F10.1,A)')  'Overall Heat Transfer Coeff U= ', U_overall, ' W/(m2.K)'
    write(*,'(A,F10.2,A)')  'Product Outlet Temperature   = ', Tout_C, ' deg C'
    write(*,'(A,F10.2,A)')  'Log Mean Temp Difference LMTD= ', LMTD_C, ' deg C'
    write(*,'(A,F10.2,A)')  'Shaft Motor Power Demand     = ', P_shaft_kW, ' kW'
    write(*,'(A)') '============================================================'

end program scraped_surface_heat_exchanger_sshe