💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Gas Radiation (Leckner)
Core Numerical Engine in Fortran 90 • 30 total downloads
gas_radiation.f90
! =========================================================================
! Source File: gas_radiation.f90
! =========================================================================
program gas_radiation
implicit none
! Inputs
double precision :: Tg_C, Ts_C
double precision :: pc, pw, L, Pt
! Internal temperatures
double precision :: Tg_K, Ts_K
double precision :: tg, ts
! Leckner matrices
double precision :: c(0:3, 0:2), d(0:3, 0:2)
! Emissivity variables
double precision :: L_cm, pL_c, pL_w, pL_sum
double precision :: x_g, y_c, y_w
double precision :: eps_0c, eps_0w
double precision :: Pec, Pew, Cc, Cw
double precision :: eps_c, eps_w
double precision :: zeta, overlap_eps, eps_total
! Absorptivity variables
double precision :: pL_c_s, pL_w_s, pL_sum_s
double precision :: x_s, y_c_s, y_w_s
double precision :: eps_0c_s, eps_0w_s
double precision :: alpha_c, alpha_w, overlap_alpha, alpha_total, temp_term
integer :: i, j, iostat_val
! Initialize Leckner matrices
! CO2 Emissivity Matrix (c_ij)
c(0,0) = -3.5430d0; c(0,1) = 1.6370d0; c(0,2) = -0.5050d0
c(1,0) = 0.3800d0; c(1,1) = -0.2760d0; c(1,2) = 0.1060d0
c(2,0) = -0.1180d0; c(2,1) = 0.0540d0; c(2,2) = -0.0160d0
c(3,0) = 0.0170d0; c(3,1) = -0.0070d0; c(3,2) = 0.0020d0
! H2O Emissivity Matrix (d_ij)
d(0,0) = -2.2118d0; d(0,1) = -1.1987d0; d(0,2) = 0.035596d0
d(1,0) = 0.85667d0; d(1,1) = 0.93048d0; d(1,2) = -0.14391d0
d(2,0) = -0.10833d0; d(2,1) = -0.17156d0; d(2,2) = 0.045915d0
d(3,0) = -0.016333d0; d(3,1) = 0.01988d0; d(3,2) = -0.005830d0
! Read inputs
read(*,*,iostat=iostat_val) Tg_C
read(*,*,iostat=iostat_val) Ts_C
read(*,*,iostat=iostat_val) pc
read(*,*,iostat=iostat_val) pw
read(*,*,iostat=iostat_val) L
read(*,*,iostat=iostat_val) Pt
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read gas calculation parameters.'
stop
end if
! Basic validations
Tg_K = Tg_C + 273.15d0
Ts_K = Ts_C + 273.15d0
if (Tg_K <= 0.0d0 .or. Ts_K <= 0.0d0) then
write(*,*) 'ERROR: Temperatures must be above absolute zero.'
stop
end if
if (pc < 0.0d0 .or. pw < 0.0d0 .or. Pt <= 0.0d0) then
write(*,*) 'ERROR: Gas pressures must be positive.'
stop
end if
if ((pc + pw) > Pt) then
write(*,*) 'ERROR: Sum of partial pressures cannot exceed total pressure.'
stop
end if
if (L <= 0.0d0) then
write(*,*) 'ERROR: Mean beam length must be positive.'
stop
end if
! Units conversion: path length to cm
L_cm = L * 100.0d0
pL_c = pc * L_cm
pL_w = pw * L_cm
pL_sum = (pc + pw) * L_cm
tg = Tg_K / 1000.0d0
x_g = log(tg)
! =======================================================
! 1. EMISSIVITY CALCULATIONS (at gas temperature Tg)
! =======================================================
! Base Emissivity CO2
y_c = log(max(1.0001d0, pL_c))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + c(i, j) * (x_g**i) * (y_c**j)
end do
end do
eps_0c = exp(temp_term)
! Base Emissivity H2O
y_w = log(max(1.0001d0, pL_w))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + d(i, j) * (x_g**i) * (y_w**j)
end do
end do
eps_0w = exp(temp_term)
! Pressure corrections
Pec = Pt + 0.28d0 * pc
Cc = 1.0d0 + ((Pec - 1.0d0) / (Pec - 1.0d0 + 0.28d0)) * 0.15d0 * &
(1.0d0 - tanh(log10(max(1.0001d0, pL_c))))
Pew = Pt + 0.85d0 * pw
Cw = 1.0d0 + ((Pew - 1.0d0) / (Pew - 1.0d0 + 0.5d0)) * 0.35d0 * &
(1.0d0 - tanh(log10(max(1.0001d0, pL_w))))
eps_c = eps_0c * Cc
eps_w = eps_0w * Cw
! Overlap correction (zeta)
overlap_eps = 0.0d0
if ((pc + pw) > 1.0d-8) then
zeta = pw / (pc + pw)
if (pL_sum >= 1.0d0) then
overlap_eps = zeta * (10.7d0 / (10.7d0 + 101.0d0 * zeta) - 0.0089d0 * zeta**10.4d0) * &
(log10(pL_sum))**2.76d0
end if
end if
eps_total = eps_c + eps_w - overlap_eps
if (eps_total < 0.0d0) eps_total = 0.0d0
! =======================================================
! 2. ABSORPTIVITY CALCULATIONS (for surface at Ts)
! =======================================================
ts = Ts_K / 1000.0d0
x_s = log(ts)
! CO2 Absorptivity (evaluated at Ts and scaled pL)
pL_c_s = pc * L_cm * (Ts_K / Tg_K)
y_c_s = log(max(1.0001d0, pL_c_s))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + c(i, j) * (x_s**i) * (y_c_s**j)
end do
end do
eps_0c_s = exp(temp_term)
alpha_c = Cc * eps_0c_s * (Tg_K / Ts_K)**0.65d0
! H2O Absorptivity
pL_w_s = pw * L_cm * (Ts_K / Tg_K)
y_w_s = log(max(1.0001d0, pL_w_s))
temp_term = 0.0d0
do i = 0, 3
do j = 0, 2
temp_term = temp_term + d(i, j) * (x_s**i) * (y_w_s**j)
end do
end do
eps_0w_s = exp(temp_term)
alpha_w = Cw * eps_0w_s * (Tg_K / Ts_K)**0.45d0
! Overlap correction for absorptivity
overlap_alpha = 0.0d0
pL_sum_s = (pc + pw) * L_cm * (Ts_K / Tg_K)
if ((pc + pw) > 1.0d-8) then
zeta = pw / (pc + pw)
if (pL_sum_s >= 1.0d0) then
overlap_alpha = zeta * (10.7d0 / (10.7d0 + 101.0d0 * zeta) - 0.0089d0 * zeta**10.4d0) * &
(log10(pL_sum_s))**2.76d0
end if
end if
alpha_total = alpha_c + alpha_w - overlap_alpha
if (alpha_total < 0.0d0) alpha_total = 0.0d0
! =======================================================
! OUTPUT DATA
! =======================================================
write(*,*) '============================================================'
write(*,*) ' LECKNER MIXTURE GAS EMISSION ENGINE'
write(*,*) '============================================================'
write(*,*)
write(*,'(A,F12.2,A)') ' Gas Temperature (Tg) = ', Tg_C, ' deg-C'
write(*,'(A,F12.2,A)') ' Wall Temperature (Ts) = ', Ts_C, ' deg-C'
write(*,'(A,F12.4,A)') ' Path Length (L) = ', L, ' m'
write(*,'(A,F12.4,A)') ' Total Pressure (Pt) = ', Pt, ' bar'
write(*,'(A,F12.4,A)') ' CO2 Partial Pressure = ', pc, ' bar'
write(*,'(A,F12.4,A)') ' H2O Partial Pressure = ', pw, ' bar'
write(*,*)
write(*,*) '--- GAS EMISSIVITY PROPERTIES ------------------------------'
write(*,'(A,F12.6)') ' CO2 Base Emissivity = ', eps_0c
write(*,'(A,F12.6)') ' CO2 Pressure Coeff Cc = ', Cc
write(*,'(A,F12.6)') ' CO2 Corrected Emiss = ', eps_c
write(*,'(A,F12.6)') ' H2O Base Emissivity = ', eps_0w
write(*,'(A,F12.6)') ' H2O Pressure Coeff Cw = ', Cw
write(*,'(A,F12.6)') ' H2O Corrected Emiss = ', eps_w
write(*,'(A,F12.6)') ' Mixture Overlap Corr = ', overlap_eps
write(*,'(A,F12.6)') ' Total Gas Emissivity = ', eps_total
write(*,*)
write(*,*) '--- GAS ABSORPTIVITY PROPERTIES ----------------------------'
write(*,'(A,F12.6)') ' CO2 Absorptivity = ', alpha_c
write(*,'(A,F12.6)') ' H2O Absorptivity = ', alpha_w
write(*,'(A,F12.6)') ' Mixture Overlap Corr = ', overlap_alpha
write(*,'(A,F12.6)') ' Total Gas Absorptivity = ', alpha_total
write(*,*)
end program gas_radiation
Solver Description
Participating gases (such as $CO_2$ and $H_2O$) emit and absorb radiation in discrete spectral bands. Unlike solids, they do not act as black or gray bodies. The total gas mixture emissivity ($\varepsilon_g$) is modeled using Leckner's bivariate correlations (1972) as a function of temperature ($T_g$) and partial pressure-length parameter ($p_g L_e$):
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 gas_radiation.f90 -o gas_radiation
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
gas_radiation < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Gas temperature Tg [°C]
1000.0
! Wall temperature Ts [°C]
300.0
! CO2 partial pressure pc [bar]
0.10
! H2O partial pressure pw [bar]
0.15
! Mean beam length L [m]
2.0
! Total pressure Pt [bar]
1.0
1000.0
! Wall temperature Ts [°C]
300.0
! CO2 partial pressure pc [bar]
0.10
! H2O partial pressure pw [bar]
0.15
! Mean beam length L [m]
2.0
! Total pressure Pt [bar]
1.0