🔬
Solver Purpose & Physical Scope
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.
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 499 times
📄 Source File:
iapws_if97.f90
📁 calcul/Thermodynamics /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 iapws_if97.f90 -o iapws_if97_calc
2. Execution with input.txt redirection:
iapws_if97_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 5.0 150.0
Parameter Description:
Input Type (1=PT, 2=Px, 3=Tx, 4=Ph, 5=Ps) Val1 (P [bar] or T [C]) Val2 (T [C], x, h [kJ/kg], or s [kJ/kg-K])