program air_cooled_hx
  implicit none
  integer :: i,fint,nrows,ntubes
  double precision :: Qd,Tpi,Tpo,Tai,D_t,fp,fh,pt,pl,Lt,h_proc
  double precision :: Tao,cp_air,rho_air,LMTD,Fc,dT1,dT2
  double precision :: A_bare,A_fin,A_tot,eta_f,eta_o,h_air,U,A_req
  double precision :: m_air,V_face,A_face,dP,P_fan,kfin,tfin
  double precision :: pi,Tas,Qs
  pi=3.14159265358979d0
  read(*,*) Qd; read(*,*) Tpi; read(*,*) Tpo; read(*,*) Tai
  read(*,*) fint; read(*,*) nrows; read(*,*) D_t; read(*,*) fp
  read(*,*) fh; read(*,*) pt; read(*,*) pl; read(*,*) ntubes
  read(*,*) Lt; read(*,*) h_proc
  D_t=D_t/1000d0; fh=fh/1000d0; pt=pt/1000d0; pl=pl/1000d0
  cp_air=1007d0; rho_air=1.2d0; kfin=200d0; tfin=0.0004d0
  Tao=Tai+Qd/(rho_air*5d0*ntubes*Lt*pt*cp_air)
  if(Tao>Tpi) Tao=Tpi-5d0
  m_air=Qd/(cp_air*(Tao-Tai))
  dT1=Tpi-Tao; dT2=Tpo-Tai
  if(abs(dT1-dT2)<0.01d0) then; LMTD=(dT1+dT2)/2d0
  else; LMTD=(dT1-dT2)/log(dT1/dT2); endif
  Fc=0.95d0
  A_bare=ntubes*nrows*pi*D_t*Lt
  A_fin=ntubes*nrows*Lt*fp*2d0*pi*((D_t/2d0+fh)**2-(D_t/2d0)**2)
  A_tot=A_bare*(1d0-fp*tfin)+A_fin
  if(A_tot<0.01d0) A_tot=A_bare
  if(fint==1) h_air=40d0
  if(fint==2) h_air=55d0
  if(fint==3) h_air=35d0
  eta_f=tanh(sqrt(2d0*h_air/(kfin*tfin))*fh)/(sqrt(2d0*h_air/(kfin*tfin))*fh)
  eta_o=1d0-A_fin/A_tot*(1d0-eta_f)
  U=1d0/(1d0/(eta_o*h_air)+D_t/(2d0*kfin)+1d0/h_proc)
  A_req=Qd/(U*Fc*LMTD)
  A_face=ntubes*pt*Lt
  V_face=m_air/(rho_air*A_face)
  dP=50d0*nrows
  P_fan=m_air/rho_air*dP/0.6d0
  write(*,'(A)') '============================================'
  write(*,'(A)') '  AIR-COOLED HEAT EXCHANGER (FIN-FAN)'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,F12.1,A)') '  Heat duty Q             = ',Qd,' W'
  write(*,'(A,F10.2,A)') '  Process T_in            = ',Tpi,' C'
  write(*,'(A,F10.2,A)') '  Process T_out           = ',Tpo,' C'
  write(*,'(A,F10.2,A)') '  Air T_in                = ',Tai,' C'
  if(fint==1) write(*,'(A)') '  Fin type                = Plain'
  if(fint==2) write(*,'(A)') '  Fin type                = Serrated'
  if(fint==3) write(*,'(A)') '  Fin type                = Studded'
  write(*,'(A,I4)')       '  Tube rows               = ',nrows
  write(*,'(A,I4)')       '  Tubes per row           = ',ntubes
  write(*,'(A,F10.2,A)') '  Tube length             = ',Lt,' m'
  write(*,'(A)') ''
  write(*,'(A)') '--- THERMAL RESULTS ---'
  write(*,'(A,F10.2,A)') '  Air outlet temp         = ',Tao,' C'
  write(*,'(A,F10.3,A)') '  Air mass flow           = ',m_air,' kg/s'
  write(*,'(A,F10.2,A)') '  LMTD                    = ',LMTD,' C'
  write(*,'(A,F10.4)')    '  F correction factor     = ',Fc
  write(*,'(A,F10.2,A)') '  h_air                   = ',h_air,' W/m2K'
  write(*,'(A,F10.4)')    '  Fin efficiency eta_f    = ',eta_f
  write(*,'(A,F10.4)')    '  Surface efficiency eta_o= ',eta_o
  write(*,'(A,F10.2,A)') '  Overall U               = ',U,' W/m2K'
  write(*,'(A,F12.2,A)') '  Area required           = ',A_req,' m2'
  write(*,'(A,F12.2,A)') '  Area available          = ',A_tot,' m2'
  write(*,'(A)') ''
  write(*,'(A)') '--- FAN SIZING ---'
  write(*,'(A,F10.2,A)') '  Face area               = ',A_face,' m2'
  write(*,'(A,F10.2,A)') '  Face velocity           = ',V_face,' m/s'
  write(*,'(A,F10.1,A)') '  Pressure drop est       = ',dP,' Pa'
  write(*,'(A,F10.1,A)') '  Fan power est           = ',P_fan,' W'
  write(*,'(A)') ''
  write(*,'(A)') '--- AIR TEMP SWEEP ---'
  write(*,'(A)') '  T_air[C]  LMTD[C]   A_req[m2]  P_fan[W]'
  write(*,'(A)') '  -------------------------------------------'
  do i=1,25
    Tas=10d0+30d0*dble(i-1)/24d0
    dT1=Tpi-Tas-Qd/(rho_air*5d0*ntubes*Lt*pt*cp_air)
    dT2=Tpo-Tas; if(dT2<1d0) dT2=1d0; if(dT1<1d0) dT1=1d0
    if(abs(dT1-dT2)<0.01d0) then; LMTD=(dT1+dT2)/2d0
    else; LMTD=(dT1-dT2)/log(dT1/dT2); endif
    Qs=Qd/(U*Fc*LMTD)
    write(*,'(2X,F6.1,4X,F8.2,4X,F10.2,4X,F8.1)') Tas,LMTD,Qs,P_fan
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  LMTD crossflow with F correction factor'
  write(*,'(A)') '  Fin eff: eta_f = tanh(mH)/(mH), m=sqrt(2h/(k*t))'
  write(*,'(A)') '  Overall: 1/U = 1/(eta_o*h_air) + R_wall + 1/h_proc'
  write(*,'(A)') '  Ref: API 661, ESCOA correlations, HEDH'
end program air_cooled_hx
