program hx_finned
    implicit none

    ! Inputs
    integer :: FinType
    double precision :: m_a, T_in_a
    double precision :: rho_a, Cp_a, mu_a, k_a
    double precision :: m_t, T_in_t
    double precision :: rho_t, Cp_t, mu_t, k_t
    double precision :: do_val, di, L
    integer :: Nrows, Nt_row
    double precision :: Pt, Pl
    double precision :: tf, pf, kf, hf
    double precision :: ktube
    integer :: Npasses

    ! Calculated properties
    integer :: Nt, iostat_val
    double precision :: Nf, s, df, Ab, Af, A_total, r_o, r_eq
    double precision :: Amin, v_max, G_max, Re_a, Pr_a, Nu_a, h_air
    double precision :: m_val, h_term, eta_f, eta_o
    double precision :: Aat, v_t, Re_t, Pr_t, Nu_t, h_tube
    double precision :: Rw, U_inv, U_o
    double precision :: C_air, C_tube, C_min, C_max, Cr, NTU, epsilon, Q
    double precision :: T_out_a, T_out_t, fa, dP_air

    double precision, parameter :: pi = 3.141592653589793d0

    ! Read all variables from stdin in order
    read(*,*,iostat=iostat_val) FinType
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inputs'
        stop
    end if
    read(*,*,iostat=iostat_val) m_a
    read(*,*,iostat=iostat_val) T_in_a
    read(*,*,iostat=iostat_val) rho_a
    read(*,*,iostat=iostat_val) Cp_a
    read(*,*,iostat=iostat_val) mu_a
    read(*,*,iostat=iostat_val) k_a
    read(*,*,iostat=iostat_val) m_t
    read(*,*,iostat=iostat_val) T_in_t
    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) di
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) Nrows
    read(*,*,iostat=iostat_val) Nt_row
    read(*,*,iostat=iostat_val) Pt
    read(*,*,iostat=iostat_val) Pl
    read(*,*,iostat=iostat_val) tf
    read(*,*,iostat=iostat_val) pf
    read(*,*,iostat=iostat_val) kf
    read(*,*,iostat=iostat_val) hf
    read(*,*,iostat=iostat_val) ktube
    read(*,*,iostat=iostat_val) Npasses

    ! Basic sanity checks
    if (do_val <= 0.0d0 .or. di <= 0.0d0 .or. L <= 0.0d0 .or. pf <= 0.0d0 .or. tf <= 0.0d0) then
        write(*,*) 'ERROR: Geometrical tube/fin values must be positive and non-zero.'
        stop
    end if

    Nt = Nrows * Nt_row
    Nf = L / pf
    s = pf - tf
    r_o = do_val / 2.0d0

    ! Bare exposed area
    Ab = pi * do_val * (L - Nf * tf) * dble(Nt)

    ! Fin diameter and Area
    df = do_val + 2.0d0 * hf
    if (FinType == 1) then
        ! Plate Fins: Continuous
        Af = 2.0d0 * Nf * (Pt * Pl - pi * do_val**2 / 4.0d0) * dble(Nt)
        r_eq = sqrt(Pt * Pl / pi)
    else
        ! Spiral Fins: Helical
        Af = dble(Nt) * Nf * pi * ((df**2 - do_val**2) / 2.0d0 + df * tf)
        r_eq = df / 2.0d0
    end if

    A_total = Af + Ab

    ! Minimum air-flow area
    if (FinType == 1) then
        Amin = dble(Nt_row) * (Pt * L - do_val * (L - Nf * tf) - Pt * Nf * tf)
    else
        Amin = dble(Nt_row) * (Pt * L - do_val * (L - Nf * tf) - df * Nf * tf)
    end if
    if (Amin <= 0.0d0) Amin = 1.0d-6

    ! Air flow properties
    v_max = m_a / (rho_a * Amin)
    G_max = rho_a * v_max
    Re_a = G_max * do_val / mu_a
    Pr_a = mu_a * Cp_a / k_a

    ! Briggs & Young Nusselt correlation
    if (Re_a > 0.0d0) then
        Nu_a = 0.134d0 * (Re_a**0.681d0) * (Pr_a**(1.0d0/3.0d0)) * (s / tf)**0.2d0 * (s / do_val)**0.1134d0
    else
        Nu_a = 0.0d0
    end if
    h_air = Nu_a * k_a / do_val

    ! Fin efficiency
    m_val = sqrt(2.0d0 * h_air / (kf * tf))
    h_term = m_val * (r_eq - r_o)
    if (h_term > 1.0d-8) then
        eta_f = tanh(h_term) / h_term
    else
        eta_f = 1.0d0
    end if
    eta_o = 1.0d0 - (Af / A_total) * (1.0d0 - eta_f)

    ! Tube-side film coefficient
    Aat = dble(Nt) * pi * di**2 / (4.0d0 * dble(Npasses))
    if (Aat <= 0.0d0) Aat = 1.0d-6
    v_t = m_t / (rho_t * Aat)
    Re_t = rho_t * v_t * di / mu_t
    Pr_t = mu_t * Cp_t / k_t

    if (Re_t >= 2100.0d0) then
        Nu_t = 0.027d0 * (Re_t**0.8d0) * (Pr_t**(1.0d0/3.0d0))
    else
        if (Re_t > 0.0d0) then
            Nu_t = 1.86d0 * (Re_t * Pr_t * di / L)**(1.0d0/3.0d0)
        else
            Nu_t = 3.66d0
        end if
        if (Nu_t < 3.66d0) Nu_t = 3.66d0
    end if
    h_tube = Nu_t * k_t / di

    ! Tube wall resistance
    Rw = (do_val * log(do_val / di)) / (2.0d0 * ktube)

    ! Overall HTC (referred to outside area)
    U_inv = 1.0d0 / (eta_o * h_air) + Rw + do_val / (di * h_tube)
    U_o = 1.0d0 / U_inv

    ! Heat capacities
    C_air = m_a * Cp_a
    C_tube = m_t * Cp_t

    if (C_air < C_tube) then
        C_min = C_air
        C_max = C_tube
    else
        C_min = C_tube
        C_max = C_air
    end if
    Cr = C_min / C_max
    if (Cr <= 0.0d0) Cr = 1.0d-6

    ! NTU & Effectiveness for unmixed-mixed crossflow
    NTU = U_o * A_total / C_min
    if (C_air <= C_tube) then
        ! Air (mixed) has C_min
        epsilon = 1.0d0 - exp(-(1.0d0 / Cr) * (1.0d0 - exp(-Cr * NTU)))
    else
        ! Tubes (unmixed) has C_min
        epsilon = (1.0d0 / Cr) * (1.0d0 - exp(-Cr * (1.0d0 - exp(-NTU))))
    end if

    Q = epsilon * C_min * abs(T_in_t - T_in_a)

    ! Outlet temperatures
    if (T_in_t > T_in_a) then
        ! Cooling tube fluid
        T_out_t = T_in_t - Q / C_tube
        T_out_a = T_in_a + Q / C_air
    else
        ! Heating tube fluid
        T_out_t = T_in_t + Q / C_tube
        T_out_a = T_in_a - Q / C_air
    end if

    ! Air-side pressure drop (Robinson & Briggs)
    if (Re_a > 0.0d0) then
        fa = 18.93d0 * Re_a**(-0.316d0) * (Pt / do_val)**(-0.927d0)
    else
        fa = 0.0d0
    end if
    dP_air = fa * G_max**2 * dble(Nrows) / (2.0d0 * rho_a)

    ! Output results to stdout in KEY=value format
    write(*,'(A,I2)') 'FINTYPE=', FinType
    write(*,'(A,F18.4)') 'AF=', Af
    write(*,'(A,F18.4)') 'AB=', Ab
    write(*,'(A,F18.4)') 'A_TOTAL=', A_total
    write(*,'(A,F18.4)') 'AMIN=', Amin
    write(*,'(A,F18.4)') 'V_MAX_AIR=', v_max
    write(*,'(A,F18.4)') 'RE_AIR=', Re_a
    write(*,'(A,F18.4)') 'H_AIR=', h_air
    write(*,'(A,F18.4)') 'ETA_F=', eta_f
    write(*,'(A,F18.4)') 'ETA_O=', eta_o
    write(*,'(A,F18.4)') 'V_TUBE=', v_t
    write(*,'(A,F18.4)') 'RE_TUBE=', Re_t
    write(*,'(A,F18.4)') 'H_TUBE=', h_tube
    write(*,'(A,F18.4)') 'U_OVERALL=', U_o
    write(*,'(A,F18.4)') 'EPSILON=', epsilon * 100.0d0
    write(*,'(A,F18.4)') 'Q=', Q
    write(*,'(A,F18.4)') 'T_OUT_AIR=', T_out_a
    write(*,'(A,F18.4)') 'T_OUT_TUBE=', T_out_t
    write(*,'(A,F18.4)') 'DP_AIR=', dP_air
    write(*,'(A,F18.6)') 'FA=', fa

end program hx_finned
