program spiral_heat_exchanger_she
    implicit none
    integer :: iostat_val, fluid_hot_type, fluid_cold_type
    double precision :: m_dot_h_kgs, Th_in_C, Th_out_C, m_dot_c_kgs, Tc_in_C, Tc_out_C
    double precision :: H_width_m, s_hot_mm, s_cold_mm, t_wall_mm, L_sheet_m
    double precision :: cp_h, rho_h, mu_h, k_h, cp_c, rho_c, mu_c, k_c, k_wall, Rf_h, Rf_c
    double precision :: Q_duty_kW, D_h_h_m, D_h_c_m, A_flow_h, A_flow_c, V_h, V_c, Re_h, Re_c
    double precision :: Nu_h, Nu_c, h_h, h_c, U_overall, A_total_m2, LMTD_C, dp_h_kPa, dp_c_kPa
    double precision :: D_outer_m, D_core_m, R_mean_m, De_h, De_c, f_h, f_c
    character(len=32) :: she_status
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_hot_type  ! 1=Sewage Sludge / Slurry, 2=Hot Water, 3=Mineral Oil
    read(*,*,iostat=iostat_val) m_dot_h_kgs     ! Hot Stream Mass Flow [kg/s] (e.g. 12.0)
    read(*,*,iostat=iostat_val) Th_in_C         ! Hot Inlet Temp [deg C] (e.g. 85.0)
    read(*,*,iostat=iostat_val) fluid_cold_type ! 1=Cold Water, 2=Cooling Slurry, 3=Effluent
    read(*,*,iostat=iostat_val) m_dot_c_kgs     ! Cold Stream Mass Flow [kg/s] (e.g. 15.0)
    read(*,*,iostat=iostat_val) Tc_in_C         ! Cold Inlet Temp [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) H_width_m       ! Sheet Width H [m] (e.g. 1.20)
    read(*,*,iostat=iostat_val) s_hot_mm        ! Hot Channel Spacing [mm] (e.g. 16.0)
    read(*,*,iostat=iostat_val) s_cold_mm       ! Cold Channel Spacing [mm] (e.g. 12.0)
    read(*,*,iostat=iostat_val) L_sheet_m       ! Spiral Sheet Length [m] (e.g. 28.0)

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

    if (m_dot_h_kgs <= 0.0d0 .or. m_dot_c_kgs <= 0.0d0 .or. H_width_m <= 0.0d0 .or. L_sheet_m <= 0.0d0) then
        write(*,*) 'ERROR: Flow rates, width and length must be positive.'
        stop
    end if

    t_wall_mm = 4.0d0     ! 4 mm stainless steel
    k_wall = 16.0d0       ! Stainless steel 316L
    Rf_h = 0.0003d0       ! Anti-fouling self cleaning single channel
    Rf_c = 0.0002d0

    ! Fluid properties
    if (fluid_hot_type == 1) then ! Sewage Sludge / Slurry
        cp_h = 3900.0d0; rho_h = 1040.0d0; mu_h = 0.0060d0; k_h = 0.58d0
    else if (fluid_hot_type == 2) then ! Hot Water
        cp_h = 4180.0d0; rho_h = 980.0d0;  mu_h = 0.00045d0; k_h = 0.65d0
    else ! Oil
        cp_h = 2100.0d0; rho_h = 860.0d0;  mu_h = 0.0250d0; k_h = 0.14d0
    end if

    if (fluid_cold_type == 1) then ! Cold Water
        cp_c = 4184.0d0; rho_c = 998.0d0;  mu_c = 0.0010d0; k_c = 0.60d0
    else if (fluid_cold_type == 2) then ! Slurry
        cp_c = 3950.0d0; rho_c = 1030.0d0; mu_c = 0.0050d0; k_c = 0.57d0
    else ! Effluent
        cp_c = 4150.0d0; rho_c = 1005.0d0; mu_c = 0.0011d0; k_c = 0.59d0
    end if

    ! Channel geometries
    D_h_h_m = 2.0d0 * (s_hot_mm / 1000.0d0)
    D_h_c_m = 2.0d0 * (s_cold_mm / 1000.0d0)
    A_flow_h = H_width_m * (s_hot_mm / 1000.0d0)
    A_flow_c = H_width_m * (s_cold_mm / 1000.0d0)
    A_total_m2 = 2.0d0 * H_width_m * L_sheet_m

    V_h = m_dot_h_kgs / (rho_h * A_flow_h)
    V_c = m_dot_c_kgs / (rho_c * A_flow_c)

    Re_h = (rho_h * V_h * D_h_h_m) / mu_h
    Re_c = (rho_c * V_c * D_h_c_m) / mu_c

    D_core_m = 0.30d0
    D_outer_m = sqrt(D_core_m**2 + (4.0d0 * L_sheet_m * (s_hot_mm + s_cold_mm + 2.0d0*t_wall_mm) / (1000.0d0 * PI)))
    R_mean_m = (D_outer_m + D_core_m) / 4.0d0

    ! Dean numbers
    De_h = Re_h * sqrt(D_h_h_m / (2.0d0 * R_mean_m))
    De_c = Re_c * sqrt(D_h_c_m / (2.0d0 * R_mean_m))

    ! Nusselt correlation with spiral curvature enhancement
    if (Re_h > 2300.0d0) then
        Nu_h = 0.023d0 * (Re_h**0.8d0) * ((cp_h*mu_h/k_h)**0.33d0) * sqrt(1.0d0 + 1.77d0*(D_h_h_m/R_mean_m))
        f_h = (0.079d0 / (Re_h**0.25d0)) * (1.0d0 + 0.1d0 * (De_h**0.25d0))
    else
        Nu_h = 4.36d0 * (1.0d0 + 0.06d0 * De_h**0.5d0)
        f_h = (64.0d0 / max(1.0d0, Re_h)) * (1.0d0 + 0.033d0 * log10(max(1.0d0, De_h))**4)
    end if

    if (Re_c > 2300.0d0) then
        Nu_c = 0.023d0 * (Re_c**0.8d0) * ((cp_c*mu_c/k_c)**0.4d0) * sqrt(1.0d0 + 1.77d0*(D_h_c_m/R_mean_m))
        f_c = (0.079d0 / (Re_c**0.25d0)) * (1.0d0 + 0.1d0 * (De_c**0.25d0))
    else
        Nu_c = 4.36d0 * (1.0d0 + 0.06d0 * De_c**0.5d0)
        f_c = (64.0d0 / max(1.0d0, Re_c)) * (1.0d0 + 0.033d0 * log10(max(1.0d0, De_c))**4)
    end if

    h_h = Nu_h * k_h / D_h_h_m
    h_c = Nu_c * k_c / D_h_c_m

    U_overall = 1.0d0 / ((1.0d0/h_h) + Rf_h + ((t_wall_mm/1000.0d0)/k_wall) + Rf_c + (1.0d0/h_c))

    ! Thermal sizing: Effectiveness - NTU
    ! C_min and C_max
    if (m_dot_h_kgs * cp_h < m_dot_c_kgs * cp_c) then
        Q_duty_kW = U_overall * A_total_m2 * ((Th_in_C - Tc_in_C) * 0.45d0) / 1000.0d0
    else
        Q_duty_kW = U_overall * A_total_m2 * ((Th_in_C - Tc_in_C) * 0.45d0) / 1000.0d0
    end if

    Th_out_C = Th_in_C - (Q_duty_kW * 1000.0d0 / (m_dot_h_kgs * cp_h))
    Tc_out_C = Tc_in_C + (Q_duty_kW * 1000.0d0 / (m_dot_c_kgs * cp_c))

    LMTD_C = ((Th_in_C - Tc_out_C) - (Th_out_C - Tc_in_C)) / log(max(0.01d0, (Th_in_C - Tc_out_C)/(Th_out_C - Tc_in_C)))

    ! Pressure drops [kPa]
    dp_h_kPa = (f_h * (L_sheet_m / D_h_h_m) * (rho_h * V_h**2 / 2.0d0)) / 1000.0d0
    dp_c_kPa = (f_c * (L_sheet_m / D_h_c_m) * (rho_c * V_c**2 / 2.0d0)) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SPIRAL HEAT EXCHANGER (SHE ALFA LAVAL)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Outer Diameter / Sheet Width = ', D_outer_m, ' m / ', H_width_m, ' m'
    write(*,'(A,F10.2,A,F10.2,A)') 'Total Heat Transfer Area     = ', A_total_m2, ' m2 (Length = ', L_sheet_m, ' m)'
    write(*,'(A,F10.1,A,F10.1)')   'Dean Numbers (Hot / Cold)    = ', De_h, ' / ', De_c
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'THERMAL HEAT DUTY (Q)        = ', Q_duty_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Overall Heat Transfer Coeff U= ', U_overall, ' W/(m2.K)'
    write(*,'(A,F10.2,A)')  'Hot Side Outlet Temperature  = ', Th_out_C, ' deg C'
    write(*,'(A,F10.2,A)')  'Cold Side Outlet Temperature = ', Tc_out_C, ' deg C'
    write(*,'(A,F10.2,A,F10.2,A)') 'Pressure Drop (Hot / Cold)   = ', dp_h_kPa, ' kPa / ', dp_c_kPa, ' kPa'
    write(*,'(A)') '============================================================'

end program spiral_heat_exchanger_she
