💻 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.
Psychrometric Processes (HVAC)
Core Numerical Engine in Fortran 90 • 34 total downloads
psychrometric.f90
! =========================================================================
! Source File: psychrometric.f90
! =========================================================================
program psychrometric
implicit none
integer :: proc_type, iostat_val, i, n_sweep
double precision :: T_in, RH_in, P_atm, V_dot
double precision :: T_out, RH_out
double precision :: T_in2, RH_in2, V_dot2
double precision :: Ps_in, W_in, h_in, v_in, Tdp_in, Twb_in
double precision :: Ps_out, W_out, h_out, v_out, Tdp_out
double precision :: mdot_da, mdot_da2, Q_total, Q_sens, Q_lat
double precision :: m_cond, SHR, BPF, T_coil
double precision :: W_mix, h_mix, T_mix, RH_mix, Ps_mix
double precision :: frac, T_sw, Q_sw, Ps_sw, W_sw, h_sw, mdot_sw
double precision :: Pw_in, Pw_out, Pw_in2
character(len=50) :: proc_name
read(*,*,iostat=iostat_val) proc_type
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid process type.'; stop; end if
read(*,*,iostat=iostat_val) T_in
read(*,*,iostat=iostat_val) RH_in
read(*,*,iostat=iostat_val) P_atm
read(*,*,iostat=iostat_val) V_dot
read(*,*,iostat=iostat_val) T_out
read(*,*,iostat=iostat_val) RH_out
read(*,*,iostat=iostat_val) T_in2
read(*,*,iostat=iostat_val) RH_in2
read(*,*,iostat=iostat_val) V_dot2
if (iostat_val /= 0) then; write(*,*) 'ERROR: Failed to read all inputs.'; stop; end if
if (P_atm <= 0.0d0) P_atm = 101.325d0
if (V_dot <= 0.0d0) V_dot = 1.0d0
if (RH_in < 0.0d0) RH_in = 0.0d0
if (RH_in > 100.0d0) RH_in = 100.0d0
if (RH_out < 0.0d0) RH_out = 0.0d0
if (RH_out > 100.0d0) RH_out = 100.0d0
if (proc_type < 1 .or. proc_type > 5) proc_type = 1
select case(proc_type)
case(1); proc_name='Sensible Heating'
case(2); proc_name='Sensible Cooling (above dew point)'
case(3); proc_name='Cooling & Dehumidification'
case(4); proc_name='Adiabatic Humidification (Evap Cooling)'
case(5); proc_name='Mixing of Two Airstreams'
end select
! ── Psychrometric functions (inline) ──────────────────────
! P_sat (Magnus formula, kPa, T in C)
Ps_in = 0.6105d0 * exp(17.27d0*T_in/(T_in+237.3d0))
Pw_in = RH_in/100.0d0 * Ps_in
W_in = 0.622d0 * Pw_in / max(P_atm - Pw_in, 0.01d0)
h_in = 1.006d0*T_in + W_in*(2501.0d0 + 1.86d0*T_in)
v_in = 0.2871d0*(T_in+273.15d0)/(P_atm - Pw_in) * (1.0d0+1.6078d0*W_in)
! dew point (inverse Magnus)
if (Pw_in > 0.001d0) then
Tdp_in = 237.3d0*log(Pw_in/0.6105d0)/(17.27d0 - log(Pw_in/0.6105d0))
else
Tdp_in = -40.0d0
end if
! wet bulb (Stull approximation 2011)
Twb_in = T_in * atan(0.151977d0*sqrt(RH_in+8.313659d0)) + &
atan(T_in+RH_in) - atan(RH_in-1.676331d0) + &
0.00391838d0*RH_in**1.5d0*atan(0.023101d0*RH_in) - 4.686035d0
mdot_da = V_dot / max(v_in, 0.01d0)
! ── Process calculations ──────────────────────────────────
Q_total=0; Q_sens=0; Q_lat=0; m_cond=0; SHR=0; BPF=0; T_coil=0
if (proc_type == 1 .or. proc_type == 2) then
! Sensible heating/cooling: W constant
W_out = W_in
Ps_out = 0.6105d0*exp(17.27d0*T_out/(T_out+237.3d0))
if (Ps_out > 0.0d0) then
RH_out = W_out * P_atm / (0.622d0 + W_out) / Ps_out * 100.0d0
else
RH_out = 0.0d0
end if
if (RH_out > 100.0d0) RH_out = 100.0d0
h_out = 1.006d0*T_out + W_out*(2501.0d0+1.86d0*T_out)
Q_total = mdot_da * (h_out - h_in)
Q_sens = Q_total
Q_lat = 0.0d0
SHR = 1.0d0
else if (proc_type == 3) then
! Cooling & dehumidification
Ps_out = 0.6105d0*exp(17.27d0*T_out/(T_out+237.3d0))
Pw_out = RH_out/100.0d0 * Ps_out
W_out = 0.622d0*Pw_out/max(P_atm-Pw_out,0.01d0)
if (W_out > W_in) W_out = W_in ! can't add moisture in cooling
h_out = 1.006d0*T_out + W_out*(2501.0d0+1.86d0*T_out)
Q_total = mdot_da * (h_in - h_out) ! positive = heat removed
Q_sens = mdot_da * 1.006d0 * (T_in - T_out)
Q_lat = Q_total - Q_sens
if (Q_total > 0.0d0) SHR = Q_sens / Q_total
m_cond = mdot_da * max(W_in - W_out, 0.0d0)
! Coil temperature estimate (ADP)
if (RH_out > 90.0d0) then
T_coil = T_out - 2.0d0
else
T_coil = T_out - 5.0d0
end if
if (abs(T_in - T_coil) > 0.01d0) then
BPF = (T_out - T_coil) / (T_in - T_coil)
end if
else if (proc_type == 4) then
! Adiabatic humidification (evaporative cooling)
! h ≈ constant, follow wet-bulb line
! T_out ≈ T_wb + (T_in - T_wb)*(1-eta), eta~0.85
T_out = Twb_in + (T_in - Twb_in)*0.15d0 ! 85% effectiveness
h_out = h_in ! adiabatic
W_out = (h_out - 1.006d0*T_out) / (2501.0d0 + 1.86d0*T_out)
if (W_out < 0.0d0) W_out = 0.0d0
Ps_out = 0.6105d0*exp(17.27d0*T_out/(T_out+237.3d0))
Pw_out = W_out * P_atm / (0.622d0 + W_out)
if (Ps_out > 0.0d0) then
RH_out = Pw_out / Ps_out * 100.0d0
else
RH_out = 100.0d0
end if
if (RH_out > 100.0d0) RH_out = 100.0d0
Q_total = 0.0d0
Q_sens = mdot_da * 1.006d0 * (T_out - T_in)
Q_lat = -Q_sens ! equal and opposite
else if (proc_type == 5) then
! Mixing
Ps_out = 0.6105d0*exp(17.27d0*T_in2/(T_in2+237.3d0))
Pw_in2 = RH_in2/100.0d0 * Ps_out
if (V_dot2 <= 0.0d0) V_dot2 = V_dot * 0.5d0
! Stream 2 properties
W_out = 0.622d0*Pw_in2/max(P_atm - Pw_in2, 0.01d0) ! W of stream 2
v_out = 0.2871d0*(T_in2+273.15d0)/(P_atm-Pw_in2)*(1.0d0+1.6078d0*W_out)
mdot_da2 = V_dot2 / max(v_out, 0.01d0)
h_out = 1.006d0*T_in2 + W_out*(2501.0d0+1.86d0*T_in2)
! Mixed state
W_mix = (mdot_da*W_in + mdot_da2*W_out) / max(mdot_da+mdot_da2, 1.0d-10)
h_mix = (mdot_da*h_in + mdot_da2*h_out) / max(mdot_da+mdot_da2, 1.0d-10)
! T_mix from h_mix = 1.006*T + W_mix*(2501+1.86*T)
T_mix = (h_mix - W_mix*2501.0d0) / (1.006d0 + W_mix*1.86d0)
Ps_mix = 0.6105d0*exp(17.27d0*T_mix/(T_mix+237.3d0))
Pw_out = W_mix*P_atm/(0.622d0+W_mix)
if (Ps_mix > 0.0d0) then
RH_mix = Pw_out/Ps_mix*100.0d0
else
RH_mix = 0.0d0
end if
if (RH_mix > 100.0d0) RH_mix = 100.0d0
! Set "out" = mixed for display
T_out = T_mix; RH_out = RH_mix; W_out = W_mix; h_out = h_mix
Q_total = 0.0d0 ! adiabatic mixing
end if
! Outlet dew point
Pw_out = W_out * P_atm / max(0.622d0 + W_out, 1.0d-10)
if (Pw_out > 0.001d0) then
Tdp_out = 237.3d0*log(Pw_out/0.6105d0)/(17.27d0-log(Pw_out/0.6105d0))
else
Tdp_out = -40.0d0
end if
v_out = 0.2871d0*(T_out+273.15d0)/(P_atm-Pw_out)*(1.0d0+1.6078d0*W_out)
! ── Output ──────────────────────────────────────────────────
write(*,'(A)') '============================================================'
write(*,'(A)') ' PSYCHROMETRIC PROCESS ANALYSIS'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Process Type = ', trim(proc_name)
write(*,'(A,F10.2,A)') ' Atmospheric Pressure = ', P_atm, ' kPa'
write(*,'(A,F10.4,A)') ' Airflow Rate = ', V_dot, ' m3/s'
write(*,'(A,F10.4,A)') ' Dry Air Mass Flow = ', mdot_da, ' kg/s'
write(*,*)
write(*,'(A)') '--- INLET AIR STATE -----------------------------------------'
write(*,'(A,F10.2,A)') ' T_db (dry bulb) = ', T_in, ' C'
write(*,'(A,F10.2,A)') ' RH = ', RH_in, ' percent'
write(*,'(A,F10.6,A)') ' W (humidity ratio) = ', W_in, ' kg/kg_da'
write(*,'(A,F10.4,A)') ' W = ', W_in*1000.0d0, ' g/kg_da'
write(*,'(A,F10.4,A)') ' h (enthalpy) = ', h_in, ' kJ/kg_da'
write(*,'(A,F10.4,A)') ' v (specific volume) = ', v_in, ' m3/kg_da'
write(*,'(A,F10.2,A)') ' T_dp (dew point) = ', Tdp_in, ' C'
write(*,'(A,F10.2,A)') ' T_wb (wet bulb) = ', Twb_in, ' C'
write(*,'(A,F10.4,A)') ' P_sat = ', Ps_in, ' kPa'
write(*,*)
write(*,'(A)') '--- OUTLET AIR STATE ----------------------------------------'
write(*,'(A,F10.2,A)') ' T_db (dry bulb) = ', T_out, ' C'
write(*,'(A,F10.2,A)') ' RH = ', RH_out, ' percent'
write(*,'(A,F10.6,A)') ' W (humidity ratio) = ', W_out, ' kg/kg_da'
write(*,'(A,F10.4,A)') ' W = ', W_out*1000.0d0, ' g/kg_da'
write(*,'(A,F10.4,A)') ' h (enthalpy) = ', h_out, ' kJ/kg_da'
write(*,'(A,F10.4,A)') ' v (specific volume) = ', v_out, ' m3/kg_da'
write(*,'(A,F10.2,A)') ' T_dp (dew point) = ', Tdp_out, ' C'
write(*,*)
write(*,'(A)') '--- PROCESS ENERGY ANALYSIS ---------------------------------'
write(*,'(A,F12.4,A)') ' Q_total = ', Q_total, ' kW'
write(*,'(A,F12.4,A)') ' Q_sensible = ', Q_sens, ' kW'
write(*,'(A,F12.4,A)') ' Q_latent = ', Q_lat, ' kW'
write(*,'(A,F10.4)') ' SHR (sensible heat ratio) = ', SHR
if (proc_type == 3) then
write(*,'(A,F12.6,A)') ' Condensate Rate = ', m_cond, ' kg/s'
write(*,'(A,F12.4,A)') ' Condensate Rate = ', m_cond*3600.0d0, ' kg/h'
write(*,'(A,F10.4)') ' Bypass Factor BPF = ', BPF
write(*,'(A,F10.2,A)') ' Coil Temp (est. ADP) = ', T_coil, ' C'
end if
if (proc_type == 5) then
write(*,'(A,F10.2,A)') ' Stream 2 T_db = ', T_in2, ' C'
write(*,'(A,F10.2,A)') ' Stream 2 RH = ', RH_in2, ' percent'
write(*,'(A,F10.4,A)') ' Stream 2 flow = ', V_dot2, ' m3/s'
write(*,'(A,F10.4,A)') ' Stream 2 mdot_da = ', mdot_da2, ' kg/s'
end if
write(*,*)
! ── Sweep ─────────────────────────────────────────────────
n_sweep = 40
if (proc_type /= 5) then
write(*,'(A)') '--- SENSITIVITY: Q VS OUTLET TEMPERATURE --------------------'
write(*,'(A)') ' T_out[C] Q_total[kW] W_out[g/kg]'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, n_sweep
T_sw = (T_in-30.0d0) + dble(i-1)*60.0d0/dble(n_sweep-1)
Ps_sw = 0.6105d0*exp(17.27d0*T_sw/(T_sw+237.3d0))
if (proc_type <= 2) then
W_sw = W_in
else if (proc_type == 3) then
Pw_out = RH_out/100.0d0 * Ps_sw
W_sw = 0.622d0*Pw_out/max(P_atm-Pw_out,0.01d0)
if (W_sw > W_in) W_sw = W_in
else
h_sw = h_in
W_sw = (h_sw - 1.006d0*T_sw)/(2501.0d0+1.86d0*T_sw)
if (W_sw < 0.0d0) W_sw = 0.0d0
end if
h_sw = 1.006d0*T_sw + W_sw*(2501.0d0+1.86d0*T_sw)
if (proc_type == 3) then
Q_sw = mdot_da*(h_in - h_sw)
else
Q_sw = mdot_da*(h_sw - h_in)
end if
write(*,'(F10.2,4X,F12.4,4X,F10.4)') T_sw, Q_sw, W_sw*1000.0d0
end do
else
write(*,'(A)') '--- SENSITIVITY: MIXING FRACTION ----------------------------'
write(*,'(A)') ' Fraction_2 T_mix[C] RH_mix[%] W_mix[g/kg]'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, n_sweep
frac = dble(i-1)/dble(n_sweep-1)
W_sw = (1.0d0-frac)*W_in + frac*W_out
h_sw = (1.0d0-frac)*h_in + frac*h_out
T_sw = (h_sw - W_sw*2501.0d0)/(1.006d0+W_sw*1.86d0)
Ps_sw = 0.6105d0*exp(17.27d0*T_sw/(T_sw+237.3d0))
Pw_out = W_sw*P_atm/(0.622d0+W_sw)
if (Ps_sw > 0.0d0) then
mdot_sw = Pw_out/Ps_sw*100.0d0
else
mdot_sw = 0.0d0
end if
if (mdot_sw > 100.0d0) mdot_sw = 100.0d0
write(*,'(F10.4,4X,F10.2,4X,F10.2,4X,F10.4)') frac, T_sw, mdot_sw, W_sw*1000.0d0
end do
end if
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' P_sat: Magnus formula 0.6105*exp(17.27*T/(T+237.3)) [kPa]'
write(*,'(A)') ' W = 0.622*Pw/(P_atm - Pw) [kg_w/kg_da]'
write(*,'(A)') ' h = 1.006*T + W*(2501 + 1.86*T) [kJ/kg_da]'
write(*,'(A)') ' T_dp: inverse Magnus equation'
write(*,'(A)') ' T_wb: Stull (2011) approximation'
write(*,'(A)') ' SHR = Q_sensible / Q_total'
end program psychrometric
Solver Description
Model atmospheric air psychrometric properties and HVAC conditioning processes. Solves sensible heating, sensible cooling, cooling and dehumidification, evaporative cooling, and adiabatic mixing of two air streams. Computes wet-bulb temperature, dew point, humidity ratio, enthalpy, specific volume, total/sensible/latent heating or cooling load requirements, sensible heat ratio (SHR), and condensate removal rates.
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 psychrometric.f90 -o psychrometric
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
psychrometric < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Process type (1=Sensible Heating, 2=Sensible Cooling, 3=Cooling & Dehumidification, 4=Evaporative Cooling, 5=Adiabatic Mixing)
1
! Inlet dry-bulb temp [°C]
25.0
! Inlet relative humidity [%]
50.0
! Atmospheric pressure [kPa]
101.325
! Air volumetric flow rate [m3/s]
2.0
! Target outlet temperature [°C]
40.0
! Target outlet relative humidity [%]
30.0
! Stream 2 dry-bulb temp [°C]
5.0
! Stream 2 relative humidity [%]
80.0
! Stream 2 volumetric flow rate [m3/s]
1.0
1
! Inlet dry-bulb temp [°C]
25.0
! Inlet relative humidity [%]
50.0
! Atmospheric pressure [kPa]
101.325
! Air volumetric flow rate [m3/s]
2.0
! Target outlet temperature [°C]
40.0
! Target outlet relative humidity [%]
30.0
! Stream 2 dry-bulb temp [°C]
5.0
! Stream 2 relative humidity [%]
80.0
! Stream 2 volumetric flow rate [m3/s]
1.0