program shell_tube_tema
    implicit none

    ! Inputs
    integer :: N_shell, N_p, layout, hot_side
    double precision :: T_hot_in, T_hot_out, t_cold_in, t_cold_out
    double precision :: m_s, m_t
    double precision :: rho_s, Cp_s, mu_s, k_s
    double precision :: rho_t, Cp_t, mu_t, k_t
    double precision :: do_val, tw, L, Pt, B
    double precision :: Rfs, Rft, kw

    ! Calculated properties
    double precision :: di, K1, n1, Ds, C, As, Gs, Res, Prs, Nus, hs, fs
    double precision :: Aat, vt, Ret, Prt, Nut, ht, ft
    double precision :: Rw, U_clean, U_fouled, Q_s, Q_t, Q
    double precision :: dT1, dT2, LMTD_cf, F, P, R, K, P_1
    double precision :: denom_val, val_num, val_den, num_val
    double precision :: A_req, A_prov, dP_shell, dP_tube, De
    integer :: Nt, Nb, iostat_val, i
    logical :: has_cross_error

    ! Shell and tube temperatures mapping
    double precision :: T_hi, T_ho, T_ci, T_co
    double precision :: T_s_in, T_s_out, T_t_in, T_t_out

    double precision, parameter :: pi = 3.141592653589793d0
    double precision :: U_clean_inv, U_fouled_inv

    ! Read all input variables
    read(*,*,iostat=iostat_val) N_shell
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid N_shell.'
        stop
    end if
    read(*,*,iostat=iostat_val) N_p
    read(*,*,iostat=iostat_val) layout
    read(*,*,iostat=iostat_val) T_hot_in
    read(*,*,iostat=iostat_val) T_hot_out
    read(*,*,iostat=iostat_val) t_cold_in
    read(*,*,iostat=iostat_val) t_cold_out
    read(*,*,iostat=iostat_val) m_s
    read(*,*,iostat=iostat_val) m_t
    read(*,*,iostat=iostat_val) rho_s
    read(*,*,iostat=iostat_val) Cp_s
    read(*,*,iostat=iostat_val) mu_s
    read(*,*,iostat=iostat_val) k_s
    read(*,*,iostat=iostat_val) rho_t
    read(*,*,iostat=iostat_val) Cp_t
    read(*,*,iostat=iostat_val) mu_t
    read(*,*,iostat=iostat_val) k_t
    read(*,*,iostat=iostat_val) do_val
    read(*,*,iostat=iostat_val) tw
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) Pt
    read(*,*,iostat=iostat_val) B
    read(*,*,iostat=iostat_val) Rfs
    read(*,*,iostat=iostat_val) Rft
    read(*,*,iostat=iostat_val) kw
    read(*,*,iostat=iostat_val) hot_side

    ! Basic sanity checks
    if (N_shell < 1) N_shell = 1
    if (N_p < 1) N_p = 2
    if (layout /= 1 .and. layout /= 2) layout = 1
    if (hot_side /= 1 .and. hot_side /= 2) hot_side = 1

    ! Map process temperatures to hot/cold streams
    T_hi = T_hot_in
    T_ho = T_hot_out
    T_ci = t_cold_in
    T_co = t_cold_out

    ! Map to shell and tube sides
    if (hot_side == 1) then
        T_s_in = T_hi
        T_s_out = T_ho
        T_t_in = T_ci
        T_t_out = T_co
    else
        T_s_in = T_ci
        T_s_out = T_co
        T_t_in = T_hi
        T_t_out = T_ho
    end if

    ! Compute heat duty (average of shell and tube sides)
    Q_s = m_s * Cp_s * abs(T_s_in - T_s_out)
    Q_t = m_t * Cp_t * abs(T_t_in - T_t_out)
    Q = (Q_s + Q_t) / 2.0d0

    ! Counter-flow LMTD
    dT1 = T_hi - T_co
    dT2 = T_ho - T_ci

    has_cross_error = .false.
    if (dT1 <= 0.0d0 .or. dT2 <= 0.0d0) then
        LMTD_cf = 0.0d0
        has_cross_error = .true.
    else if (abs(dT1 - dT2) < 1.0d-6) then
        LMTD_cf = dT1
    else
        LMTD_cf = (dT1 - dT2) / log(dT1 / dT2)
    end if

    ! P and R calculations
    if (T_hi - T_ci <= 0.0d0) then
        has_cross_error = .true.
    else
        P = (T_co - T_ci) / (T_hi - T_ci)
        if (abs(T_co - T_ci) < 1.0d-8) then
            R = 1.0d0
        else
            R = (T_hi - T_ho) / (T_co - T_ci)
        end if
    end if

    F = -1.0d0
    if (.not. has_cross_error) then
        if (P < 0.0d0 .or. P >= 1.0d0 .or. R < 0.0d0) then
            has_cross_error = .true.
        else
            if (abs(R - 1.0d0) < 1.0d-5) then
                R = 1.0001d0
            end if

            if (1.0d0 - P * R <= 0.0d0 .or. 1.0d0 - P <= 0.0d0) then
                has_cross_error = .true.
            else
                K = ((1.0d0 - P * R) / (1.0d0 - P))**(1.0d0 / dble(N_shell))
                if (abs(K - R) < 1.0d-8) then
                    P_1 = P / dble(N_shell)
                else
                    P_1 = (1.0d0 - K) / (R - K)
                end if

                denom_val = 2.0d0 - P_1 * (R + 1.0d0 + sqrt(R**2 + 1.0d0))
                num_val = 2.0d0 - P_1 * (R + 1.0d0 - sqrt(R**2 + 1.0d0))

                if (denom_val <= 0.0d0 .or. num_val <= 0.0d0 .or. (1.0d0 - P_1) <= 0.0d0 .or. (1.0d0 - R*P_1) <= 0.0d0) then
                    has_cross_error = .true.
                else
                    val_num = sqrt(R**2 + 1.0d0) * log((1.0d0 - P_1) / (1.0d0 - R * P_1))
                    val_den = (R - 1.0d0) * log(num_val / denom_val)
                    if (val_den /= 0.0d0) then
                        F = val_num / val_den
                        if (F <= 0.0d0 .or. F > 1.0d0) then
                            has_cross_error = .true.
                        end if
                    else
                        has_cross_error = .true.
                    end if
                end if
            end if
        end if
    end if

    ! Iterative design logic
    if (.not. has_cross_error) then
        ! Lookup K1 and n1
        call get_bundle_constants(layout, N_p, K1, n1)

        di = do_val - 2.0d0 * tw

        ! De: Equivalent diameter
        if (layout == 1) then
            De = (2.0d0 * sqrt(3.0d0) * Pt**2 - pi * do_val**2) / (pi * do_val)
        else
            De = 4.0d0 * (Pt**2 - pi * do_val**2 / 4.0d0) / (pi * do_val)
        end if

        ! Start sizing loop
        Nt = N_p
        do while (Nt <= 100000)
            ! Shell diameter
            Ds = do_val * (dble(Nt) / K1)**(1.0d0 / n1)

            ! Shell cross-flow area
            C = Pt - do_val
            As = (Ds * C * B) / Pt

            ! Shell-side velocity, Re, Pr
            Gs = m_s / As
            Res = (Gs * De) / mu_s
            Prs = (mu_s * Cp_s) / k_s

            ! Shell film coefficient
            Nus = 0.36d0 * (Res**0.55d0) * (Prs**(1.0d0/3.0d0))
            hs = (Nus * k_s) / De

            ! Tube-side area, velocity, Re, Pr
            Aat = (dble(Nt) * pi * di**2) / (4.0d0 * dble(N_p))
            vt = m_t / (rho_t * Aat)
            Ret = (rho_t * vt * di) / mu_t
            Prt = (mu_t * Cp_t) / k_t

            ! Tube film coefficient
            if (Ret >= 2100.0d0) then
                Nut = 0.027d0 * (Ret**0.8d0) * (Prt**(1.0d0/3.0d0))
            else
                Nut = 1.86d0 * (Ret * Prt * di / L)**(1.0d0/3.0d0)
                if (Nut < 3.66d0) Nut = 3.66d0
            end if
            ht = (Nut * k_t) / di

            ! Wall conduction resistance
            Rw = (do_val * log(do_val / di)) / (2.0d0 * kw)

            ! Overall heat transfer coefficients (based on Ao)
            U_clean_inv = 1.0d0/hs + Rw + do_val / (di * ht)
            U_clean = 1.0d0 / U_clean_inv

            U_fouled_inv = 1.0d0/hs + Rfs + Rw + (do_val/di)*Rft + do_val / (di * ht)
            U_fouled = 1.0d0 / U_fouled_inv

            ! Required Area based on fouled U
            A_req = Q / (U_fouled * F * LMTD_cf)

            ! Provided Area
            A_prov = dble(Nt) * pi * do_val * L

            if (A_prov >= A_req) then
                exit
            end if

            Nt = Nt + N_p
        end do

        ! Final outputs computations
        if (Ret >= 2100.0d0) then
            ft = 0.184d0 * (Ret**(-0.2d0))
        else
            ft = 64.0d0 / Ret
        end if
        dP_tube = (ft * (L * dble(N_p) / di) + 4.0d0 * dble(N_p)) * (rho_t * vt**2) / 2.0d0

        fs = 1.44d0 * (Res**(-0.3d0))
        Nb = idint(L / B) - 1
        if (Nb < 0) Nb = 0
        dP_shell = (fs * Gs**2 * Ds * dble(Nb + 1)) / (2.0d0 * rho_s * De)
    else
        Nt = 0
        Ds = 0.0d0
        De = 0.0d0
        hs = 0.0d0
        ht = 0.0d0
        U_clean = 0.0d0
        U_fouled = 0.0d0
        dP_shell = 0.0d0
        dP_tube = 0.0d0
        A_req = 0.0d0
        A_prov = 0.0d0
        Nb = 0
        Res = 0.0d0
        Ret = 0.0d0
        Gs = 0.0d0
        vt = 0.0d0
    end if

    ! Print results to stdout in KEY=value format for easy parsing
    write(*,'(A,I2)') 'N_SHELL=', N_shell
    write(*,'(A,I2)') 'N_PASSES=', N_p
    write(*,'(A,I2)') 'LAYOUT=', layout
    write(*,'(A,F14.4)') 'T_HI=', T_hi
    write(*,'(A,F14.4)') 'T_HO=', T_ho
    write(*,'(A,F14.4)') 'T_CI=', T_ci
    write(*,'(A,F14.4)') 'T_CO=', T_co
    write(*,'(A,F14.4)') 'M_S=', m_s
    write(*,'(A,F14.4)') 'M_T=', m_t
    write(*,'(A,F14.2)') 'Q=', Q
    write(*,'(A,I6)') 'N_T=', Nt
    write(*,'(A,F14.6)') 'D_S=', Ds
    write(*,'(A,F14.4)') 'H_S=', hs
    write(*,'(A,F14.4)') 'H_T=', ht
    write(*,'(A,F14.4)') 'U_CLEAN=', U_clean
    write(*,'(A,F14.4)') 'U_FOULED=', U_fouled
    write(*,'(A,F14.4)') 'DP_SHELL=', dP_shell
    write(*,'(A,F14.4)') 'DP_TUBE=', dP_tube
    write(*,'(A,F14.4)') 'A_REQ=', A_req
    write(*,'(A,F14.4)') 'A_PROV=', A_prov
    write(*,'(A,F14.4)') 'LMTD_CF=', LMTD_cf
    write(*,'(A,F14.6)') 'F_FACTOR=', F
    write(*,'(A,F14.4)') 'RE_S=', Res
    write(*,'(A,F14.4)') 'RE_T=', Ret
    write(*,'(A,F14.4)') 'V_S=', Gs / rho_s
    write(*,'(A,F14.4)') 'V_T=', vt
    write(*,'(A,I4)') 'N_BAFFLES=', Nb
    if (has_cross_error) then
        write(*,'(A,I1)') 'HAS_CROSS_ERROR=', 1
    else
        write(*,'(A,I1)') 'HAS_CROSS_ERROR=', 0
    end if

    ! Print temperature profiles for visualization
    write(*,'(A)') '--- TEMPERATURE PROFILE ALONG EXCHANGER --------------------'
    write(*,'(A)') '  Position %    T_hot [C]    T_cold [C]'
    write(*,'(A)') '  --------------------------------------------------'
    do i = 0, 20
        write(*,'(F10.1,4X,F10.2,4X,F10.2)') dble(i)/20.0d0*100.0d0, &
            T_hi - (T_hi - T_ho) * dble(i)/20.0d0, &
            T_co - (T_co - T_ci) * dble(i)/20.0d0
    end do

contains

    subroutine get_bundle_constants(lay, Np, K1, n1)
        integer, intent(in) :: lay
        integer, intent(in) :: Np
        double precision, intent(out) :: K1, n1

        if (lay == 1) then ! Triangular 30
            select case (Np)
            case (1)
                K1 = 0.319d0; n1 = 2.142d0
            case (2)
                K1 = 0.249d0; n1 = 2.207d0
            case (4)
                K1 = 0.175d0; n1 = 2.285d0
            case (6)
                K1 = 0.0743d0; n1 = 2.499d0
            case default ! 8 or more
                K1 = 0.0365d0; n1 = 2.675d0
            end select
        else ! Square 90
            select case (Np)
            case (1)
                K1 = 0.215d0; n1 = 2.207d0
            case (2)
                K1 = 0.156d0; n1 = 2.291d0
            case (4)
                K1 = 0.158d0; n1 = 2.311d0
            case (6)
                K1 = 0.081d0; n1 = 2.473d0
            case default ! 8 or more
                K1 = 0.0321d0; n1 = 2.721d0
            end select
        end if
    end subroutine get_bundle_constants

end program shell_tube_tema
