program combustion_products
    implicit none
    integer::ftype,iostat_val,i,ns_sw
    double precision::ea_pct,fuel_flow,C_f,H_f,O_f,N_f,S_f,MW_f,Hc
    double precision::O2_stoich,Air_stoich,Air_actual,O2_actual,O2_excess
    double precision::n_CO2,n_H2O,n_SO2,n_N2,n_O2,n_total
    double precision::x_CO2,x_H2O,x_SO2,x_N2,x_O2
    double precision::AFR_mass,T_air,T_ad,sum_nCp
    double precision::ea_sw,T_sw,O2_pct_sw,P1
    T_air=298.15d0
    read(*,*,iostat=iostat_val) ftype
    read(*,*,iostat=iostat_val) ea_pct
    read(*,*,iostat=iostat_val) fuel_flow
    read(*,*,iostat=iostat_val) C_f
    read(*,*,iostat=iostat_val) H_f
    read(*,*,iostat=iostat_val) O_f
    read(*,*,iostat=iostat_val) N_f
    read(*,*,iostat=iostat_val) S_f
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    select case(ftype)
    case(1);C_f=0.749d0;H_f=0.251d0;O_f=0;N_f=0;S_f=0;MW_f=16.04d0;Hc=50000.0d0
    case(2);C_f=0.817d0;H_f=0.183d0;O_f=0;N_f=0;S_f=0;MW_f=44.10d0;Hc=46350.0d0
    case(3);C_f=0.799d0;H_f=0.201d0;O_f=0;N_f=0;S_f=0;MW_f=30.07d0;Hc=47800.0d0
    case(4);C_f=0;H_f=1.0d0;O_f=0;N_f=0;S_f=0;MW_f=2.016d0;Hc=120000.0d0
    case(5);C_f=0.429d0;H_f=0;O_f=0.571d0;N_f=0;S_f=0;MW_f=28.01d0;Hc=10100.0d0
    case default;ftype=6
      if(Hc<=0) Hc=30000.0d0
      MW_f=100.0d0/(C_f/12.0d0+H_f/1.0d0+O_f/16.0d0+N_f/14.0d0+S_f/32.0d0+1.0d-30)
    end select
    if(ea_pct<0)ea_pct=0;if(fuel_flow<=0)fuel_flow=100.0d0
    O2_stoich=C_f/12.0d0+H_f/4.0d0-O_f/32.0d0+S_f/32.0d0
    Air_stoich=O2_stoich/0.21d0
    Air_actual=Air_stoich*(1.0d0+ea_pct/100.0d0)
    O2_actual=O2_stoich*(1.0d0+ea_pct/100.0d0)
    O2_excess=O2_stoich*ea_pct/100.0d0
    n_CO2=C_f/12.0d0;n_H2O=H_f/2.0d0;n_SO2=S_f/32.0d0
    n_N2=Air_actual*0.79d0+N_f/28.0d0;n_O2=O2_excess
    n_total=n_CO2+n_H2O+n_SO2+n_N2+n_O2
    if(n_total>0)then
      x_CO2=n_CO2/n_total;x_H2O=n_H2O/n_total
      x_SO2=n_SO2/n_total;x_N2=n_N2/n_total;x_O2=n_O2/n_total
    end if
    AFR_mass=Air_actual*28.97d0
    sum_nCp=n_CO2*44.2d0+n_H2O*36.0d0+n_N2*30.5d0+n_O2*32.0d0+n_SO2*50.0d0
    if(sum_nCp>0)then;T_ad=T_air+Hc*1000.0d0/sum_nCp;else;T_ad=T_air;end if
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   COMBUSTION PRODUCTS CALCULATOR'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,I4)')      '  Fuel Type                 = ',ftype
    write(*,'(A,F10.2,A)') '  Excess Air                = ',ea_pct,' percent'
    write(*,'(A,F10.2,A)') '  Fuel Flow Rate            = ',fuel_flow,' kg/h'
    write(*,'(A,F10.4)')   '  C (mass fraction)         = ',C_f
    write(*,'(A,F10.4)')   '  H (mass fraction)         = ',H_f
    write(*,'(A,F10.4)')   '  O (mass fraction)         = ',O_f
    write(*,'(A,F10.4)')   '  N (mass fraction)         = ',N_f
    write(*,'(A,F10.4)')   '  S (mass fraction)         = ',S_f
    write(*,*)
    write(*,'(A)')'--- STOICHIOMETRY -------------------------------------------'
    write(*,'(A,F12.6,A)') '  O2 Stoichiometric         = ',O2_stoich,' kmol/kg_fuel'
    write(*,'(A,F12.6,A)') '  Air Stoichiometric        = ',Air_stoich,' kmol/kg_fuel'
    write(*,'(A,F12.6,A)') '  Air Actual                = ',Air_actual,' kmol/kg_fuel'
    write(*,'(A,F12.2,A)') '  AFR (mass)                = ',AFR_mass,' kg_air/kg_fuel'
    write(*,*)
    write(*,'(A)')'--- PRODUCT COMPOSITION (per kg fuel) -----------------------'
    write(*,'(A)')'  Species   kmol/kg_fuel   mol%        mass_flow[kg/h]'
    write(*,'(A)')'  -----------------------------------------------------------'
    write(*,'(A,F12.6,2X,F8.2,2X,F12.2)') '  CO2     ',n_CO2,x_CO2*100,n_CO2*44.0d0*fuel_flow
    write(*,'(A,F12.6,2X,F8.2,2X,F12.2)') '  H2O     ',n_H2O,x_H2O*100,n_H2O*18.0d0*fuel_flow
    write(*,'(A,F12.6,2X,F8.2,2X,F12.2)') '  N2      ',n_N2,x_N2*100,n_N2*28.0d0*fuel_flow
    write(*,'(A,F12.6,2X,F8.2,2X,F12.2)') '  O2 exc  ',n_O2,x_O2*100,n_O2*32.0d0*fuel_flow
    if(n_SO2>0)write(*,'(A,F12.6,2X,F8.2,2X,F12.2)') '  SO2     ',n_SO2,x_SO2*100,n_SO2*64.0d0*fuel_flow
    write(*,*)
    write(*,'(A)')'--- ADIABATIC FLAME TEMPERATURE -----------------------------'
    write(*,'(A,F12.1,A)') '  T_adiabatic               = ',T_ad,' K'
    write(*,'(A,F12.1,A)') '  T_adiabatic               = ',T_ad-273.15d0,' C'
    write(*,*)
    ns_sw=25
    write(*,'(A)')'--- SWEEP: EXCESS AIR VS T_AD AND O2% -----------------------'
    write(*,'(A)')'  ExcessAir%  T_ad[K]     T_ad[C]     O2_dry%'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,ns_sw
      ea_sw=dble(i-1)*200.0d0/dble(ns_sw-1)
      O2_excess=O2_stoich*ea_sw/100.0d0
      n_N2=Air_stoich*(1.0d0+ea_sw/100.0d0)*0.79d0+N_f/28.0d0
      n_total=n_CO2+n_H2O+n_SO2+n_N2+O2_excess
      if(n_total-n_H2O>0)then;O2_pct_sw=O2_excess/(n_total-n_H2O)*100.0d0
      else;O2_pct_sw=0;end if
      sum_nCp=n_CO2*44.2d0+n_H2O*36.0d0+n_N2*30.5d0+O2_excess*32.0d0
      if(sum_nCp>0)then;T_sw=T_air+Hc*1000.0d0/sum_nCp;else;T_sw=T_air;end if
      write(*,'(F8.1,4X,F10.1,4X,F10.1,4X,F8.2)') ea_sw,T_sw,T_sw-273.15d0,O2_pct_sw
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  Stoichiometric O2 = C/12 + H/4 - O/32 + S/32'
    write(*,'(A)')'  T_ad = T_air + Hc/(sum ni*Cpi)'
    write(*,'(A)')'  Approximate Cp: CO2=44.2, H2O=36.0, N2=30.5 J/(mol.K)'
end program combustion_products