program cylinder_crossflow
  implicit none
  double precision :: D,V,Tinf,Ts,rho,mu,k,Pr,cp
  double precision :: Tf,Re,Nu_CB,Nu_H,h,Q_L,Cd,Fd
  double precision :: C_h,m_h,Vs,Res,Nus,hs,Qs
  integer :: ftype, i
  ! Read inputs
  read(*,*) D
  read(*,*) V
  read(*,*) Tinf
  read(*,*) Ts
  read(*,*) ftype
  read(*,*) rho
  read(*,*) mu
  read(*,*) k
  read(*,*) Pr
  read(*,*) cp
  Tf = (Tinf+Ts)/2.0d0
  ! Fluid defaults
  if(ftype==1) then; rho=1.177d0; mu=1.85d-5; k=0.0263d0; Pr=0.71d0; cp=1007d0
  elseif(ftype==2) then; rho=997d0; mu=8.9d-4; k=0.613d0; Pr=6.13d0; cp=4180d0
  elseif(ftype==3) then; rho=870d0; mu=0.05d0; k=0.14d0; Pr=500d0; cp=2000d0
  endif
  ! Reynolds
  Re = rho*V*D/mu
  ! Churchill-Bernstein
  Nu_CB = 0.3d0 + 0.62d0*Re**0.5d0*Pr**(1d0/3d0) &
       / (1d0+(0.4d0/Pr)**(2d0/3d0))**0.25d0 &
       * (1d0+(Re/282000d0)**(5d0/8d0))**(4d0/5d0)
  ! Hilpert
  if(Re<4d0) then; C_h=0.989d0; m_h=0.330d0
  elseif(Re<40d0) then; C_h=0.911d0; m_h=0.385d0
  elseif(Re<4000d0) then; C_h=0.683d0; m_h=0.466d0
  elseif(Re<40000d0) then; C_h=0.193d0; m_h=0.618d0
  else; C_h=0.027d0; m_h=0.805d0
  endif
  Nu_H = C_h * Re**m_h * Pr**(1d0/3d0)
  h = Nu_CB * k / D
  Q_L = h * 3.14159265d0 * D * (Ts - Tinf)
  ! Drag
  if(Re<1d0) then; Cd=24d0/Re
  elseif(Re<1000d0) then; Cd=1d0+10d0/Re**(2d0/3d0)
  else; Cd=0.4d0+6d0/Re**0.5d0
  endif
  Fd = 0.5d0*Cd*rho*V*V*D
  write(*,'(A)') '============================================'
  write(*,'(A)') '  FLOW OVER CYLINDER IN CROSS FLOW'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,F12.6,A)') '  Cylinder Diameter D     = ',D,' m'
  write(*,'(A,F12.4,A)') '  Free Stream Velocity V  = ',V,' m/s'
  write(*,'(A,F10.2,A)') '  Free Stream Temp T_inf  = ',Tinf,' C'
  write(*,'(A,F10.2,A)') '  Surface Temp T_s        = ',Ts,' C'
  write(*,'(A,F10.2,A)') '  Film Temperature T_f    = ',Tf,' C'
  write(*,'(A)') ''
  write(*,'(A)') '--- FLUID PROPERTIES ---'
  write(*,'(A,F12.4,A)') '  Density rho             = ',rho,' kg/m3'
  write(*,'(A,ES12.4,A)') '  Viscosity mu            = ',mu,' Pa.s'
  write(*,'(A,F12.6,A)') '  Conductivity k          = ',k,' W/mK'
  write(*,'(A,F12.4)')    '  Prandtl Pr              = ',Pr
  write(*,'(A)') ''
  write(*,'(A)') '--- RESULTS ---'
  write(*,'(A,ES14.4)')   '  Reynolds Number Re_D    = ',Re
  write(*,'(A,F12.2)')    '  Nu Churchill-Bernstein  = ',Nu_CB
  write(*,'(A,F12.2)')    '  Nu Hilpert              = ',Nu_H
  write(*,'(A,F12.4,A)')  '  Convection Coeff h      = ',h,' W/m2K'
  write(*,'(A,F12.2,A)')  '  Heat Transfer Q/L       = ',Q_L,' W/m'
  write(*,'(A,F10.4)')    '  Drag Coefficient Cd     = ',Cd
  write(*,'(A,F12.4,A)')  '  Drag Force F/L          = ',Fd,' N/m'
  if(Re>2d5) then
    write(*,'(A)') '  BL Separation           ~ 140 deg (turbulent)'
  else
    write(*,'(A)') '  BL Separation           ~ 80 deg (laminar)'
  endif
  write(*,'(A)') ''
  write(*,'(A)') '--- VELOCITY SWEEP ---'
  write(*,'(A)') '  V[m/s]     Re_D         Nu_CB      h[W/m2K]   Q/L[W/m]'
  write(*,'(A)') '  -----------------------------------------------------------'
  do i=1,30
    Vs = 0.1d0 + (V*5d0 - 0.1d0)*dble(i-1)/29d0
    Res = rho*Vs*D/mu
    Nus = 0.3d0 + 0.62d0*Res**0.5d0*Pr**(1d0/3d0) &
         / (1d0+(0.4d0/Pr)**(2d0/3d0))**0.25d0 &
         * (1d0+(Res/282000d0)**(5d0/8d0))**(4d0/5d0)
    hs = Nus*k/D
    Qs = hs*3.14159265d0*D*(Ts-Tinf)
    write(*,'(2X,F8.3,2X,ES11.3,2X,F10.2,2X,F10.3,2X,F12.2)') Vs,Res,Nus,hs,Qs
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Churchill-Bernstein (Re*Pr>0.2):'
  write(*,'(A)') '  Nu = 0.3 + 0.62*Re^0.5*Pr^(1/3)/[1+(0.4/Pr)^(2/3)]^0.25'
  write(*,'(A)') '       * [1+(Re/282000)^(5/8)]^(4/5)'
  write(*,'(A)') '  Ref: Incropera Ch.7 Eq.7.54'
end program cylinder_crossflow