! ============================================================================ ! 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