program hx_pressure_drop
    implicit none

    ! Inputs
    double precision :: m_s, m_t
    double precision :: rho_s, rho_t
    double precision :: mu_s, mu_t
    double precision :: mu_ws, mu_wt
    double precision :: Ds, B, Pt
    integer :: layout
    double precision :: di, do_val, L
    integer :: Nt, Np

    ! Outputs & Intermediate variables
    double precision :: De, As, Gs, Res, vs, fs, dP_shell
    double precision :: Aat, vt, Ret, ft, dP_friction, dP_return, dP_tube_total
    integer :: Nb, iostat_val
    double precision, parameter :: pi = 3.141592653589793d0

    ! Read all variables from stdin in order
    read(*,*,iostat=iostat_val) m_s
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inputs'
        stop
    end if
    read(*,*,iostat=iostat_val) m_t
    read(*,*,iostat=iostat_val) rho_s
    read(*,*,iostat=iostat_val) rho_t
    read(*,*,iostat=iostat_val) mu_s
    read(*,*,iostat=iostat_val) mu_t
    read(*,*,iostat=iostat_val) mu_ws
    read(*,*,iostat=iostat_val) mu_wt
    read(*,*,iostat=iostat_val) Ds
    read(*,*,iostat=iostat_val) B
    read(*,*,iostat=iostat_val) Pt
    read(*,*,iostat=iostat_val) layout
    read(*,*,iostat=iostat_val) di
    read(*,*,iostat=iostat_val) do_val
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) Nt
    read(*,*,iostat=iostat_val) Np

    ! Default viscosity check
    if (mu_ws <= 0.0d0) mu_ws = mu_s
    if (mu_wt <= 0.0d0) mu_wt = mu_t

    ! Standard sanity checks to prevent division by zero
    if (Ds <= 0.0d0 .or. B <= 0.0d0 .or. Pt <= 0.0d0 .or. di <= 0.0d0 .or. do_val <= 0.0d0 .or. Nt <= 0 .or. Np <= 0) then
        write(*,*) 'ERROR: Geometrical values must be positive and non-zero.'
        stop
    end if

    ! ----------------------------------------------------
    ! A. Shell-side Pressure Drop (Kern Method)
    ! ----------------------------------------------------
    if (layout == 1) then
        ! Triangular layout 30
        De = 4.0d0 * (0.866d0 * Pt**2 - pi * do_val**2 / 8.0d0) / (pi * do_val / 2.0d0)
    else
        ! Square layout 90
        De = 4.0d0 * (Pt**2 - pi * do_val**2 / 4.0d0) / (pi * do_val)
    end if

    As = Ds * (Pt - do_val) * B / Pt
    if (As <= 0.0d0) As = 1.0d-6

    Gs = m_s / As
    vs = Gs / rho_s
    Res = Gs * De / mu_s

    if (Res > 0.0d0) then
        fs = 1.44d0 * Res**(-0.3d0)
    else
        fs = 0.0d0
    end if

    Nb = idint(L / B) - 1
    if (Nb < 0) Nb = 0

    dP_shell = fs * Gs**2 * Ds * dble(Nb + 1) / (2.0d0 * rho_s * De * (mu_s / mu_ws)**0.14d0)

    ! ----------------------------------------------------
    ! B. Tube-side Pressure Drop (TEMA Standard)
    ! ----------------------------------------------------
    Aat = dble(Nt) * pi * di**2 / (4.0d0 * dble(Np))
    if (Aat <= 0.0d0) Aat = 1.0d-6

    vt = m_t / (rho_t * Aat)
    Ret = rho_t * vt * di / mu_t

    if (Ret < 2100.0d0) then
        if (Ret > 0.0d0) then
            ft = 64.0d0 / Ret
        else
            ft = 0.0d0
        end if
    else
        ft = 0.184d0 * Ret**(-0.2d0)
    end if

    dP_friction = ft * (L * dble(Np) / di) * (rho_t * vt**2 / 2.0d0)
    dP_return = 4.0d0 * dble(Np) * (rho_t * vt**2 / 2.0d0)

    dP_tube_total = (dP_friction + dP_return) * (mu_t / mu_wt)**(-0.14d0)

    ! Print results to stdout in KEY=value format
    write(*,'(A,F18.4)') 'DE=', De
    write(*,'(A,F18.6)') 'A_S=', As
    write(*,'(A,F18.4)') 'G_S=', Gs
    write(*,'(A,F18.4)') 'RE_S=', Res
    write(*,'(A,F18.4)') 'V_S=', vs
    write(*,'(A,F18.6)') 'F_S=', fs
    write(*,'(A,I5)') 'N_BAFFLES=', Nb
    write(*,'(A,F18.4)') 'DP_SHELL=', dP_shell

    write(*,'(A,F18.6)') 'A_AT=', Aat
    write(*,'(A,F18.4)') 'V_T=', vt
    write(*,'(A,F18.4)') 'RE_T=', Ret
    write(*,'(A,F18.6)') 'F_T=', ft
    write(*,'(A,F18.4)') 'DP_FRICTION=', dP_friction
    write(*,'(A,F18.4)') 'DP_RETURN=', dP_return
    write(*,'(A,F18.4)') 'DP_TUBE_TOTAL=', dP_tube_total

end program hx_pressure_drop