program hx_plate
    implicit none

    ! Inputs
    integer :: mode, Nplates
    double precision :: T_hi, T_ho, T_ci, T_co
    double precision :: m_h, m_c
    double precision :: rho_h, Cp_h, mu_h, k_h, mu_wh
    double precision :: rho_c, Cp_c, mu_c, k_c, mu_wc
    double precision :: W, Lp, t_plate, b_depth, pc_pitch, beta_deg, d_port_mm
    double precision :: dP_max_kpa
    double precision :: Rfh, Rfc, kw

    ! Calculated properties
    double precision :: X, phi, Dh, Q_h, Q_c, Q, dT1, dT2, LMTD_cf
    double precision :: hs_h, hs_c, dP_h, dP_c, dP_port_h, dP_port_c, U_f, U_c, A_prov, A_req
    double precision :: Re_h, Re_c, vch_h, vch_c
    integer :: Nplates_curr, iostat_val
    logical :: has_cross_error

    double precision, parameter :: pi = 3.141592653589793d0
    double precision :: dP_max

    ! Read all inputs from stdin
    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid inputs'
        stop
    end if
    read(*,*,iostat=iostat_val) Nplates
    read(*,*,iostat=iostat_val) T_hi
    read(*,*,iostat=iostat_val) T_ho
    read(*,*,iostat=iostat_val) T_ci
    read(*,*,iostat=iostat_val) T_co
    read(*,*,iostat=iostat_val) m_h
    read(*,*,iostat=iostat_val) m_c
    read(*,*,iostat=iostat_val) rho_h
    read(*,*,iostat=iostat_val) Cp_h
    read(*,*,iostat=iostat_val) mu_h
    read(*,*,iostat=iostat_val) k_h
    read(*,*,iostat=iostat_val) mu_wh
    read(*,*,iostat=iostat_val) rho_c
    read(*,*,iostat=iostat_val) Cp_c
    read(*,*,iostat=iostat_val) mu_c
    read(*,*,iostat=iostat_val) k_c
    read(*,*,iostat=iostat_val) mu_wc
    read(*,*,iostat=iostat_val) W
    read(*,*,iostat=iostat_val) Lp
    read(*,*,iostat=iostat_val) t_plate
    read(*,*,iostat=iostat_val) b_depth
    read(*,*,iostat=iostat_val) pc_pitch
    read(*,*,iostat=iostat_val) beta_deg
    read(*,*,iostat=iostat_val) d_port_mm
    read(*,*,iostat=iostat_val) dP_max_kpa
    read(*,*,iostat=iostat_val) Rfh
    read(*,*,iostat=iostat_val) Rfc
    read(*,*,iostat=iostat_val) kw

    ! Defaults handling
    if (mu_wh <= 0.0d0) mu_wh = mu_h
    if (mu_wc <= 0.0d0) mu_wc = mu_c

    ! Basic sanity checks
    if (W <= 0.0d0 .or. Lp <= 0.0d0 .or. b_depth <= 0.0d0 .or. pc_pitch <= 0.0d0 .or. beta_deg <= 0.0d0) then
        write(*,*) 'ERROR: Geometrical plate values must be positive and non-zero.'
        stop
    end if

    ! Compute heat duty (average of hot and cold sides)
    Q_h = m_h * Cp_h * abs(T_hi - T_ho)
    Q_c = m_c * Cp_c * abs(T_co - T_ci)
    Q = (Q_h + Q_c) / 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

    if (has_cross_error) then
        write(*,*) 'ERROR: Temperature cross error detected. LMTD is undefined.'
        stop
    end if

    ! Convert dP_max from kPa to Pa
    dP_max = dP_max_kpa * 1000.0d0

    ! Initial geometrical constants
    X = (pi * b_depth) / pc_pitch
    phi = (1.0d0 / 6.0d0) * (1.0d0 + sqrt(1.0d0 + X**2) + 4.0d0 * sqrt(1.0d0 + X**2 / 2.0d0))
    Dh = 2.0d0 * b_depth / phi

    if (mode == 1) then
        ! Verification mode: evaluate at input plate count
        call calculate_phe(Nplates, hs_h, hs_c, dP_h, dP_c, &
            dP_port_h, dP_port_c, U_f, U_c, A_prov, A_req, &
            Re_h, Re_c, vch_h, vch_c)
    else
        ! Sizing mode: iterate plate count to meet area and pressure limits
        Nplates_curr = Nplates
        if (mod(Nplates_curr, 2) == 0) Nplates_curr = Nplates_curr + 1
        if (Nplates_curr < 5) Nplates_curr = 5

        do while (Nplates_curr <= 50000)
            call calculate_phe(Nplates_curr, hs_h, hs_c, dP_h, dP_c, &
                dP_port_h, dP_port_c, U_f, U_c, A_prov, A_req, &
                Re_h, Re_c, vch_h, vch_c)
            if (A_prov >= A_req .and. dP_h <= dP_max .and. dP_c <= dP_max) then
                Nplates = Nplates_curr
                exit
            end if
            Nplates_curr = Nplates_curr + 2
        end do

        if (Nplates_curr > 50000) then
            write(*,*) 'ERROR: Plate sizing loop did not converge within 50,000 plates.'
            stop
        end if
    end if

    ! Print outputs to stdout
    write(*,'(A,I5)') 'N_PLATES=', Nplates
    write(*,'(A,F18.4)') 'ENLARGEMENT_PHI=', phi
    write(*,'(A,F18.6)') 'DH=', Dh
    write(*,'(A,F18.4)') 'Q=', Q
    write(*,'(A,F18.4)') 'LMTD=', LMTD_cf
    write(*,'(A,F18.4)') 'A_PROV=', A_prov
    write(*,'(A,F18.4)') 'A_REQ=', A_req
    write(*,'(A,F18.4)') 'U_CLEAN=', U_c
    write(*,'(A,F18.4)') 'U_FOULED=', U_f
    write(*,'(A,F18.4)') 'H_HOT=', hs_h
    write(*,'(A,F18.4)') 'H_COLD=', hs_c
    write(*,'(A,F18.4)') 'RE_HOT=', Re_h
    write(*,'(A,F18.4)') 'RE_COLD=', Re_c
    write(*,'(A,F18.4)') 'V_HOT=', vch_h
    write(*,'(A,F18.4)') 'V_COLD=', vch_c
    write(*,'(A,F18.4)') 'DP_HOT_TOTAL=', dP_h
    write(*,'(A,F18.4)') 'DP_COLD_TOTAL=', dP_c
    write(*,'(A,F18.4)') 'DP_HOT_PORT=', dP_port_h
    write(*,'(A,F18.4)') 'DP_COLD_PORT=', dP_port_c
    write(*,'(A,F18.4)') 'DP_HOT_CHAN=', dP_h - dP_port_h
    write(*,'(A,F18.4)') 'DP_COLD_CHAN=', dP_c - dP_port_c

