💻 Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Psychrometrics (Moist Air Solver)
Core Numerical Engine in Fortran 90 • 0 total downloads
psychrometrics.f90
! =========================================================================
! Source File: psychrometrics.f90
! =========================================================================
program psychrometrics
implicit none
! Constants
double precision, parameter :: R_a = 0.28705d0 ! kJ/(kg.K) gas constant for dry air
double precision, parameter :: Cp_da = 1.006d0 ! kJ/(kg.K) specific heat of dry air
double precision, parameter :: Cp_wv = 1.86d0 ! kJ/(kg.K) specific heat of water vapor
double precision, parameter :: h_fg = 2501.0d0 ! kJ/kg latent heat at 0 C
! Inputs
integer :: mode ! 1 = T & RH, 2 = T & Twb, 3 = T & Tdp, 4 = T & w
double precision :: T, P, val
! Solved properties
double precision :: RH, T_wb, T_dp, w, h, v, Pv, Psat_T
! Solver bounds
double precision :: a, b, mid, w_calc, eps
integer :: i, max_iter
logical :: has_error
! Read inputs from stdin
read(*,*,iostat=i) mode
if (i /= 0) then
write(*,*) 'ERROR: Invalid input mode.'
stop
end if
read(*,*,iostat=i) T
if (i /= 0) then
write(*,*) 'ERROR: Invalid dry-bulb temperature.'
stop
end if
read(*,*,iostat=i) P
if (i /= 0) then
write(*,*) 'ERROR: Invalid barometric pressure.'
stop
end if
read(*,*,iostat=i) val
if (i /= 0) then
write(*,*) 'ERROR: Invalid second parameter value.'
stop
end if
! Basic validations
if (P <= 0.0d0) then
write(*,*) 'ERROR: Barometric pressure must be positive.'
stop
end if
has_error = .false.
Psat_T = psat(T)
! Compute based on combination mode
select case (mode)
case (1)
! Mode 1: Dry-Bulb T & Relative Humidity RH (%)
RH = val
if (RH < 0.0d0) RH = 0.0d0
if (RH > 100.0d0) RH = 100.0d0
Pv = (RH / 100.0d0) * Psat_T
if (Pv >= P) then
Pv = P - 1.0d-4
RH = (Pv / Psat_T) * 100.0d0
end if
w = 0.62198d0 * Pv / (P - Pv)
T_dp = T_dew(Pv)
call solve_twb(T, P, w, T_dp, T_wb)
case (2)
! Mode 2: Dry-Bulb T & Wet-Bulb T_wb (C)
T_wb = val
if (T_wb > T) then
write(*,*) 'ERROR: Wet-bulb temperature cannot exceed dry-bulb temperature.'
stop
end if
w = calculate_w_from_twb(T, P, T_wb)
if (w < 0.0d0) then
w = 0.0d0
Pv = 0.0d0
RH = 0.0d0
T_dp = -100.0d0
else
Pv = P * w / (0.62198d0 + w)
if (Pv > Psat_T) then
! Clamping to saturation if solver yields physically impossible state
Pv = Psat_T
w = 0.62198d0 * Pv / (P - Pv)
RH = 100.0d0
T_dp = T
T_wb = T
else
RH = (Pv / Psat_T) * 100.0d0
T_dp = T_dew(Pv)
end if
end if
case (3)
! Mode 3: Dry-Bulb T & Dew-Point T_dp (C)
T_dp = val
if (T_dp > T) then
write(*,*) 'ERROR: Dew-point temperature cannot exceed dry-bulb temperature.'
stop
end if
Pv = psat(T_dp)
if (Pv >= P) then
write(*,*) 'ERROR: Dew-point saturation pressure exceeds barometric pressure.'
stop
end if
w = 0.62198d0 * Pv / (P - Pv)
RH = (Pv / Psat_T) * 100.0d0
call solve_twb(T, P, w, T_dp, T_wb)
case (4)
! Mode 4: Dry-Bulb T & Humidity Ratio w (kg/kg)
w = val
if (w < 0.0d0) then
write(*,*) 'ERROR: Humidity ratio cannot be negative.'
stop
end if
Pv = P * w / (0.62198d0 + w)
if (Pv > Psat_T) then
write(*,*) 'ERROR: Specified humidity ratio exceeds saturation state at T.'
stop
end if
RH = (Pv / Psat_T) * 100.0d0
T_dp = T_dew(Pv)
call solve_twb(T, P, w, T_dp, T_wb)
case default
write(*,*) 'ERROR: Unknown input combination mode.'
stop
end select
! Calculate Enthalpy (h in kJ/kg dry air)
h = Cp_da * T + w * (h_fg + Cp_wv * T)
! Calculate Specific Volume (v in m3/kg dry air)
v = R_a * (T + 273.15d0) / (P - Pv)
! Print outputs to stdout in clear key-value pairs
write(*,'(A,I2)') 'MODE=', mode
write(*,'(A,F14.4)') 'T=', T
write(*,'(A,F14.4)') 'P=', P
write(*,'(A,F14.4)') 'T_WB=', T_wb
write(*,'(A,F14.4)') 'T_DP=', T_dp
write(*,'(A,F14.4)') 'RH=', RH
write(*,'(A,F14.8)') 'W=', w
write(*,'(A,F14.4)') 'H=', h
write(*,'(A,F14.6)') 'V=', v
write(*,'(A,F14.6)') 'PV=', Pv
write(*,'(A,F14.6)') 'PSAT=', Psat_T
contains
! Saturation pressure of liquid water or ice (Buck equation) in kPa
double precision function psat(T_c)
double precision, intent(in) :: T_c
if (T_c >= 0.0d0) then
psat = 0.61121d0 * exp((17.67d0 * T_c) / (T_c + 243.5d0))
else
psat = 0.61115d0 * exp((22.452d0 * T_c) / (T_c + 272.62d0))
end if
end function psat
! Dew point calculation from partial vapor pressure (Buck equation inverse) in C
double precision function T_dew(Pv_val)
double precision, intent(in) :: Pv_val
double precision :: y
if (Pv_val <= 1.0d-8) then
T_dew = -100.0d0
return
end if
if (Pv_val >= 0.61121d0) then
y = log(Pv_val / 0.61121d0)
T_dew = (243.5d0 * y) / (17.67d0 - y)
else
y = log(Pv_val / 0.61115d0)
T_dew = (272.62d0 * y) / (22.452d0 - y)
end if
end function T_dew
! adiabatic saturation relation for humidity ratio w given dry-bulb T and wet-bulb Twb
double precision function calculate_w_from_twb(T_db, P_bar, T_wb_val)
double precision, intent(in) :: T_db, P_bar, T_wb_val
double precision :: psat_twb, w_sat_twb
psat_twb = psat(T_wb_val)
if (P_bar > psat_twb) then
w_sat_twb = 0.62198d0 * psat_twb / (P_bar - psat_twb)
else
w_sat_twb = 9999.0d0
end if
calculate_w_from_twb = ((2501.0d0 - 2.381d0 * T_wb_val) * w_sat_twb - Cp_da * (T_db - T_wb_val)) / &
(2501.0d0 + Cp_wv * T_db - 4.186d0 * T_wb_val)
end function calculate_w_from_twb
! Numerical Bisection solver for wet-bulb temperature
subroutine solve_twb(T_db, P_bar, w_act, T_dp_val, T_wb_res)
double precision, intent(in) :: T_db, P_bar, w_act, T_dp_val
double precision, intent(out) :: T_wb_res
double precision :: t_low, t_high, t_m, w_m
integer :: iter
t_low = T_dp_val
t_high = T_db
! Check if already saturated or bounds close
if (abs(t_high - t_low) < 1.0d-4) then
T_wb_res = T_db
return
end if
do iter = 1, 80
t_m = 0.5d0 * (t_low + t_high)
w_m = calculate_w_from_twb(T_db, P_bar, t_m)
if (w_m > w_act) then
t_high = t_m
else
t_low = t_m
end if
if (abs(t_high - t_low) < 1.0d-5) exit
end do
T_wb_res = 0.5d0 * (t_low + t_high)
end subroutine solve_twb
end program psychrometrics
Solver Description
Calculates psychrometric moist air properties using Dalton's law of partial pressures and Buck's high-precision saturation vapor pressure formulations. It computes humidity ratio ($\omega$), relative humidity ($RH$), dew-point ($T_{dp}$), wet-bulb ($T_{wb}$), enthalpy ($h$), specific volume ($v$), and vapor pressures ($P_v, P_{sat}$). An iterative bisection numerical method is implemented to solve for wet-bulb temperature by matching thermodynamic psychrometric wet-bulb equations.
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 -ffree-line-length-none psychrometrics.f90 -o psychro_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
psychro_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Input mode (1=T_db/RH, 2=T_db/T_wb, 3=T_db/T_dp, 4=T_db/W)
1
! Dry-Bulb Temperature ($T_{db}$) [°C]
25.0
! Barometric Pressure ($P$) [kPa]
101.325
! Second Value (RH [%] or T_wb [°C] or T_dp [°C] or W [kg/kg])
50.0
1
! Dry-Bulb Temperature ($T_{db}$) [°C]
25.0
! Barometric Pressure ($P$) [kPa]
101.325
! Second Value (RH [%] or T_wb [°C] or T_dp [°C] or W [kg/kg])
50.0