program evaporator_design
  implicit none
  integer :: i,ne
  double precision :: mf,xf,xp,Ts,BPEc,U0,Tc
  double precision :: m_evap,m_prod,hfg,Q_total,A_total
  double precision :: dT_avail,dT_eff,BPE_tot,economy
  double precision :: m_steam,spec_en,Q_e,A_e,x_e,T_e,dT_e
  double precision :: eco_s,A_s,Q_s
  read(*,*) mf; read(*,*) xf; read(*,*) xp; read(*,*) Ts
  read(*,*) ne; read(*,*) BPEc; read(*,*) U0; read(*,*) Tc
  if(BPEc<1d-6) BPEc=1.78d0
  if(U0<1d-6) U0=2500d0
  hfg=2260d3
  m_prod=mf*xf/xp
  m_evap=mf-m_prod
  BPE_tot=BPEc*(xf+xp)/2d0*dble(ne)
  dT_avail=Ts-Tc-BPE_tot
  if(dT_avail<1d0) dT_avail=1d0
  dT_eff=dT_avail/dble(ne)
  Q_total=m_evap*hfg
  A_total=0d0
  m_steam=m_evap/dble(ne)
  economy=m_evap/m_steam
  spec_en=m_steam*hfg/m_evap
  do i=1,ne
    Q_e=m_evap/dble(ne)*hfg
    A_e=Q_e/(U0*dT_eff)
    A_total=A_total+A_e
  enddo
  write(*,'(A)') '============================================'
  write(*,'(A)') '  MULTI-EFFECT EVAPORATOR DESIGN'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,F10.4,A)') '  Feed flow rate          = ',mf,' kg/s'
  write(*,'(A,F10.4)')    '  Feed concentration      = ',xf
  write(*,'(A,F10.4)')    '  Product concentration   = ',xp
  write(*,'(A,F10.2,A)') '  Steam temperature       = ',Ts,' C'
  write(*,'(A,I4)')       '  Number of effects       = ',ne
  write(*,'(A,F10.4)')    '  BPE coefficient         = ',BPEc
  write(*,'(A,F10.1,A)') '  Overall U               = ',U0,' W/m2K'
  write(*,'(A,F10.2,A)') '  Condenser temp          = ',Tc,' C'
  write(*,'(A)') ''
  write(*,'(A)') '--- RESULTS ---'
  write(*,'(A,F10.4,A)') '  Product flow            = ',m_prod,' kg/s'
  write(*,'(A,F10.4,A)') '  Total evaporation       = ',m_evap,' kg/s'
  write(*,'(A,F10.4,A)') '  Steam consumption       = ',m_steam,' kg/s'
  write(*,'(A,F10.2)')    '  Steam economy           = ',economy
  write(*,'(A,F12.1,A)') '  Total heat duty Q       = ',Q_total,' W'
  write(*,'(A,F12.1,A)') '  Total heat duty Q       = ',Q_total/1d6,' MW'
  write(*,'(A,F12.2,A)') '  Total area required     = ',A_total,' m2'
  write(*,'(A,F10.2,A)') '  dT per effect           = ',dT_eff,' C'
  write(*,'(A,F10.2,A)') '  BPE total               = ',BPE_tot,' C'
  write(*,'(A,F10.1,A)') '  Specific energy         = ',spec_en/1d3,' kJ/kg'
  write(*,'(A)') ''
  write(*,'(A)') '--- EFFECT-BY-EFFECT ---'
  write(*,'(A)') '  Effect  Q[kW]      A[m2]      x_out      T_boil[C]'
  write(*,'(A)') '  -------------------------------------------------------'
  do i=1,ne
    Q_e=m_evap/dble(ne)*hfg
    A_e=Q_e/(U0*dT_eff)
    x_e=xf+(xp-xf)*dble(i)/dble(ne)
    T_e=Ts-dT_eff*dble(i)
    write(*,'(2X,I4,2X,F10.1,2X,F10.2,2X,F10.4,2X,F10.2)') i,Q_e/1000d0,A_e,x_e,T_e
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- EFFECTS SWEEP (1 to 6) ---'
  write(*,'(A)') '  N_eff  Economy   A_total[m2]  Spec.En[kJ/kg]'
  write(*,'(A)') '  -----------------------------------------------'
  do i=1,6
    eco_s=dble(i)
    Q_s=m_evap*hfg
    dT_e=(Ts-Tc-BPEc*(xf+xp)/2d0*dble(i))/dble(i)
    if(dT_e<0.1d0) dT_e=0.1d0
    A_s=dble(i)*Q_s/(dble(i)*U0*dT_e)
    write(*,'(2X,I4,2X,F8.2,4X,F10.2,4X,F10.1)') i,eco_s,A_s,hfg/eco_s/1000d0
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Economy = m_evap_total / m_steam (approaches N for ideal)'
  write(*,'(A)') '  BPE = BPE_coeff * concentration (linearized Duhring)'
  write(*,'(A)') '  Q = U * A * dT_eff per effect'
  write(*,'(A)') '  Ref: McCabe, Smith & Harriott, Unit Operations Ch.16'
  write(*,'(A)') '       Geankoplis, Transport Processes, Ch.8'
end program evaporator_design