🔬
Solver Purpose & Physical Scope
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.
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 455 times
📄 Source File:
psychrometric.f90
📁 calcul/Thermodynamics /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 psychrometric.f90 -o psychrometric
2. Execution with input.txt redirection:
psychrometric < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 25.0 50.0 101.325 2.0 40.0 30.0 5.0 80.0 1.0
Parameter Description:
Process type (1=Sensible Heating, 2=Sensible Cooling, 3=Cooling & Dehumidification, 4=Evaporative Cooling, 5=Adiabatic Mixing) Inlet dry-bulb temp [°C] Inlet relative humidity [%] Atmospheric pressure [kPa] Air volumetric flow rate [m3/s] Target outlet temperature [°C] Target outlet relative humidity [%] Stream 2 dry-bulb temp [°C] Stream 2 relative humidity [%] Stream 2 volumetric flow rate [m3/s]