program combined_convection
    implicit none
    double precision :: D, V, Ti, Ts, Lc, rho, mu, kf, Pr, cp, beta
    integer :: or_type, fl_type
    
    double precision :: Tf, dT, nu, alpha, Re, Gr, Ra, Ri, Nuf, Nun, Nuc, h, g
    double precision :: temp_val
    integer :: iostat_val, i
    double precision :: Vs, Res, Grs, Ris, Nufs, Nuns, Nucs, hs
    
    g = 9.80665d0
    
    ! Read inputs sequentially
    read(*,*,iostat=iostat_val) D
    read(*,*,iostat=iostat_val) V
    read(*,*,iostat=iostat_val) Ti
    read(*,*,iostat=iostat_val) Ts
    read(*,*,iostat=iostat_val) Lc
    read(*,*,iostat=iostat_val) or_type
    read(*,*,iostat=iostat_val) fl_type
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) kf
    read(*,*,iostat=iostat_val) Pr
    read(*,*,iostat=iostat_val) cp
    read(*,*,iostat=iostat_val) temp_val ! read the 13th value which is beta from php, or custom beta
    
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all input values.'
        stop
    end if
    
    Tf = (Ti + Ts) / 2.0d0
    dT = abs(Ts - Ti)
    
    ! Select defaults based on fluid type
    if (fl_type == 1) then
        ! Air
        rho = 1.177d0
        mu = 1.85d-5
        kf = 0.0263d0
        Pr = 0.71d0
        cp = 1007d0
        beta = 1.0d0 / (Tf + 273.15d0)
    else if (fl_type == 2) then
        ! Water
        rho = 997d0
        mu = 8.9d-4
        kf = 0.613d0
        Pr = 6.13d0
        cp = 4180d0
        beta = 2.1d-4
    else
        ! Engine Oil
        rho = 870d0
        mu = 0.05d0
        kf = 0.14d0
        Pr = 500d0
        cp = 2000d0
        beta = 7.0d-4
    end if
    
    ! Override if user provided custom beta
    if (temp_val > 0.0d0) beta = temp_val
    
    nu = mu / rho
    alpha = kf / (rho * cp)
    
    Re = rho * V * D / mu
    Gr = g * beta * dT * Lc**3 / nu**2
    Ra = Gr * Pr
    
    if (Re > 1.0d-10) then
        Ri = Gr / Re**2
    else
        Ri = 9999.0d0
    end if
    
    ! Forced Convection Nusselt Nuf (Flat Plate model)
    if (Re < 5.0d5) then
        Nuf = 0.664d0 * Re**0.5d0 * Pr**(1.0d0/3.0d0)
    else
        Nuf = (0.037d0 * Re**0.8d0 - 871.0d0) * Pr**(1.0d0/3.0d0)
    end if
    if (Nuf < 0.0d0) Nuf = 0.0d0
    
    ! Natural Convection Nusselt Nun (Vertical Plate model)
    Nun = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
           ((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
           
    ! Combined Nusselt Nuc
    if (or_type == 1) then
        ! Assisting
        Nuc = (Nuf**3 + Nun**3)**(1.0d0/3.0d0)
    else if (or_type == 2) then
        ! Opposing
        Nuc = abs(Nuf**3 - Nun**3)**(1.0d0/3.0d0)
    else
        ! Transverse / Cross Flow
        Nuc = (Nuf**3 + Nun**3)**(1.0d0/3.0d0)
    end if
    
    h = Nuc * kf / Lc
    
    write(*,'(A)') '============================================'
    write(*,'(A)') '  COMBINED FREE & FORCED CONVECTION ENGINE'
    write(*,'(A)') '============================================'
    write(*,*)
    write(*,'(A,ES14.4)') '  Reynolds Re             = ', Re
    write(*,'(A,ES14.4)') '  Grashof Number          = ', Gr
    write(*,'(A,ES14.4)') '  Richardson Number Ri    = ', Ri
    write(*,'(A,F12.4)')    '  Nu Forced               = ', Nuf
    write(*,'(A,F12.4)')    '  Nu Natural              = ', Nun
    write(*,'(A,F12.4)')    '  Nu Combined             = ', Nuc
    write(*,'(A,F12.4,A)')  '  Coeff h                 = ', h, ' W/m2K'
    write(*,*)
    
    write(*,'(A)') '--- VELOCITY SWEEP ---'
    write(*,'(A)') '  V[m/s]     Re           Gr           Ri           Nu_f       Nu_n       Nu_c       h[W/m2K]'
    write(*,'(A)') '  ------------------------------------------------------------------------------------------------'
    do i=1,25
        Vs = 0.01d0 + (max(V,0.2d0)*3.0d0 - 0.01d0)*dble(i-1)/24.0d0
        Res = rho * Vs * D / mu
        Grs = g * beta * dT * Lc**3 / nu**2
        if (Res > 1.0d-10) then
            Ris = Grs / Res**2
        else
            Ris = 9999.0d0
        end if
        
        if (Res < 5.0d5) then
            Nufs = 0.664d0 * Res**0.5d0 * Pr**(1.0d0/3.0d0)
        else
            Nufs = (0.037d0 * Res**0.8d0 - 871.0d0) * Pr**(1.0d0/3.0d0)
        end if
        if (Nufs < 0.0d0) Nufs = 0.0d0
        
        Nuns = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
               ((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
               
        if (or_type == 1) then
            Nucs = (Nufs**3 + Nuns**3)**(1.0d0/3.0d0)
        else if (or_type == 2) then
            Nucs = abs(Nufs**3 - Nuns**3)**(1.0d0/3.0d0)
        else
            Nucs = (Nufs**3 + Nuns**3)**(1.0d0/3.0d0)
        end if
        
        hs = Nucs * kf / Lc
        write(*,'(2X,F8.3,2X,ES10.3,2X,ES10.3,2X,ES10.3,2X,F9.2,2X,F9.2,2X,F9.2,2X,F10.3)') Vs, Res, Grs, Ris, Nufs, Nuns, Nucs, hs
    end do
    write(*,*)
end program combined_convection
