💻 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.
Evaporation & Humidification
Core Numerical Engine in Fortran 90 • 31 total downloads
evaporation_humidification.f90
! =========================================================================
! Source File: evaporation_humidification.f90
! =========================================================================
program evaporation_humidification
implicit none
integer :: geom_type, i, n_points, iostat_val
double precision :: L, W, Dhyd, u, Tdb, Ts, P_kPa, RH, Dab, rho, mu, kair, cp, MWA, activity, flowArea
double precision :: area, Lc, nu, alpha, Le, Re, Sc, Sh, hm, P, Rgas, Ctot
double precision :: psat_s, psat_db, pA_s, pA_inf, y_s, y_inf, omega_in, omega_s, omega_out, RH_out
double precision :: NA, mass_flux, evap_rate, mdot_air, T_wb, wb_dep, hfg, h_in, h_out_sat
double precision :: RHp, pA_in_p, yinfp, omegap, NAp, ratep, Twbp, depp, RHoutp, omegaoutp
character(len=80) :: geom_name, corr_name
read(*,*,iostat=iostat_val) geom_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid geometry type input.'
stop
end if
read(*,*,iostat=iostat_val) L
read(*,*,iostat=iostat_val) W
read(*,*,iostat=iostat_val) Dhyd
read(*,*,iostat=iostat_val) u
read(*,*,iostat=iostat_val) Tdb
read(*,*,iostat=iostat_val) Ts
read(*,*,iostat=iostat_val) P_kPa
read(*,*,iostat=iostat_val) RH
read(*,*,iostat=iostat_val) Dab
read(*,*,iostat=iostat_val) rho
read(*,*,iostat=iostat_val) mu
read(*,*,iostat=iostat_val) kair
read(*,*,iostat=iostat_val) cp
read(*,*,iostat=iostat_val) MWA
read(*,*,iostat=iostat_val) activity
read(*,*,iostat=iostat_val) flowArea
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all evaporation/humidification parameters.'
stop
end if
if (L <= 0.0d0 .or. W <= 0.0d0 .or. Dhyd <= 0.0d0) then
write(*,*) 'ERROR: Geometry dimensions must be positive.'
stop
end if
if (P_kPa <= 0.0d0 .or. Dab <= 0.0d0 .or. rho <= 0.0d0 .or. mu <= 0.0d0) then
write(*,*) 'ERROR: Pressure and fluid properties must be positive.'
stop
end if
if (RH < 0.0d0 .or. RH > 100.0d0) then
write(*,*) 'ERROR: Relative humidity must be between 0 and 100 percent.'
stop
end if
if (activity < 0.0d0) then
write(*,*) 'ERROR: Activity/Raoult factor cannot be negative.'
stop
end if
if (geom_type == 1) then
geom_name = 'Open Surface / Flat Plate'
Lc = L
area = L * W
else
geom_name = 'Humidifier Duct / Wetted Pad'
Lc = Dhyd
area = L * W
end if
P = P_kPa * 1000.0d0
Rgas = 8.314462618d0
Ctot = P / (Rgas * (Tdb + 273.15d0))
nu = mu / rho
alpha = kair / (rho * cp)
Le = alpha / Dab
Re = rho * u * Lc / mu
Sc = mu / (rho * Dab)
call calc_sh(geom_type, Re, Sc, Sh, corr_name)
hm = Sh * Dab / Lc
psat_s = saturation_pressure_water(Ts)
psat_db = saturation_pressure_water(Tdb)
pA_s = min(activity * psat_s, 0.995d0 * P)
pA_inf = min((RH/100.0d0) * psat_db, 0.995d0 * P)
y_s = pA_s / P
y_inf = pA_inf / P
omega_in = humidity_ratio(P, pA_inf)
omega_s = humidity_ratio(P, pA_s)
if (y_s > y_inf) then
NA = hm * Ctot * log((1.0d0 - y_inf)/(1.0d0 - y_s))
else
NA = hm * Ctot * (y_s - y_inf)
end if
mass_flux = NA * MWA / 1000.0d0
evap_rate = mass_flux * area
if (flowArea <= 0.0d0) flowArea = area
mdot_air = rho * u * flowArea
if (mdot_air > 0.0d0) then
omega_out = omega_in + evap_rate / mdot_air
else
omega_out = omega_in
end if
RH_out = relative_humidity_from_omega(P, Tdb, omega_out) * 100.0d0
if (RH_out > 100.0d0) RH_out = 100.0d0
T_wb = wet_bulb_temperature(Tdb, P, omega_in)
wb_dep = Tdb - T_wb
hfg = latent_heat_water(Ts)
h_in = moist_air_enthalpy(Tdb, omega_in)
h_out_sat = moist_air_enthalpy(T_wb, humidity_ratio(P, saturation_pressure_water(T_wb)))
write(*,'(A)') '============================================================'
write(*,'(A)') ' EVAPORATION AND HUMIDIFICATION ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Geometry Name = ', trim(geom_name)
write(*,'(A,A)') ' Correlation = ', trim(corr_name)
write(*,'(A,ES12.4,A)') ' Transfer Area = ', area, ' m2'
write(*,'(A,ES12.4,A)') ' Characteristic Length = ', Lc, ' m'
write(*,*)
write(*,'(A)') '--- AIR / VAPOR STATE ---------------------------------------'
write(*,'(A,F12.4,A)') ' Dry Bulb Temperature = ', Tdb, ' deg-C'
write(*,'(A,F12.4,A)') ' Surface Temperature = ', Ts, ' deg-C'
write(*,'(A,F12.4,A)') ' Total Pressure = ', P_kPa, ' kPa'
write(*,'(A,F12.4,A)') ' Inlet Relative Humidity = ', RH, ' percent'
write(*,'(A,ES12.4,A)') ' Surface Vapor Pressure = ', pA_s, ' Pa'
write(*,'(A,ES12.4,A)') ' Bulk Vapor Pressure = ', pA_inf, ' Pa'
write(*,'(A,ES12.4)') ' Surface Mole Fraction = ', y_s
write(*,'(A,ES12.4)') ' Bulk Mole Fraction = ', y_inf
write(*,'(A,ES12.4)') ' Inlet Humidity Ratio = ', omega_in
write(*,'(A,ES12.4)') ' Outlet Humidity Ratio = ', omega_out
write(*,'(A,F12.4,A)') ' Outlet Relative Humidity = ', RH_out, ' percent'
write(*,*)
write(*,'(A)') '--- DIMENSIONLESS GROUPS ------------------------------------'
write(*,'(A,ES12.4)') ' Reynolds Number (Re) = ', Re
write(*,'(A,ES12.4)') ' Schmidt Number (Sc) = ', Sc
write(*,'(A,ES12.4)') ' Sherwood Number (Sh) = ', Sh
write(*,'(A,ES12.4)') ' Lewis Number (Le) = ', Le
write(*,*)
write(*,'(A)') '--- EVAPORATION / HUMIDIFICATION RESULTS --------------------'
write(*,'(A,ES12.4,A)') ' Mass Transfer Coeff (h_m) = ', hm, ' m/s'
write(*,'(A,ES12.4,A)') ' Molar Flux (N_A) = ', NA, ' mol/m2.s'
write(*,'(A,ES12.4,A)') ' Mass Flux = ', mass_flux, ' kg/m2.s'
write(*,'(A,ES12.4,A)') ' Evaporation Rate = ', evap_rate, ' kg/s'
write(*,'(A,F12.4,A)') ' Wet Bulb Temperature = ', T_wb, ' deg-C'
write(*,'(A,F12.4,A)') ' Wet Bulb Depression = ', wb_dep, ' K'
write(*,'(A,ES12.4,A)') ' Latent Heat h_fg = ', hfg, ' J/kg'
write(*,'(A,ES12.4,A)') ' Psychrometric h_in = ', h_in, ' J/kg dry-air'
write(*,'(A,ES12.4,A)') ' Psychrometric h_sat_wb = ', h_out_sat, ' J/kg dry-air'
write(*,*)
write(*,'(A)') '--- PARAMETRIC HUMIDIFICATION PROFILE -----------------------'
write(*,'(A)') ' RH_in[%] omega_in evap_rate[kg/s] T_wb[C] depression[K] RH_out[%]'
write(*,'(A)') ' ---------------------------------------------------------------------------'
n_points = 41
do i = 0, n_points-1
RHp = 5.0d0 + 90.0d0 * dble(i) / dble(n_points-1)
pA_in_p = min((RHp/100.0d0) * psat_db, 0.995d0 * P)
yinfp = pA_in_p / P
omegap = humidity_ratio(P, pA_in_p)
if (y_s > yinfp) then
NAp = hm * Ctot * log((1.0d0 - yinfp)/(1.0d0 - y_s))
else
NAp = hm * Ctot * (y_s - yinfp)
end if
ratep = NAp * MWA / 1000.0d0 * area
Twbp = wet_bulb_temperature(Tdb, P, omegap)
depp = Tdb - Twbp
if (mdot_air > 0.0d0) then
omegaoutp = omegap + ratep / mdot_air
else
omegaoutp = omegap
end if
RHoutp = min(100.0d0, relative_humidity_from_omega(P, Tdb, omegaoutp) * 100.0d0)
write(*,'(F9.3,2X,ES12.4,2X,ES14.4,2X,F9.3,2X,F12.4,2X,F10.3)') RHp, omegap, ratep, Twbp, depp, RHoutp
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Raoult law: p_A,s = activity * p_sat(T_s).'
write(*,'(A)') ' Stefan diffusion: N_A = h_m C ln[(1-y_A,inf)/(1-y_A,s)].'
write(*,'(A)') ' Lewis number: Le = alpha/D_AB = k/(rho cp D_AB).'
write(*,'(A)') ' Wet-bulb estimated by moist-air enthalpy equality.'
contains
double precision function saturation_pressure_water(Tc)
implicit none
double precision, intent(in) :: Tc
! Buck/Magnus-style equation over liquid water, Pa.
saturation_pressure_water = 611.21d0 * exp((18.678d0 - Tc/234.5d0) * (Tc/(257.14d0 + Tc)))
end function saturation_pressure_water
double precision function humidity_ratio(Pa, pv)
implicit none
double precision, intent(in) :: Pa, pv
humidity_ratio = 0.621945d0 * pv / max(1.0d-9, Pa - pv)
end function humidity_ratio
double precision function vapor_pressure_from_omega(Pa, omega)
implicit none
double precision, intent(in) :: Pa, omega
vapor_pressure_from_omega = Pa * omega / (0.621945d0 + omega)
end function vapor_pressure_from_omega
double precision function relative_humidity_from_omega(Pa, Tc, omega)
implicit none
double precision, intent(in) :: Pa, Tc, omega
double precision :: pv, ps
pv = vapor_pressure_from_omega(Pa, omega)
ps = saturation_pressure_water(Tc)
relative_humidity_from_omega = pv / ps
end function relative_humidity_from_omega
double precision function moist_air_enthalpy(Tc, omega)
implicit none
double precision, intent(in) :: Tc, omega
moist_air_enthalpy = 1006.0d0*Tc + omega*(2501000.0d0 + 1860.0d0*Tc)
end function moist_air_enthalpy
double precision function latent_heat_water(Tc)
implicit none
double precision, intent(in) :: Tc
latent_heat_water = (2501.0d0 - 2.361d0*Tc) * 1000.0d0
end function latent_heat_water
double precision function wet_bulb_temperature(Tdry, Pa, omega)
implicit none
double precision, intent(in) :: Tdry, Pa, omega
integer :: it
double precision :: lo, hi, mid, htarget, hmid, ws
htarget = moist_air_enthalpy(Tdry, omega)
lo = -50.0d0
hi = Tdry
do it = 1, 80
mid = 0.5d0*(lo+hi)
ws = humidity_ratio(Pa, saturation_pressure_water(mid))
hmid = moist_air_enthalpy(mid, ws)
if (hmid > htarget) then
hi = mid
else
lo = mid
end if
end do
wet_bulb_temperature = 0.5d0*(lo+hi)
end function wet_bulb_temperature
subroutine calc_sh(g, Re_in, Sc_in, Sh_out, corr)
implicit none
integer, intent(in) :: g
double precision, intent(in) :: Re_in, Sc_in
double precision, intent(out) :: Sh_out
character(len=80), intent(out) :: corr
if (g == 1) then
if (Re_in < 5.0d5) then
corr = 'Flat plate laminar: Sh = 0.664 Re^0.5 Sc^(1/3)'
Sh_out = 0.664d0 * sqrt(max(Re_in,0.0d0)) * Sc_in**(1.0d0/3.0d0)
else
corr = 'Flat plate turbulent: Sh = (0.037 Re^0.8 - 871) Sc^(1/3)'
Sh_out = (0.037d0*max(Re_in,0.0d0)**0.8d0 - 871.0d0) * Sc_in**(1.0d0/3.0d0)
end if
else
if (Re_in < 2300.0d0) then
corr = 'Duct laminar approximate: Sh = 3.66'
Sh_out = 3.66d0
else
corr = 'Duct turbulent mass analogy: Sh = 0.023 Re^0.8 Sc^(1/3)'
Sh_out = 0.023d0 * max(Re_in,0.0d0)**0.8d0 * Sc_in**(1.0d0/3.0d0)
end if
end if
if (Sh_out < 0.0d0) Sh_out = 0.0d0
end subroutine calc_sh
end program evaporation_humidification
Solver Description
Model convective evaporation rates of water pools or droplets and air humidification kinetics.
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 evaporation_humidification.f90 -o evaporation_humidification
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
evaporation_humidification < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry (1=Flat Plate/Pool, 2=Circular Duct, 3=Droplet/Sphere)
1
! Length/Diameter L [m]
1.0
! Width W [m]
1.0
! Hydraulic Diameter Dhyd [m]
1.0
! Air Velocity u [m/s]
2.0
! Air Dry-Bulb Temperature Tdb [C]
25.0
! Water Surface Temperature Ts [C]
15.0
! Total Pressure P [kPa]
101.325
! Relative Humidity RH [%]
50.0
! Binary Diffusivity Dab [m2/s]
2.6e-5
! Air Density rho [kg/m3]
1.184
! Air Viscosity mu [Pa-s]
1.85e-5
! Air Thermal Conductivity k [W/m-K]
0.026
! Air Specific Heat Cp [J/kg-K]
1007.0
! Water Molecular Weight MWA [g/mol]
18.015
! Liquid Activity Coefficient
1.0
! Flow Cross-Section Area [m2]
0.5
1
! Length/Diameter L [m]
1.0
! Width W [m]
1.0
! Hydraulic Diameter Dhyd [m]
1.0
! Air Velocity u [m/s]
2.0
! Air Dry-Bulb Temperature Tdb [C]
25.0
! Water Surface Temperature Ts [C]
15.0
! Total Pressure P [kPa]
101.325
! Relative Humidity RH [%]
50.0
! Binary Diffusivity Dab [m2/s]
2.6e-5
! Air Density rho [kg/m3]
1.184
! Air Viscosity mu [Pa-s]
1.85e-5
! Air Thermal Conductivity k [W/m-K]
0.026
! Air Specific Heat Cp [J/kg-K]
1007.0
! Water Molecular Weight MWA [g/mol]
18.015
! Liquid Activity Coefficient
1.0
! Flow Cross-Section Area [m2]
0.5