💻 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.
Steam Tables Advanced (IAPWS-IF97)
Core Numerical Engine in Fortran 90 • 45 total downloads
iapws_if97.f90
! =========================================================================
! Source File: iapws_if97.f90
! =========================================================================
! ============================================================================
! ThermoFluidCalc — IAPWS-IF97 Water and Steam Properties Calculator
! Self-contained Fortran 90 implementation generated programmatically
! ============================================================================
module iapws_if97_mod
implicit none
double precision, parameter :: R = 0.461526D0 ! kJ/(kg K)
double precision, parameter :: Tc = 647.096D0 ! K
double precision, parameter :: Pc = 22.064D0 ! MPa
double precision, parameter :: rhoc = 322.0D0 ! kg/m3
double precision, parameter :: Pmin = 0.000611212677444D0 ! MPa
double precision, parameter :: Ps_623 = 16.5291642526D0 ! MPa
! Region 1 coefficients
double precision, dimension(34) :: r1_n = (/ &
1.463297121316700D-01, -8.454818716911400D-01, -3.756360367204000D+00, &
3.385516916838500D+00, -9.579196338787200D-01, 1.577203851322800D-01, &
-1.661641719950100D-02, 8.121462998356800D-04, 2.831908012380400D-04, &
-6.070630156587400D-04, -1.899006821841900D-02, -3.252974877050500D-02, &
-2.184171717541400D-02, -5.283835796993000D-05, -4.718432107326700D-04, &
-3.000178079302600D-04, 4.766139390698700D-05, -4.414184533084600D-06, &
-7.269499629759400D-16, -3.167964484505400D-05, -2.827079798531200D-06, &
-8.520512812010300D-10, -2.242528190800000D-06, -6.517122289560100D-07, &
-1.434172993792400D-13, -4.051699686011700D-07, -1.273430174164100D-09, &
-1.742487123063400D-10, -6.876213129553100D-19, 1.447830782852100D-20, &
2.633578166279500D-23, -1.194762264007100D-23, 1.822809458140400D-24, &
-9.353708729245800D-26 /)
integer, dimension(34) :: r1_li = (/ &
0, 0, 0, 0, 0, 0, 0, 0, &
1, 1, 1, 1, 1, 1, 2, 2, &
2, 2, 2, 3, 3, 3, 4, 4, &
4, 5, 8, 8, 21, 23, 29, 30, &
31, 32 /)
integer, dimension(34) :: r1_lj = (/ &
-2, -1, 0, 1, 2, 3, 4, 5, &
-9, -7, -1, 0, 1, 3, -3, 0, &
1, 3, 17, -4, 0, 6, -5, -2, &
10, -8, -11, -6, -29, -31, -38, -39, &
-40, -41 /)
! Region 2 coefficients
double precision, dimension(43) :: r2_n = (/ &
-1.773174247321300D-03, -1.783486229235800D-02, -4.599601369636500D-02, &
-5.758125908343200D-02, -5.032527872793000D-02, -3.303264167020300D-05, &
-1.894898751631500D-04, -3.939277724335500D-03, -4.379729565057300D-02, &
-2.667454791408700D-05, 2.048173769230900D-08, 4.387066728443500D-07, &
-3.227767723857000D-05, -1.503392454214800D-03, -4.066825356264900D-02, &
-7.884730955936700D-10, 1.279071785228500D-08, 4.822537271850700D-07, &
2.292207633766100D-06, -1.671476645106100D-11, -2.117147232135500D-03, &
-2.389574193410400D+01, -5.905956432427000D-18, -1.262180889910100D-06, &
-3.894684243573900D-02, 1.125621136045900D-11, -8.231134089799800D+00, &
1.980971280208800D-08, 1.040696521017400D-19, -1.023474709592900D-13, &
-1.001817937951100D-09, -8.088290864698500D-11, 1.069303187940900D-01, &
-3.366225057417100D-01, 8.918584535542100D-25, 3.062931687623200D-13, &
-4.200246769820800D-06, -5.905602968563900D-26, 3.782694761345700D-06, &
-1.276860893468100D-15, 7.308761059506100D-29, 5.541471535077800D-17, &
-9.436970724121000D-07 /)
integer, dimension(43) :: r2_li = (/ &
1, 1, 1, 1, 1, 2, 2, 2, &
2, 2, 3, 3, 3, 3, 3, 4, &
4, 4, 5, 6, 6, 6, 7, 7, &
7, 8, 8, 9, 10, 10, 10, 16, &
16, 18, 20, 20, 20, 21, 22, 23, &
24, 24, 24 /)
integer, dimension(43) :: r2_lj = (/ &
0, 1, 2, 3, 6, 1, 2, 4, &
7, 36, 0, 1, 3, 6, 35, 1, &
2, 3, 7, 3, 16, 35, 0, 11, &
25, 8, 36, 13, 4, 10, 14, 29, &
50, 57, 20, 35, 48, 21, 53, 39, &
26, 40, 58 /)
! Region 2 Ideal gas cp0 coefficients
double precision, dimension(9) :: r2_cp0_no = (/ &
-9.692768650021700D+00, 1.008665596801800D+01, -5.608791128302000D-03, &
7.145273808145500D-02, -4.071049822392800D-01, 1.424081917144400D+00, &
-4.383951131945000D+00, -2.840863246077200D-01, 2.126846375330700D-02 /)
integer, dimension(9) :: r2_cp0_jo = (/ &
0, 1, -5, -4, -3, -2, -1, 2, &
3 /)
! Region 3 coefficients
double precision, dimension(39) :: r3_n = (/ &
-1.573284529023900D+01, 2.094439697430700D+01, -7.686770787871600D+00, &
2.618594778795400D+00, -2.808078114862000D+00, 1.205336969651700D+00, &
-8.456681281250201D-03, -1.265431547771400D+00, -1.152440780668100D+00, &
8.852104398431800D-01, -6.420776518160700D-01, 3.849346018667100D-01, &
-8.521470882420600D-01, 4.897228154187700D+00, -3.050261725696500D+00, &
3.942053687915400D-02, 1.255840842430800D-01, -2.799932969871000D-01, &
1.389979956946000D+00, -2.018991502357000D+00, -8.214763717396300D-03, &
-4.759603573492300D-01, 4.398407447350000D-02, -4.447643542873900D-01, &
9.057207071973300D-01, 7.052245008796700D-01, 1.077051262633200D-01, &
-3.291362325895400D-01, -5.087106204115800D-01, -2.217540087309600D-02, &
9.426075166509200D-02, 1.643627844796100D-01, -1.350337224134800D-02, &
-1.483434535247200D-02, 5.792295362808399D-04, 3.230890470371100D-03, &
8.096480299621501D-05, -1.655767979503700D-04, -4.492389906181500D-05 /)
integer, dimension(39) :: r3_li = (/ &
0, 0, 0, 0, 0, 0, 0, 1, &
1, 1, 1, 2, 2, 2, 2, 2, &
2, 3, 3, 3, 3, 3, 4, 4, &
4, 4, 5, 5, 5, 6, 6, 6, &
7, 8, 9, 9, 10, 10, 11 /)
integer, dimension(39) :: r3_lj = (/ &
0, 1, 2, 7, 10, 12, 23, 2, &
6, 15, 17, 0, 2, 6, 7, 22, &
26, 0, 2, 4, 16, 26, 0, 2, &
4, 26, 1, 3, 26, 0, 2, 26, &
2, 26, 2, 26, 0, 1, 26 /)
contains
! -------------------------------------------------------------------------
! Saturation Pressure: P = f(T) in MPa
! -------------------------------------------------------------------------
double precision function psat(T)
double precision, intent(in) :: T
double precision :: n(0:10)
double precision :: tita, A, B, C
n(0) = 0.0D0
n(1) = 0.11670521452767D+04
n(2) = -0.72421316703206D+06
n(3) = -0.17073846940092D+02
n(4) = 0.12020824702470D+05
n(5) = -0.32325550322333D+07
n(6) = 0.14915108613530D+02
n(7) = -0.48232657361591D+04
n(8) = 0.40511340542057D+06
n(9) = -0.23855557567849D+00
n(10) = 0.65017534844798D+03
if (T < 273.15D0 .or. T > Tc) then
psat = -1.0D0
return
end if
tita = T + n(9) / (T - n(10))
A = tita**2 + n(1)*tita + n(2)
B = n(3)*tita**2 + n(4)*tita + n(5)
C = n(6)*tita**2 + n(7)*tita + n(8)
psat = (2.0D0 * C / (-B + sqrt(B**2 - 4.0D0*A*C)))**4
end function psat
! -------------------------------------------------------------------------
! Saturation Temperature: T = f(P) in K
! -------------------------------------------------------------------------
double precision function tsat(P)
double precision, intent(in) :: P
double precision :: n(0:10)
double precision :: beta, E, F, G, D
n(0) = 0.0D0
n(1) = 0.11670521452767D+04
n(2) = -0.72421316703206D+06
n(3) = -0.17073846940092D+02
n(4) = 0.12020824702470D+05
n(5) = -0.32325550322333D+07
n(6) = 0.14915108613530D+02
n(7) = -0.48232657361591D+04
n(8) = 0.40511340542057D+06
n(9) = -0.23855557567849D+00
n(10) = 0.65017534844798D+03
if (P < Pmin .or. P > Pc) then
tsat = -1.0D0
return
end if
beta = P**0.25D0
E = beta**2 + n(3)*beta + n(6)
F = n(1)*beta**2 + n(4)*beta + n(7)
G = n(2)*beta**2 + n(5)*beta + n(8)
D = 2.0D0 * G / (-F - sqrt(F**2 - 4.0D0*E*G))
tsat = (n(10) + D - sqrt((n(10) + D)**2 - 4.0D0 * (n(9) + n(10)*D))) / 2.0D0
end function tsat
! -------------------------------------------------------------------------
! Boundary Region 2 / 3: P = f(T) in MPa
! -------------------------------------------------------------------------
double precision function p23(T)
double precision, intent(in) :: T
p23 = 0.34805185628969D+03 - 0.11671859879975D+01 * T + 0.10192970039326D-02 * T**2
end function p23
! -------------------------------------------------------------------------
! Region 1 Properties Solver
! -------------------------------------------------------------------------
subroutine solve_region1(T, P, v, h, s, cp, cv, w)
double precision, intent(in) :: T, P
double precision, intent(out) :: v, h, s, cp, cv, w
double precision :: Tr, Pr, g, gp, gpp, gt, gtt, gpt
integer :: i
Tr = 1386.0D0 / T
Pr = P / 16.53D0
g = 0.0D0
gp = 0.0D0
gpp = 0.0D0
gt = 0.0D0
gtt = 0.0D0
gpt = 0.0D0
do i = 1, 34
g = g + r1_n(i) * (7.1D0 - Pr)**r1_li(i) * (Tr - 1.222D0)**r1_lj(i)
if (r1_li(i) > 0) then
gp = gp - r1_n(i) * r1_li(i) * (7.1D0 - Pr)**(r1_li(i)-1) * (Tr - 1.222D0)**r1_lj(i)
end if
if (r1_li(i) > 1) then
gpp = gpp + r1_n(i) * r1_li(i) * (r1_li(i)-1) * (7.1D0 - Pr)**(r1_li(i)-2) * (Tr - 1.222D0)**r1_lj(i)
end if
if (r1_lj(i) /= 0) then
gt = gt + r1_n(i) * r1_lj(i) * (7.1D0 - Pr)**r1_li(i) * (Tr - 1.222D0)**(r1_lj(i)-1)
end if
if (r1_lj(i) /= 0 .and. r1_lj(i) /= 1) then
gtt = gtt + r1_n(i) * r1_lj(i) * (r1_lj(i)-1) * (7.1D0 - Pr)**r1_li(i) * (Tr - 1.222D0)**(r1_lj(i)-2)
end if
if (r1_li(i) > 0 .and. r1_lj(i) /= 0) then
gpt = gpt - r1_n(i) * r1_li(i) * r1_lj(i) * (7.1D0 - Pr)**(r1_li(i)-1) * (Tr - 1.222D0)**(r1_lj(i)-1)
end if
end do
v = Pr * gp * R * T / P / 1000.0D0
h = Tr * gt * R * T
s = R * (Tr * gt - g)
cp = -R * Tr**2 * gtt
cv = R * (-Tr**2 * gtt + (gp - Tr * gpt)**2 / gpp)
w = sqrt(R * T * 1000.0D0 * gp**2 / ((gp - Tr * gpt)**2 / (Tr**2 * gtt) - gpp))
end subroutine solve_region1
! -------------------------------------------------------------------------
! Region 2cp0 Helper
! -------------------------------------------------------------------------
subroutine region2_cp0(Tr, Pr, go, gop, gopp, got, gott, gopt)
double precision, intent(in) :: Tr, Pr
double precision, intent(out) :: go, gop, gopp, got, gott, gopt
integer :: i
go = log(Pr)
gop = 1.0D0 / Pr
gopp = -1.0D0 / Pr**2
gopt = 0.0D0
got = 0.0D0
gott = 0.0D0
do i = 1, 9
go = go + r2_cp0_no(i) * Tr**r2_cp0_jo(i)
got = got + r2_cp0_no(i) * r2_cp0_jo(i) * Tr**(r2_cp0_jo(i)-1)
gott = gott + r2_cp0_no(i) * r2_cp0_jo(i) * (r2_cp0_jo(i)-1) * Tr**(r2_cp0_jo(i)-2)
end do
end subroutine region2_cp0
! -------------------------------------------------------------------------
! Region 2 Properties Solver
! -------------------------------------------------------------------------
subroutine solve_region2(T, P, v, h, s, cp, cv, w)
double precision, intent(in) :: T, P
double precision, intent(out) :: v, h, s, cp, cv, w
double precision :: Tr, Pr, go, gop, gopp, got, gott, gopt
double precision :: gr, grp, grpp, grt, grtt, grpt
integer :: i
Tr = 540.0D0 / T
Pr = P / 1.0D0
call region2_cp0(Tr, Pr, go, gop, gopp, got, gott, gopt)
gr = 0.0D0
grp = 0.0D0
grpp = 0.0D0
grt = 0.0D0
grtt = 0.0D0
grpt = 0.0D0
do i = 1, 43
gr = gr + r2_n(i) * Pr**r2_li(i) * (Tr - 0.5D0)**r2_lj(i)
if (r2_li(i) > 0) then
grp = grp + r2_n(i) * r2_li(i) * Pr**(r2_li(i)-1) * (Tr - 0.5D0)**r2_lj(i)
end if
if (r2_li(i) > 1) then
grpp = grpp + r2_n(i) * r2_li(i) * (r2_li(i)-1) * Pr**(r2_li(i)-2) * (Tr - 0.5D0)**r2_lj(i)
end if
if (r2_lj(i) > 0) then
grt = grt + r2_n(i) * r2_lj(i) * Pr**r2_li(i) * (Tr - 0.5D0)**(r2_lj(i)-1)
end if
if (r2_lj(i) > 1) then
grtt = grtt + r2_n(i) * r2_lj(i) * (r2_lj(i)-1) * Pr**r2_li(i) * (Tr - 0.5D0)**(r2_lj(i)-2)
end if
if (r2_li(i) > 0 .and. r2_lj(i) > 0) then
grpt = grpt + r2_n(i) * r2_li(i) * r2_lj(i) * Pr**(r2_li(i)-1) * (Tr - 0.5D0)**(r2_lj(i)-1)
end if
end do
v = Pr * (gop + grp) * R * T / P / 1000.0D0
h = Tr * (got + grt) * R * T
s = R * (Tr * (got + grt) - (go + gr))
cp = -R * Tr**2 * (gott + grtt)
cv = R * (-Tr**2 * (gott + grtt) - (1.0D0 + Pr*grp - Tr*Pr*grpt)**2 / (1.0D0 - Pr**2 * grpp))
w = (R*T*1000.0D0*(1.0D0 + 2.0D0*Pr*grp + Pr**2 * grp**2) / &
(1.0D0 - Pr**2 * grpp + (1.0D0 + Pr*grp - Tr*Pr*grpt)**2 / (Tr**2 * (gott + grtt))))**0.5D0
end subroutine solve_region2
! -------------------------------------------------------------------------
! Region 3 Properties Solver (Inputs: density rho in kg/m3 and T in K)
! -------------------------------------------------------------------------
subroutine solve_region3(rho, T, P, v, h, s, cp, cv, w)
double precision, intent(in) :: rho, T
double precision, intent(out) :: P, v, h, s, cp, cv, w
double precision :: d, Tr, g, gd, gdd, gt, gtt, gdt
integer :: i
d = rho / rhoc
Tr = Tc / T
! 1.0658070028513D0 is the coefficient of log(d) in standard Region 3 Helmholtz
g = 1.0658070028513D0 * log(d)
gd = 1.0658070028513D0 / d
gdd = -1.0658070028513D0 / d**2
gt = 0.0D0
gtt = 0.0D0
gdt = 0.0D0
do i = 1, 39
g = g + r3_n(i) * d**r3_li(i) * Tr**r3_lj(i)
if (r3_li(i) > 0) then
gd = gd + r3_n(i) * r3_li(i) * d**(r3_li(i)-1) * Tr**r3_lj(i)
end if
if (r3_li(i) > 1) then
gdd = gdd + r3_n(i) * r3_li(i) * (r3_li(i)-1) * d**(r3_li(i)-2) * Tr**r3_lj(i)
end if
if (r3_lj(i) > 0) then
gt = gt + r3_n(i) * r3_lj(i) * d**r3_li(i) * Tr**(r3_lj(i)-1)
end if
if (r3_lj(i) > 1) then
gtt = gtt + r3_n(i) * r3_lj(i) * (r3_lj(i)-1) * d**r3_li(i) * Tr**(r3_lj(i)-2)
end if
if (r3_li(i) > 0 .and. r3_lj(i) > 0) then
gdt = gdt + r3_n(i) * r3_li(i) * r3_lj(i) * d**(r3_li(i)-1) * Tr**(r3_lj(i)-1)
end if
end do
P = d * gd * R * T * rho / 1000.0D0
v = 1.0D0 / rho
h = R * T * (Tr * gt + d * gd)
s = R * (Tr * gt - g)
cp = R * (-Tr**2 * gtt + (d*gd - d*Tr*gdt)**2 / (2.0D0*d*gd + d**2 * gdd))
cv = -R * Tr**2 * gtt
w = sqrt(R * T * 1000.0D0 * (2.0D0*d*gd + d**2 * gdd - (d*gd - d*Tr*gdt)**2 / (Tr**2 * gtt)))
end subroutine solve_region3
! -------------------------------------------------------------------------
! Region 3 Solver from P and T (finds density first using bisection)
! -------------------------------------------------------------------------
subroutine solve_region3_PT(T, P, v, h, s, cp, cv, w)
double precision, intent(in) :: T, P
double precision, intent(out) :: v, h, s, cp, cv, w
double precision :: rho_min, rho_max, rho_mid, P_calc
double precision :: d_dummy, h_dummy, s_dummy, cp_dummy, cv_dummy, w_dummy
integer :: iter
! Bisection search for density rho
rho_min = 1.0D0
rho_max = 1200.0D0
do iter = 1, 50
rho_mid = 0.5D0 * (rho_min + rho_max)
call solve_region3(rho_mid, T, P_calc, d_dummy, h_dummy, s_dummy, cp_dummy, cv_dummy, w_dummy)
if (P_calc < P) then
if (T < Tc) then
! Saturated liquid region has higher density, so higher density yields higher pressure
rho_min = rho_mid
else
rho_min = rho_mid
end if
else
rho_max = rho_mid
end if
end do
rho_mid = 0.5D0 * (rho_min + rho_max)
call solve_region3(rho_mid, T, P_calc, v, h, s, cp, cv, w)
v = 1.0D0 / rho_mid
end subroutine solve_region3_PT
! -------------------------------------------------------------------------
! Determine Region ID based on T (K) and P (MPa)
! -------------------------------------------------------------------------
integer function get_region(T, P)
double precision, intent(in) :: T, P
double precision :: ps
if (T < 623.15D0) then
ps = psat(T)
if (P >= ps) then
get_region = 1
else
get_region = 2
end if
else if (T <= 863.15D0) then
if (P >= p23(T)) then
get_region = 3
else
get_region = 2
end if
else
get_region = 2
end if
end function get_region
! -------------------------------------------------------------------------
! Master Properties Solver from Pressure and Temperature
! -------------------------------------------------------------------------
subroutine solve_PT(P, T, v, h, s, cp, cv, w, region_id)
double precision, intent(in) :: P, T
double precision, intent(out) :: v, h, s, cp, cv, w
integer, intent(out) :: region_id
region_id = get_region(T, P)
select case (region_id)
case (1)
call solve_region1(T, P, v, h, s, cp, cv, w)
case (2)
call solve_region2(T, P, v, h, s, cp, cv, w)
case (3)
call solve_region3_PT(T, P, v, h, s, cp, cv, w)
case default
call solve_region2(T, P, v, h, s, cp, cv, w)
end select
end subroutine solve_PT
end module iapws_if97_mod
! ============================================================================
! Main Program for iapws_if97_calc
! ============================================================================
program iapws_if97_calc
use iapws_if97_mod
implicit none
integer :: input_type
double precision :: val1, val2
! Conversions
double precision :: T_c, T_k, P_bar, P_mpa
double precision :: v, h, s, u, cp, cv, w, density, quality
character(len=30) :: phase_name
integer :: region_id
! Saturation checks
double precision :: t_sat, h_f, h_g, s_f, s_g, v_f, v_g, cp_f, cp_g, cv_f, cv_g, w_f, w_g
double precision :: t_min, t_max, t_mid, h_mid, s_mid
integer :: iter, reg_dummy
! Read inputs from stdin
read(*,*) input_type
read(*,*) val1
read(*,*) val2
quality = -1.0D0
phase_name = "Single Phase"
select case (input_type)
! -------------------------------------------------------------------------
! 1. Pressure & Temperature (val1 = P [bar], val2 = T [C])
! -------------------------------------------------------------------------
case (1)
P_bar = val1
P_mpa = P_bar / 10.0D0
T_c = val2
T_k = T_c + 273.15D0
if (T_k < 273.15D0 .or. T_k > 1073.15D0) then
write(*,*) 'ERROR: Temperature must be between 0 C and 800 C'
stop
end if
if (P_mpa < Pmin .or. P_mpa > 100.0D0) then
write(*,*) 'ERROR: Pressure must be between 0.006 bar and 1000 bar'
stop
end if
call solve_PT(P_mpa, T_k, v, h, s, cp, cv, w, region_id)
if (region_id == 1) then
phase_name = "Subcooled Liquid"
quality = 0.0D0
else if (region_id == 2) then
phase_name = "Superheated Vapor"
quality = 1.0D0
else if (region_id == 3) then
phase_name = "Supercritical Fluid"
quality = -1.0D0
end if
! -------------------------------------------------------------------------
! 2. Pressure & Quality (val1 = P [bar], val2 = x)
! -------------------------------------------------------------------------
case (2)
P_bar = val1
P_mpa = P_bar / 10.0D0
quality = val2
if (P_mpa < Pmin .or. P_mpa > Pc) then
write(*,*) 'ERROR: Saturated state requires P between 0.006 bar and 220.64 bar'
stop
end if
if (quality < 0.0D0 .or. quality > 1.0D0) then
write(*,*) 'ERROR: Quality x must be between 0 and 1'
stop
end if
T_k = tsat(P_mpa)
T_c = T_k - 273.15D0
region_id = 4
call solve_PT(P_mpa, T_k, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
! For vapor line, we evaluate Region 2 at Saturation Temperature and Pressure
call solve_region2(T_k, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)
! Blend properties
v = (1.0D0 - quality) * v_f + quality * v_g
h = (1.0D0 - quality) * h_f + quality * h_g
s = (1.0D0 - quality) * s_f + quality * s_g
cp = (1.0D0 - quality) * cp_f + quality * cp_g
cv = (1.0D0 - quality) * cv_f + quality * cv_g
w = (1.0D0 - quality) * w_f + quality * w_g
if (quality == 0.0D0) then
phase_name = "Saturated Liquid"
else if (quality == 1.0D0) then
phase_name = "Saturated Vapor"
else
phase_name = "Saturated Mixture"
end if
! -------------------------------------------------------------------------
! 3. Temperature & Quality (val1 = T [C], val2 = x)
! -------------------------------------------------------------------------
case (3)
T_c = val1
T_k = T_c + 273.15D0
quality = val2
if (T_k < 273.15D0 .or. T_k >= Tc) then
write(*,*) 'ERROR: Saturated state requires T between 0 C and 373.94 C'
stop
end if
if (quality < 0.0D0 .or. quality > 1.0D0) then
write(*,*) 'ERROR: Quality x must be between 0 and 1'
stop
end if
P_mpa = psat(T_k)
P_bar = P_mpa * 10.0D0
region_id = 4
call solve_PT(P_mpa, T_k, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
call solve_region2(T_k, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)
! Blend properties
v = (1.0D0 - quality) * v_f + quality * v_g
h = (1.0D0 - quality) * h_f + quality * h_g
s = (1.0D0 - quality) * s_f + quality * s_g
cp = (1.0D0 - quality) * cp_f + quality * cp_g
cv = (1.0D0 - quality) * cv_f + quality * cv_g
w = (1.0D0 - quality) * w_f + quality * w_g
if (quality == 0.0D0) then
phase_name = "Saturated Liquid"
else if (quality == 1.0D0) then
phase_name = "Saturated Vapor"
else
phase_name = "Saturated Mixture"
end if
! -------------------------------------------------------------------------
! 4. Pressure & Enthalpy (val1 = P [bar], val2 = h [kJ/kg])
! -------------------------------------------------------------------------
case (4)
P_bar = val1
P_mpa = P_bar / 10.0D0
h = val2
if (P_mpa < Pmin .or. P_mpa > 100.0D0) then
write(*,*) 'ERROR: Pressure must be between 0.006 bar and 1000 bar'
stop
end if
! 1. Check if we are in two-phase region (Region 4)
if (P_mpa <= Pc) then
T_sat = tsat(P_mpa)
call solve_PT(P_mpa, T_sat, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
call solve_region2(T_sat, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)
if (h >= h_f .and. h <= h_g) then
! Saturation mixture
quality = (h - h_f) / (h_g - h_f)
T_k = T_sat
T_c = T_k - 273.15D0
region_id = 4
v = (1.0D0 - quality) * v_f + quality * v_g
s = (1.0D0 - quality) * s_f + quality * s_g
cp = (1.0D0 - quality) * cp_f + quality * cp_g
cv = (1.0D0 - quality) * cv_f + quality * cv_g
w = (1.0D0 - quality) * w_f + quality * w_g
phase_name = "Saturated Mixture"
end if
end if
! 2. Single-phase region (either P > Pc or h is outside saturation limits)
if (quality < 0.0D0) then
! Check if it is liquid (h < h_f) or vapor (h > h_g)
if (P_mpa <= Pc .and. h < h_f) then
phase_name = "Subcooled Liquid"
t_min = 273.15D0
t_max = T_sat
else if (P_mpa <= Pc .and. h > h_g) then
phase_name = "Superheated Vapor"
t_min = T_sat
t_max = 1073.15D0
else
phase_name = "Supercritical Fluid"
t_min = 273.15D0
t_max = 1073.15D0
end if
! Bisection search for Temperature
do iter = 1, 40
t_mid = 0.5D0 * (t_min + t_max)
call solve_PT(P_mpa, t_mid, v, h_mid, s, cp, cv, w, region_id)
if (h_mid < h) then
t_min = t_mid
else
t_max = t_mid
end if
end do
T_k = 0.5D0 * (t_min + t_max)
T_c = T_k - 273.15D0
call solve_PT(P_mpa, T_k, v, h, s, cp, cv, w, region_id)
if (region_id == 1) quality = 0.0D0
if (region_id == 2) quality = 1.0D0
end if
! -------------------------------------------------------------------------
! 5. Pressure & Entropy (val1 = P [bar], val2 = s [kJ/kg-K])
! -------------------------------------------------------------------------
case (5)
P_bar = val1
P_mpa = P_bar / 10.0D0
s = val2
if (P_mpa < Pmin .or. P_mpa > 100.0D0) then
write(*,*) 'ERROR: Pressure must be between 0.006 bar and 1000 bar'
stop
end if
! 1. Check if we are in two-phase region (Region 4)
if (P_mpa <= Pc) then
T_sat = tsat(P_mpa)
call solve_PT(P_mpa, T_sat, v_f, h_f, s_f, cp_f, cv_f, w_f, reg_dummy)
call solve_region2(T_sat, P_mpa, v_g, h_g, s_g, cp_g, cv_g, w_g)
if (s >= s_f .and. s <= s_g) then
! Saturation mixture
quality = (s - s_f) / (s_g - s_f)
T_k = T_sat
T_c = T_k - 273.15D0
region_id = 4
v = (1.0D0 - quality) * v_f + quality * v_g
h = (1.0D0 - quality) * h_f + quality * h_g
cp = (1.0D0 - quality) * cp_f + quality * cp_g
cv = (1.0D0 - quality) * cv_f + quality * cv_g
w = (1.0D0 - quality) * w_f + quality * w_g
phase_name = "Saturated Mixture"
end if
end if
! 2. Single-phase region (either P > Pc or s is outside saturation limits)
if (quality < 0.0D0) then
if (P_mpa <= Pc .and. s < s_f) then
phase_name = "Subcooled Liquid"
t_min = 273.15D0
t_max = T_sat
else if (P_mpa <= Pc .and. s > s_g) then
phase_name = "Superheated Vapor"
t_min = T_sat
t_max = 1073.15D0
else
phase_name = "Supercritical Fluid"
t_min = 273.15D0
t_max = 1073.15D0
end if
! Bisection search for Temperature
do iter = 1, 40
t_mid = 0.5D0 * (t_min + t_max)
call solve_PT(P_mpa, t_mid, v, h, s_mid, cp, cv, w, region_id)
if (s_mid < s) then
t_min = t_mid
else
t_max = t_mid
end if
end do
T_k = 0.5D0 * (t_min + t_max)
T_c = T_k - 273.15D0
call solve_PT(P_mpa, T_k, v, h, s, cp, cv, w, region_id)
if (region_id == 1) quality = 0.0D0
if (region_id == 2) quality = 1.0D0
end if
case default
write(*,*) 'ERROR: Unknown input combination'
stop
end select
! Density conversion: v in m3/kg => rho = 1/v kg/m3
density = 1.0D0 / v
! Specific Internal Energy u = h - P * v
! P is in bar. 1 bar = 100 kPa.
! P * v is in (100 kPa) * (m3/kg) = 100 kJ/kg.
u = h - P_bar * 100.0D0 * v
! Print output in standard key-value format for PHP wrapper
write(*,'(A,A)') 'STATE=', trim(phase_name)
write(*,'(A,I1)') 'REGION=', region_id
write(*,'(A,F14.4)') 'T=', T_c
write(*,'(A,F14.4)') 'P=', P_bar
write(*,'(A,F14.6)') 'V=', v
write(*,'(A,F14.4)') 'RHO=', density
write(*,'(A,F14.4)') 'H=', h
write(*,'(A,F14.6)') 'S=', s
write(*,'(A,F14.4)') 'U=', u
write(*,'(A,F14.4)') 'CP=', cp
write(*,'(A,F14.4)') 'CV=', cv
write(*,'(A,F14.2)') 'SOUND=', w
write(*,'(A,F14.4)') 'QUALITY=', quality
end program iapws_if97_calc
Solver Description
Calculates complete water and steam properties using the IAPWS-IF97 formulation, including Region IDs, specific volume, density, enthalpy, entropy, internal energy, Cp, Cv, and quality.
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 iapws_if97.f90 -o iapws_if97_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
iapws_if97_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Input Type (1=PT, 2=Px, 3=Tx, 4=Ph, 5=Ps)
1
! Val1 (P [bar] or T [C])
5.0
! Val2 (T [C], x, h [kJ/kg], or s [kJ/kg-K])
150.0
1
! Val1 (P [bar] or T [C])
5.0
! Val2 (T [C], x, h [kJ/kg], or s [kJ/kg-K])
150.0