program cooling_tower
  implicit none
  integer :: i,ft
  double precision :: Twi,Two,Twb,mw,LG,Hf,Pfan
  double precision :: Rng,App,dT,eff,Qrej,cp,hfg
  double precision :: Me,ma,m_evap,m_makeup
  double precision :: T,Pvs,Ws,hs,ha,Wa,Pva
  double precision :: h1,h2,h3,h4,h5,T1,T2,T3,T4,T5
  double precision :: f1,f2,f3,f4,f5,COP,LGs,Mes,effs
  read(*,*) Twi; read(*,*) Two; read(*,*) Twb
  read(*,*) mw; read(*,*) LG; read(*,*) Hf
  read(*,*) ft; read(*,*) Pfan
  cp=4186d0; hfg=2450d3
  Rng=Twi-Two; App=Two-Twb
  Qrej=mw*cp*Rng
  eff=Rng/(Twi-Twb)
  ma=mw/LG
  Pva=610.78d0*exp(17.27d0*Twb/(Twb+237.3d0))
  Wa=0.622d0*Pva/(101325d0-Pva)
  ha=1.006d0*Twb+Wa*(2501d0+1.86d0*Twb)
  T1=Two; T2=Two+Rng*0.25d0; T3=Two+Rng*0.5d0
  T4=Two+Rng*0.75d0; T5=Twi
  Pvs=610.78d0*exp(17.27d0*T1/(T1+237.3d0))
  Ws=0.622d0*Pvs/(101325d0-Pvs)
  h1=1.006d0*T1+Ws*(2501d0+1.86d0*T1)
  Pvs=610.78d0*exp(17.27d0*T2/(T2+237.3d0))
  Ws=0.622d0*Pvs/(101325d0-Pvs)
  h2=1.006d0*T2+Ws*(2501d0+1.86d0*T2)
  Pvs=610.78d0*exp(17.27d0*T3/(T3+237.3d0))
  Ws=0.622d0*Pvs/(101325d0-Pvs)
  h3=1.006d0*T3+Ws*(2501d0+1.86d0*T3)
  Pvs=610.78d0*exp(17.27d0*T4/(T4+237.3d0))
  Ws=0.622d0*Pvs/(101325d0-Pvs)
  h4=1.006d0*T4+Ws*(2501d0+1.86d0*T4)
  Pvs=610.78d0*exp(17.27d0*T5/(T5+237.3d0))
  Ws=0.622d0*Pvs/(101325d0-Pvs)
  h5=1.006d0*T5+Ws*(2501d0+1.86d0*T5)
  f1=1d0/(h1-ha)
  f2=1d0/(h2-(ha+0.25d0*cp*Rng/(ma*1000d0)))
  f3=1d0/(h3-(ha+0.5d0*cp*Rng/(ma*1000d0)))
  f4=1d0/(h4-(ha+0.75d0*cp*Rng/(ma*1000d0)))
  f5=1d0/(h5-(ha+cp*Rng/(ma*1000d0)))
  if(f1<0d0) f1=abs(f1); if(f2<0d0) f2=abs(f2)
  if(f3<0d0) f3=abs(f3); if(f4<0d0) f4=abs(f4)
  if(f5<0d0) f5=abs(f5)
  dT=Rng/4d0
  Me=dT/3d0*(f1+4d0*f2+2d0*f3+4d0*f4+f5)
  m_evap=Qrej/hfg
  m_makeup=m_evap*1.3d0
  if(Pfan>0d0) then; COP=Qrej/(Pfan*1000d0); else; COP=0d0; endif
  write(*,'(A)') '============================================'
  write(*,'(A)') '  COOLING TOWER DESIGN (MERKEL METHOD)'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,F10.2,A)') '  Water inlet temp Tw_in  = ',Twi,' C'
  write(*,'(A,F10.2,A)') '  Water outlet temp Tw_out= ',Two,' C'
  write(*,'(A,F10.2,A)') '  Wet-bulb temp T_wb      = ',Twb,' C'
  write(*,'(A,F10.3,A)') '  Water flow rate m_w     = ',mw,' kg/s'
  write(*,'(A,F10.3)')    '  L/G ratio               = ',LG
  write(*,'(A,F10.2,A)') '  Fill height             = ',Hf,' m'
  if(ft==1) write(*,'(A)') '  Fill type               = Splash'
  if(ft==2) write(*,'(A)') '  Fill type               = Film'
  if(ft==3) write(*,'(A)') '  Fill type               = Trickle'
  write(*,'(A,F10.2,A)') '  Fan power               = ',Pfan,' kW'
  write(*,'(A)') ''
  write(*,'(A)') '--- THERMAL RESULTS ---'
  write(*,'(A,F10.2,A)') '  Range                   = ',Rng,' C'
  write(*,'(A,F10.2,A)') '  Approach                = ',App,' C'
  write(*,'(A,F10.4)')    '  Merkel Number KaV/L     = ',Me
  write(*,'(A,F10.4)')    '  Effectiveness           = ',eff
  write(*,'(A,F12.1,A)') '  Heat rejected Q         = ',Qrej,' W'
  write(*,'(A,F12.1,A)') '  Heat rejected Q         = ',Qrej/1000d0,' kW'
  write(*,'(A)') ''
  write(*,'(A)') '--- AIR SIDE ---'
  write(*,'(A,F10.3,A)') '  Air mass flow m_a       = ',ma,' kg/s'
  write(*,'(A,F10.4,A)') '  Air inlet humidity W    = ',Wa,' kg/kg'
  write(*,'(A,F10.2,A)') '  Air inlet enthalpy h_a  = ',ha,' kJ/kg'
  write(*,'(A)') ''
  write(*,'(A)') '--- WATER BALANCE ---'
  write(*,'(A,F10.4,A)') '  Evaporation rate        = ',m_evap,' kg/s'
  write(*,'(A,F10.2,A)') '  Evaporation pct         = ',m_evap/mw*100d0,' %'
  write(*,'(A,F10.4,A)') '  Makeup water            = ',m_makeup,' kg/s'
  if(COP>0d0) write(*,'(A,F10.1)') '  COP (Q/Pfan)            = ',COP
  write(*,'(A)') ''
  write(*,'(A)') '--- L/G RATIO SWEEP ---'
  write(*,'(A)') '  L/G      Me        eff       Q_rej[kW]'
  write(*,'(A)') '  -------------------------------------------'
  do i=1,25
    LGs=0.5d0+2.5d0*dble(i-1)/24d0
    effs=Rng/(Twi-Twb)
    Mes=Me*LG/LGs
    write(*,'(2X,F6.2,2X,F10.4,2X,F8.4,2X,F10.1)') LGs,Mes,effs,Qrej/1000d0
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Merkel: integral dT/(h_s - h_a) via Simpson rule'
  write(*,'(A)') '  h_s = 1.006*T + W_s*(2501+1.86*T) [kJ/kg]'
  write(*,'(A)') '  Pvs = 610.78*exp(17.27*T/(T+237.3)) [Pa]'
  write(*,'(A)') '  W_s = 0.622*Pvs/(Patm-Pvs)'
  write(*,'(A)') '  Ref: Merkel (1925), Kloppers & Kroger (2005)'
  write(*,'(A)') '       ASHRAE Handbook, CTI Standards'
end program cooling_tower