program tube_banks
  implicit none
  double precision :: D,ST,SL,V,Tinf,Ts,rho,mu,k,Pr,cp
  double precision :: Vmax,AT,AD,SD2,Re_max,Nu,h,dP_row,dP
  double precision :: C1,m1,C2,Vs,Res,Nus,hs,dPs
  double precision :: STD,SLD,chi,f_z
  double precision :: C2_tab(8),N_tab(8)
  integer :: config,Nrows,ftype,i,j
  data N_tab/1d0,2d0,3d0,5d0,7d0,10d0,13d0,20d0/
  data C2_tab/0.64d0,0.76d0,0.84d0,0.92d0,0.95d0,0.97d0,0.99d0,1.0d0/
  read(*,*) config
  read(*,*) Nrows
  read(*,*) D
  read(*,*) ST
  read(*,*) SL
  read(*,*) V
  read(*,*) Tinf
  read(*,*) Ts
  read(*,*) ftype
  read(*,*) rho,mu,k,Pr,cp
  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
  STD = ST/D; SLD = SL/D
  AT = ST - D
  if(config==1) then
    Vmax = V * ST / AT
  else
    SD2 = sqrt(SL*SL + (ST/2d0)*(ST/2d0))
    AD = SD2 - D
    if(AD < AT) then
      Vmax = V * ST / (2d0*AD)
    else
      Vmax = V * ST / AT
    endif
  endif
  Re_max = rho*Vmax*D/mu
  ! Zukauskas C1,m
  if(config==1) then
    if(Re_max<100d0) then; C1=0.9d0; m1=0.4d0
    elseif(Re_max<1000d0) then; C1=0.52d0; m1=0.5d0
    elseif(Re_max<2d5) then; C1=0.27d0; m1=0.63d0
    else; C1=0.033d0; m1=0.8d0; endif
  else
    if(Re_max<500d0) then; C1=1.04d0; m1=0.4d0
    elseif(Re_max<1000d0) then; C1=0.71d0; m1=0.5d0
    elseif(Re_max<2d5) then; C1=0.35d0; m1=0.6d0
    else; C1=0.031d0; m1=0.8d0; endif
  endif
  ! Row correction C2
  C2 = 1.0d0
  if(Nrows < 20) then
    do j=1,7
      if(dble(Nrows)>=N_tab(j) .and. dble(Nrows)<N_tab(j+1)) then
        C2 = C2_tab(j) + (C2_tab(j+1)-C2_tab(j)) &
           *(dble(Nrows)-N_tab(j))/(N_tab(j+1)-N_tab(j))
        exit
      endif
    enddo
    if(Nrows<=1) C2=C2_tab(1)
  endif
  Nu = C1 * C2 * Re_max**m1 * Pr**0.36d0
  h = Nu * k / D
  ! Pressure drop (simplified Zukauskas)
  chi = 1.0d0
  if(config==1) then
    f_z = (0.25d0 + 0.1175d0/(STD-1d0)**1.08d0) * Re_max**(-0.16d0)
  else
    f_z = (0.25d0 + 0.1372d0/(STD-1d0)**1.13d0) * Re_max**(-0.16d0)
  endif
  dP_row = chi * f_z * rho * Vmax*Vmax / 2d0
  dP = dble(Nrows) * dP_row

  write(*,'(A)') '============================================'
  write(*,'(A)') '  FLOW ACROSS BANKS OF TUBES'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  if(config==1) then
    write(*,'(A)') '  Configuration           = INLINE'
  else
    write(*,'(A)') '  Configuration           = STAGGERED'
  endif
  write(*,'(A,I6)')       '  Number of Rows N        = ',Nrows
  write(*,'(A,F10.4,A)')  '  Tube Diameter D         = ',D,' m'
  write(*,'(A,F10.4,A)')  '  Transverse Pitch S_T    = ',ST,' m'
  write(*,'(A,F10.4,A)')  '  Longitudinal Pitch S_L  = ',SL,' m'
  write(*,'(A,F10.3,A)')  '  Approach Velocity V     = ',V,' m/s'
  write(*,'(A,F10.2,A)')  '  T_inf                   = ',Tinf,' C'
  write(*,'(A,F10.2,A)')  '  T_surface               = ',Ts,' C'
  write(*,'(A)') ''
  write(*,'(A)') '--- GEOMETRY ---'
  write(*,'(A,F8.3)')     '  S_T / D                 = ',STD
  write(*,'(A,F8.3)')     '  S_L / D                 = ',SLD
  write(*,'(A,F10.3,A)')  '  Maximum Velocity V_max  = ',Vmax,' m/s'
  write(*,'(A)') ''
  write(*,'(A)') '--- RESULTS ---'
  write(*,'(A,ES14.4)')   '  Reynolds Re_max         = ',Re_max
  write(*,'(A,F10.4)')    '  C1 (Zukauskas)          = ',C1
  write(*,'(A,F10.4)')    '  m  (Zukauskas)          = ',m1
  write(*,'(A,F10.4)')    '  C2 (row correction)     = ',C2
  write(*,'(A,F12.2)')    '  Nusselt Number Nu       = ',Nu
  write(*,'(A,F12.4,A)')  '  Convection Coeff h      = ',h,' W/m2K'
  write(*,'(A,F12.2,A)')  '  Pressure Drop dP        = ',dP,' Pa'
  write(*,'(A,F12.2,A)')  '  dP per row              = ',dP_row,' Pa'
  write(*,'(A)') ''
  write(*,'(A)') '--- VELOCITY SWEEP ---'
  write(*,'(A)') '  V[m/s]   Vmax[m/s]  Re_max       Nu        h[W/m2K]  dP[Pa]'
  write(*,'(A)') '  -----------------------------------------------------------------'
  do i=1,25
    Vs = 0.5d0 + (V*4d0-0.5d0)*dble(i-1)/24d0
    if(config==1) then
      Res = rho*Vs*ST/AT*D/mu
    else
      if(AD<AT) then
        Res = rho*Vs*ST/(2d0*AD)*D/mu
      else
        Res = rho*Vs*ST/AT*D/mu
      endif
    endif
    if(config==1) then
      if(Res<100d0) then; C1=0.9d0;m1=0.4d0
      elseif(Res<1000d0) then; C1=0.52d0;m1=0.5d0
      elseif(Res<2d5) then; C1=0.27d0;m1=0.63d0
      else; C1=0.033d0;m1=0.8d0; endif
    else
      if(Res<500d0) then; C1=1.04d0;m1=0.4d0
      elseif(Res<1000d0) then; C1=0.71d0;m1=0.5d0
      elseif(Res<2d5) then; C1=0.35d0;m1=0.6d0
      else; C1=0.031d0;m1=0.8d0; endif
    endif
    Nus = C1*C2*Res**m1*Pr**0.36d0
    hs = Nus*k/D
    if(config==1) then
      dPs=dble(Nrows)*(0.25d0+0.1175d0/(STD-1d0)**1.08d0)*Res**(-0.16d0)*rho*(Vs*ST/AT)**2/2d0
    else
      dPs=dble(Nrows)*(0.25d0+0.1372d0/(STD-1d0)**1.13d0)*Res**(-0.16d0)*rho*(Vs*ST/AT)**2/2d0
    endif
    write(*,'(2X,F7.2,3X,F8.2,3X,ES10.3,2X,F9.2,2X,F10.3,2X,F10.2)') &
      Vs,Vs*ST/AT,Res,Nus,hs,dPs
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Zukauskas: Nu = C1*C2*Re_max^m * Pr^0.36 * (Pr/Pr_s)^0.25'
  write(*,'(A)') '  Ref: Incropera Ch.7 Eq.7.58, Kothandaraman Ch.8 Sec.8.6'
end program tube_banks