contains

    subroutine calculate_phe(Np_val, hs_h, hs_c, dP_h, dP_c, &
            dP_port_h, dP_port_c, U_f, U_c, A_prov_val, A_req_val, &
            Re_h_val, Re_c_val, vch_h_val, vch_c_val)
        integer, intent(in) :: Np_val
        double precision, intent(out) :: hs_h, hs_c, dP_h, dP_c, dP_port_h, dP_port_c, U_f, U_c, A_prov_val, A_req_val
        double precision, intent(out) :: Re_h_val, Re_c_val, vch_h_val, vch_c_val

        integer :: Nch
        double precision :: Ach, Pr_h, Pr_c, Kf_h, p_h, f_turb_h, f_lam_h, f_h
        double precision :: Kf_c, p_c, f_turb_c, f_lam_c, f_c
        double precision :: beta_rad, Nu_h, Nu_c, Rw
        double precision :: U_clean_inv, U_fouled_inv, d_port, v_port_h, v_port_c
        double precision :: dP_channel_h, dP_channel_c

        Nch = (Np_val - 1) / 2
        if (Nch < 1) Nch = 1

        Ach = b_depth * W
        vch_h_val = m_h / (rho_h * Ach * dble(Nch))
        vch_c_val = m_c / (rho_c * Ach * dble(Nch))

        Re_h_val = rho_h * vch_h_val * Dh / mu_h
        Re_c_val = rho_c * vch_c_val * Dh / mu_c

        Pr_h = mu_h * Cp_h / k_h
        Pr_c = mu_c * Cp_c / k_c

        ! Chevron friction factors smooth continuous fit
        Kf_h = 0.15d0 + 0.08d0 * exp(0.06d0 * beta_deg)
        p_h = 0.30d0 - 0.002d0 * beta_deg
        f_turb_h = Kf_h * Re_h_val**(-p_h)
        f_lam_h = (10.0d0 + 0.7d0 * beta_deg) / Re_h_val
        f_h = f_turb_h
        if (f_lam_h > f_h) f_h = f_lam_h

        Kf_c = 0.15d0 + 0.08d0 * exp(0.06d0 * beta_deg)
        p_c = 0.30d0 - 0.002d0 * beta_deg
        f_turb_c = Kf_c * Re_c_val**(-p_c)
        f_lam_c = (10.0d0 + 0.7d0 * beta_deg) / Re_c_val
        f_c = f_turb_c
        if (f_lam_c > f_c) f_c = f_lam_c

        ! Nusselt numbers
        beta_rad = beta_deg * pi / 180.0d0
        Nu_h = 0.205d0 * Pr_h**(1.0d0/3.0d0) * (mu_h / mu_wh)**0.14d0 * &
            (f_h * Re_h_val**2 * sin(2.0d0 * beta_rad))**(1.0d0/3.0d0)
        Nu_c = 0.205d0 * Pr_c**(1.0d0/3.0d0) * (mu_c / mu_wc)**0.14d0 * &
            (f_c * Re_c_val**2 * sin(2.0d0 * beta_rad))**(1.0d0/3.0d0)

        hs_h = Nu_h * k_h / Dh
        hs_c = Nu_c * k_c / Dh

        Rw = t_plate / kw

        U_clean_inv = 1.0d0 / hs_h + Rw + 1.0d0 / hs_c
        U_c = 1.0d0 / U_clean_inv

        U_fouled_inv = 1.0d0 / hs_h + Rfh + Rw + Rfc + 1.0d0 / hs_c
        U_f = 1.0d0 / U_fouled_inv

        A_req_val = Q / (U_f * LMTD_cf)
        A_prov_val = dble(Np_val - 2) * W * Lp * phi

        ! Port velocity
        d_port = d_port_mm / 1000.0d0
        v_port_h = m_h / (rho_h * pi * d_port**2 / 4.0d0)
        v_port_c = m_c / (rho_c * pi * d_port**2 / 4.0d0)

        ! Pressure drops
        dP_channel_h = f_h * (Lp / Dh) * (rho_h * vch_h_val**2 / 2.0d0) * (mu_h / mu_wh)**(-0.14d0)
        dP_channel_c = f_c * (Lp / Dh) * (rho_c * vch_c_val**2 / 2.0d0) * (mu_c / mu_wc)**(-0.14d0)

        dP_port_h = 1.4d0 * 1.0d0 * (rho_h * v_port_h**2 / 2.0d0)
        dP_port_c = 1.4d0 * 1.0d0 * (rho_c * v_port_c**2 / 2.0d0)

        dP_h = dP_channel_h + dP_port_h
        dP_c = dP_channel_c + dP_port_c

    end subroutine calculate_phe

end program hx_plate
