💻 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
Core Numerical Engine in Fortran 90 • 137 total downloads
steam.f90 / steam_calc.f90
! =========================================================================
! Source File: steam.f90
! =========================================================================
subroutine base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
!*****************************************************************************80
!
!! base() calculates quantities associated with the base Helmholtz function.
!
! Discussion:
!
! The equation for the base Helmholtz function AB(T,RHO) is:
!
! AB(T,RHO) = R * T * (
! - ln ( 1 - y )
! - ( beta - 1 ) / ( 1 - y )
! + ( alpha + beta + 1 ) / ( 2 * ( 1 - y )^2 )
! + 4 * y * ( ( Bbar / b ) - gamma )
! - 0.5 * ( alpha - beta + 3 )
! + ln ( RHO * R * T / P0 ) )
! (Equation 2)
! where
!
! y = b * rho / 4,
! alpha = 11,
! beta = 133/3,
! gamma = 7/2,
! P0 = 0.101325 MegaPascals = 1 atm
!
! and
!
! b(T) = b1 * ln(T/T0) + sum(j=0,1,3,5) b(j)*(T0/T)^j (Equation 3)
!
! Bbar(T) = sum(j=0,1,2,4) B(j)*(T0/T)^j (Equation 4).
!
! where
!
! T0=647.073 K and the coefficients b(j) and B(j) are
!
! j b(j) B(j)
! -- ----------- ----------
! 0 0.7478629 1.1278334
! 1 -0.3540782 -0.5944001
! 2 0 -5.010996
! 3 0.007159876 0
! 4 0 0.63684256
! 5 -0.003528426 0
!
! For the derived quantities, the following relations are used:
!
! Pressure: PB = RHO^2 * dAB/dRHO
! Density derivative: DPDRB = 2*PB/RHO + RHO^2 * d2AB/dRHO2
! Temperature derivative: DPDTB = RHO^2 * d2AB/(dRHO dT)
! Specific entropy: SB = ( UB - AB ) / T
! Specific internal energy: UB = AB + T * SB
! Specific enthalpy: HB = UB + PB / RHO
! Specific heat capacity
! at constant volume: CVB = - T * d2AB/dT2
! Specific Gibbs function: GB = AB + PB / RHO
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 03 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Input:
!
! real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output:
!
! real ( kind = rk8 ) AB, the base value of the Helmholtz function,
! in KJ/kg.
!
! real ( kind = rk8 ) CVB, the base value of the isochoric (constant
! volume) heat capacity, in KJ/(kg degrees Kelvin).
!
! real ( kind = rk8 ) DPDRB, the base value of the partial
! derivative dP(T,RHO)/dRHO, with T held fixed, in (MegaPascals CM3)/G.
!
! real ( kind = rk8 ) DPDTB, the base value of the partial
! derivative dP(T,RHO)/dT, with RHO held fixed, in
! MegaPascals/degrees Kelvin.
!
! real ( kind = rk8 ) GB, the base value of the Gibbs free energy,
! in KJ/kg.
!
! real ( kind = rk8 ) HB, the base value of enthalpy, in KJ/kg.
!
! real ( kind = rk8 ) PB, the base pressure, in MegaPascals.
!
! real ( kind = rk8 ) SB, the base value of entropy,
! in KJ/(kg degrees Kelvin).
!
! real ( kind = rk8 ) UB, the base value of internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ), parameter :: alpha = 11.0D+00
real ( kind = rk8 ) b1
real ( kind = rk8 ) b1t
real ( kind = rk8 ) b1tt
real ( kind = rk8 ) b2
real ( kind = rk8 ) b2t
real ( kind = rk8 ) b2tt
real ( kind = rk8 ), parameter :: beta = 44.333333333333D+00
real ( kind = rk8 ) cvb
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dz
real ( kind = rk8 ) dz0
real ( kind = rk8 ), parameter :: gamma = 3.5D+00
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gb
real ( kind = rk8 ) hb
real ( kind = rk8 ), parameter :: p_zero = 0.101325D+00
real ( kind = rk8 ) pb
real ( kind = rk8 ) rho
real ( kind = rk8 ) sb
real ( kind = rk8 ) t
real ( kind = rk8 ) ub
real ( kind = rk8 ) x
real ( kind = rk8 ) y
real ( kind = rk8 ) z
real ( kind = rk8 ) z0
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BASE - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BASE - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
!
! Compute auxilliary quantities for Equation 2.
!
call bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
y = 0.25D+00 * b1 * rho
x = 1.0D+00 - y
!
! Evaluate Equation 2.
!
ab = - log ( 1.0D+00 - y ) &
- ( beta - 1.0D+00 ) / ( 1.0D+00 - y ) &
+ ( alpha + beta + 1.0D+00 ) / ( 2.0D+00 * ( 1.0D+00 - y )**2 ) &
+ 4.0D+00 * y * ( ( b2 / b1 ) - gamma ) &
- 0.5D+00 * ( alpha - beta + 3.0D+00 ) &
+ log ( rho * gascon() * t / p_zero )
!
! Determine quantities defined in terms of AB.
!
pb = ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**3 &
+ 4.0D+00 * y * ( b2 / b1 - gamma )
z0 = ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**3
z = z0 + 4.0D+00 * y * ( b2 / b1 - gamma )
dz0 = ( alpha + 2.0D+00 * beta * y ) / ( 1.0D+00 - y )**3 &
+ 3.0D+00 * ( 1.0D+00 + alpha * y + beta * y**2 ) / ( 1.0D+00 - y )**4
dz = dz0 + 4.0D+00 * ( b2 / b1 - gamma )
gb = ab + pb
ub = - t * b1t * ( pb - 1.0D+00 - rho * b2 ) / b1 - rho * t * b2t
hb = pb + ub
!
! An incorrect version of this equation began:
!
! cvb = 2.0D+00 * ub + ( pb - 1.0D+00 ) &
!
! and caused me no end of trouble. My fault, JVB, 03 February 2002
!
cvb = 2.0D+00 * ub + ( z0 - 1.0D+00 ) &
* ( ( t * b1t / b1 )**2 - t**2 * b1tt / b1 ) &
- rho * t**2 * ( b2tt - gamma * b1tt ) - ( t * b1t / b1 )**2 * y * dz0
dpdtb = pb / t + rho * ( 0.25D+00 * ( dz0 + 4.0D+00 * ( b2 / b1 - gamma ) ) &
* b1t + b2t - b2 / b1 * b1t )
sb = ub - ab
dpdrb = pb + y * ( dz0 + 4.0D+00 * ( b2 / b1 - gamma ) )
!
! Assign dimensions.
!
ab = gascon() * t * ab
cvb = gascon() * cvb
dpdrb = gascon() * t * dpdrb
dpdtb = gascon() * t * rho * dpdtb
gb = gascon() * t * gb
hb = gascon() * t * hb
pb = gascon() * t * rho * pb
sb = gascon() * sb
ub = gascon() * t * ub
return
end
subroutine bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
!*****************************************************************************80
!
!! bb() calculates the B's of equations 3 and 4.
!
! Discussion:
!
! Here
!
! b(T) = b1 * ln(T/T0) + sum(j=0,1,3,5) b(j)*(T0/T)**j (Equation 3)
!
! Bbar(T) = sum(j=0,1,2,4) B(j)*(T0/T)**j (Equation 4).
!
! where
!
! T0 = 647.073 K
!
! and the coefficients b(j) and B(j) are
!
! j b(j) B(j)
! -- ----------- ----------
! 0 0.7478629 1.1278334
! 1 -0.3540782 -0.5944001
! 2 0 -5.010996
! 3 0.007159876 0
! 4 0 0.63684256
! 5 -0.003528426 0
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) B1, the coefficient b from equation 3,
! in CM3/G.
!
! Output, real ( kind = rk8 ) B2, the coefficient Bbar from equation 4,
! in CM3/G.
!
! Output, real ( kind = rk8 ) B1T, the derivative dB1/dT,
! in (CM3)/(G Degrees Kelvin).
!
! Output, real ( kind = rk8 ) B2T, the derivative dB2/dT,
! in (CM3)/(G Degrees Kelvin).
!
! Output, real ( kind = rk8 ) B1TT, the second derivative of B1 with
! respect to T, in (CM3)/(G (Degrees Kelvin)^2 ).
!
! Output, real ( kind = rk8 ) B2TT, the second derivative of B2 with
! respect to T, in (CM3)/(G (Degrees Kelvin)^2 ).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) b1
real ( kind = rk8 ) b1t
real ( kind = rk8 ) b1tt
real ( kind = rk8 ) b2
real ( kind = rk8 ) b2t
real ( kind = rk8 ) b2tt
real ( kind = rk8 ), parameter, dimension ( 10 ) :: bp = (/ &
0.7478629D+00, -0.3540782D+00, 0.0D+00, 0.0D+00, &
0.007159876D+00, 0.0D+00, -0.003528426D+00, 0.0D+00, &
0.0D+00, 0.0D+00 /)
real ( kind = rk8 ), parameter, dimension ( 10 ) :: bq = (/ &
1.1278334D+00, 0.0D+00, -0.5944001D+00, -5.010996D+00, &
0.0D+00, 0.63684256D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00 /)
integer i
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.073D+00
real ( kind = rk8 ) v(10)
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BB - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Set V(I) = ( T_REF / T )**(I-1).
!
v(1) = 1.0D+00
do i = 2, 10
v(i) = v(i-1) * t_ref / t
end do
!
! Set B1, B1T, B1TT.
!
b1 = bp(1) + bp(2) * log ( 1.0D+00 / v(2) )
b1t = bp(2) * v(2) / t_ref
b1tt = 0.0D+00
do i = 3, 10
b1 = b1 + bp(i) * v(i-1)
b1t = b1t - dble ( i - 2 ) * bp(i) * v(i-1) / t
b1tt = b1tt + bp(i) * dble ( i - 2 )**2 * v(i-1) / t**2
end do
b1tt = b1tt - ( b1t / t )
!
! Set B2, B2T, B2TT.
!
b2 = bq(1)
b2t = 0.0D+00
b2tt = 0.0D+00
do i = 3, 10
b2 = b2 + bq(i) * v(i-1)
b2t = b2t - dble ( i - 2 ) * bq(i) * v(i-1) / t
b2tt = b2tt + bq(i) * dble ( i - 2 )**2 * v(i-1) / t**2
end do
b2tt = b2tt - ( b2t / t )
return
end
subroutine corr ( t, p, p_consistent, rhol, rhov, delg )
!*****************************************************************************80
!
!! corr() evaluates an adjustment to the Gibbs function.
!
! Discussion:
!
! CORR is given T and P at or near the vapor pressure and evaluates
! the corresponding liquid and vapor densities, and the residual
! function DELG = (GL-GV)/(R*T) where GL and GV are the Gibbs functions
! for the liquid and vapor phases, respectively.
!
! These quantities are used to calculate a correction to the vapor
! pressure or the vapor temperature.
!
! The states corresponding to the coexisting phases of liquid
! and vapor for the temperature range from the triple point
! to within 0.5 C of the critical point 0.01 <= t <= tk-0.5 C
! have been determined in exact accord with the Gibbs condition
! of phase equilibrium: DELG = G(g)-G(l) = 0, P, t constant,
! where G(g) and G(l) are the values of the Gibbs function
! for saturated gas and liquid respectively.
!
! For the region (tk-t)<=0.5 C, an exact solution for the
! Helmholtz function yields values of density for the saturated
! liquid that are shifted to lower values. Also, the isotherms
! in the pressure-density plane and the Gibbs function-density
! plane are nearly flat, so that it is difficult to obtain
! solutions. As an alternative to exact solution, the power
! law equation is used to define states:
!
! rho(gas) = 0.322 - 0.657 * (1 - T/647.126)**0.325 (g/cm3).
! rho(liq) = 0.322 + 0.657 * (1 - T/647.126)**0.325 (g/cm3).
!
! In a poor instance of programming, the input pressure was
! originally overwritten on output by a value consistent with
! the computed densities. This causes no end of misunderstandings,
! since other routines expect the value of pressure to be input
! only. The code is now revised so that there is an input P
! and an output P. In a huff, JVB 05 February 2002.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 05 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the vapor temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) P_CONSISTENT, the vapor pressure, in MegaPascals,
! consistent with RHOL and RHOV. This is equal to the input value of
! P unless 646.3 <= T.
!
! Input/output, real ( kind = rk8 ) RHOL, the liquid density, in G/CM3.
! On input, if RHOL is positive, it is used as an initial
! estimate for the iteration.
!
! Input/output, real ( kind = rk8 ) RHOV, the vapor density, in G/CM3.
! On input, if RHOV is positive, it is used as an initial
! estimate for the iteration.
!
! Output, real ( kind = rk8 ) DELG, the residual function (GL-GV)/(R*T),
! where GL is the liquid Gibbs function, GV the vapor Gibbs function,
! dimensionless. If 646.3 < T, DELG is 0.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ) ab
real ( kind = rk8 ) ar
real ( kind = rk8 ) cd
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cv
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvr
logical, parameter :: debug = .false.
real ( kind = rk8 ) delg
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) g
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gb
real ( kind = rk8 ) gl
real ( kind = rk8 ) gr
real ( kind = rk8 ) gv
real ( kind = rk8 ) h
real ( kind = rk8 ) hb
real ( kind = rk8 ) hr
real ( kind = rk8 ) p
real ( kind = rk8 ) p_consistent
real ( kind = rk8 ), parameter :: p_crit = 22.055D+00
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) rhol
real ( kind = rk8 ), parameter :: rho_min = 1.0D-08
real ( kind = rk8 ) rhov
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) s
real ( kind = rk8 ) sb
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) tau
real ( kind = rk8 ), parameter :: t_crit = 647.1260000001D+00
real ( kind = rk8 ) u
real ( kind = rk8 ) ub
real ( kind = rk8 ) ur
p_consistent = p
!
! Initialize output quantities.
!
delg = 0.0D+00
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'CORR - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'CORR - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was P = ', p
stop
end if
if ( t <= 646.3D+00 ) then
if ( rhol <= 0.0D+00 ) then
rho_start = 1.11D+00 - 0.0004D+00 * t
else
rho_start = rhol
end if
call dense ( p_consistent, t, rho_start, rho, dpdr )
call therm ( t, rho, a, cjth, cjtt, cd, cv, dpdr, dpdt, g, h, &
p_consistent, s, u )
rhol = rho
gl = g
if ( rhov <= 0.0D+00 ) then
rho_start = p_consistent / ( gascon() * t )
else
rho_start = rhov
end if
call dense ( p_consistent, t, rho_start, rho, dpdr )
rho = max ( rho, rho_min )
call therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, &
p_consistent, s, u )
rhov = rho
gv = g
delg = ( gl - gv ) / ( gascon() * t )
p_consistent = p
if ( debug ) then
write ( *, '(a,g14.6)' ) ' CORR - RHOL = ', rhol
write ( *, '(a,g14.6)' ) ' RHOV = ', rhov
end if
else if ( t <= t_crit ) then
if ( debug ) then
write ( *, '(a)' ) 'CORR - Twilight zone'
end if
delg = 0.0D+00
tau = 0.657128D+00 * ( 1.0D+00 - t / t_crit )**0.325D+00
rhol = 0.322D+00 + tau
rhov = 0.322D+00 - tau
rho = rhov
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
p_consistent = pb + pr
else
if ( debug ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'CORR - Weirdo zone'
end if
rhol = 0.322D+00
rhov = 0.322D+00
p_consistent = p_crit
delg = 0.0D+00
end if
return
end
subroutine cp_values ( n, tc, p, cp )
!*****************************************************************************80
!
!! cp_values() returns some values of the specific heat at constant pressure.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 229-237.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) CP, the specific heat at constant pressure,
! in KJ/(kg K).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 24
real ( kind = rk8 ) cp
real ( kind = rk8 ), save, dimension ( n_data ) :: cpvec = (/ &
4.228D+00, 2.042D+00, 1.975D+00, 2.013D+00, 2.040D+00, &
2.070D+00, 2.135D+00, 2.203D+00, 2.378D+00, 2.541D+00, &
2.792D+00, 2.931D+00, 4.226D+00, 4.223D+00, 4.202D+00, &
4.177D+00, 4.130D+00, 4.089D+00, 4.053D+00, 4.021D+00, &
3.909D+00, 3.844D+00, 3.786D+00, 2.89D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 5.0D+00, 10.0D+00, 50.0D+00, &
100.0D+00, 200.0D+00, 300.0D+00, 400.0D+00, 500.0D+00, &
1000.0D+00, 1500.0D+00, 2000.0D+00, 5000.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 100.0D+00, 200.0D+00, 300.0D+00, 350.0D+00, &
400.0D+00, 500.0D+00, 600.0D+00, 850.0D+00, 1100.0D+00, &
1600.0D+00, 2000.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
cp = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
cp = cpvec(n)
end if
return
end
subroutine dense ( p, t, rho_start, rho, dpdr )
!*****************************************************************************80
!
!! dense() computes the density for a given pressure and temperature.
!
! Discussion:
!
! The use of the variable RHO_START for two opposing purposes is
! poor practice and will be corrected one of these days. Meanwhile,
! the algorithm's behavior, particularly in the two-phase region,
! is very suspect.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 19 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO_START, an initial guess for the density,
! in G/CM3. The value of RHO_START also signals whether a vapor or liquid
! calculation is to be done. If DPDR is computed negative, then for
! 0.2967 <= RHO_START, liquid is assumed, otherwise gas.
!
! Output, real ( kind = rk8 ) RHO, the density for the given
! pressure and temperature, in G/CM3.
!
! Output, real ( kind = rk8 ) DPDR, the partial derivative
! dP(T,RHO)/dRHO, with T held fixed, in (MegaPascals CM3)/G.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ) ar
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dp
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) dpdx
real ( kind = rk8 ) errtol
real ( kind = rk8 ) gb
real ( kind = rk8 ) gr
real ( kind = rk8 ) hb
real ( kind = rk8 ) hr
integer it
integer, parameter :: it_max = 50
real ( kind = rk8 ) p
real ( kind = rk8 ) pb
real ( kind = rk8 ) pp
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ), parameter :: rho_max = 1.9D+00
real ( kind = rk8 ), parameter :: rho_min = 1.0D-08
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) sb
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) ub
real ( kind = rk8 ) ur
real ( kind = rk8 ) x
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DENSE - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DENSE - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was P = ', p
stop
end if
rho = rho_start
rho = max ( rho, rho_min )
rho = min ( rho, rho_max )
do it = 1, it_max
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
pp = pb + pr
dpdr = dpdrb + dpdrr
!
! Check for negative DP/DRho, which characterizes the two-phase region.
!
if ( dpdr <= 0.0D+00 ) then
if ( 0.2967D+00 <= rho_start ) then
rho = rho * 1.02D+00
else
rho = rho * 0.98D+00
end if
if ( it <= 10 ) then
cycle
end if
end if
dpdx = 1.1D+00 * dpdr
dpdx = max ( dpdx, 0.01D+00 )
dp = abs ( 1.0D+00 - pp / p )
if ( dp <= errtol .or. &
( 0.3D+00 < rho .and. dp <= errtol ) .or. &
( 0.7D+00 < rho .and. dp <= 10.0D+00 * errtol ) ) then
return
end if
x = ( p - pp ) / dpdx
if ( 0.1D+00 < abs ( x ) ) then
x = x * 0.1D+00 / abs ( x )
end if
rho = rho + x
rho = max ( rho, rho_min )
rho = min ( rho, rho_max )
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DENSE - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' Number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Last iterate was ', rho
return
end
subroutine dielectric ( t, rho, eps )
!*****************************************************************************80
!
!! dielectric() returns the static dielectric constant.
!
! Discussion:
!
! According to the IAPS, the equation used is valid in the range
!
! 273.15 degrees Kelvin <= T <= 823.15 degrees K
! 0 MegaPascals <= P <= 500 MegaPascals.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 266.
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) EPS, the dielectric constant, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: npol_c = 4
real ( kind = rk8 ), parameter, dimension ( 10 ) :: a = (/ &
7.62571D+00, 244.003D+00, -140.569D+00, 27.7841D+00, -96.2805D+00, &
41.7909D+00, -10.2099D+00, -45.2059D+00, 84.6395D+00, -35.8644D+00 /)
real ( kind = rk8 ) c(0:npol_c)
real ( kind = rk8 ) eps
real ( kind = rk8 ) rho
real ( kind = rk8 ) t
real ( kind = rk8 ) t_copy
real ( kind = rk8 ), parameter :: t_max = 823.15D+00
real ( kind = rk8 ), parameter :: t_min = 273.15D+00
real ( kind = rk8 ), parameter :: t_ref = 298.15D+00
real ( kind = rk8 ) t_star
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DIELECTRIC - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DIELECTRIC - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
t_copy = t
t_copy = min ( t_copy, t_max )
t_copy = max ( t_copy, t_min )
t_star = t_copy / t_ref
c(0) = 1.0D+00
c(1) = a(1) / t_star
c(2) = ( a(2) / t_star ) + a(3) + a(4) * t_star
c(3) = ( a(5) / t_star ) + a(6) * t_star + a(7) * t_star**2
c(4) = ( a(8) / t_star**2 ) + ( a(9) / t_star ) + a(10)
call r8poly_val_horner ( npol_c, c, rho, eps )
return
end
subroutine dielectric_values ( n, tc, p, eps )
!*****************************************************************************80
!
!! dielectric_values() returns some values of the static dielectric constant.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 03 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 266.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) EPS, the dielectric constant, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 15
real ( kind = rk8 ) eps
real ( kind = rk8 ), save, dimension ( n_data ) :: epsvec = (/ &
88.29D+00, 90.07D+00, 92.02D+00, 95.14D+00, 100.77D+00, &
78.85D+00, 70.27D+00, 62.60D+00, 55.78D+00, 44.31D+00, &
35.11D+00, 20.40D+00, 1.17D+00, 1.11D+00, 1.08D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
100.0D+00, 500.0D+00, 1000.0D+00, 2000.0D+00, 5000.0D+00, &
100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00, &
100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00, 100.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, 150.0D+00, &
200.0D+00, 300.0D+00, 400.0D+00, 500.0D+00, 600.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
eps = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
eps = epsvec(n)
end if
return
end
function gascon ( )
!*****************************************************************************80
!
!! gascon() returns the value of the specific gas constant.
!
! Discussion:
!
! The specific gas constant R is related to the universal gas
! constant R-bar = 8.31441 J/(mol degrees Kelvin) by the molar mass
! M = 18.0152 g/mol:
!
! R = R-bar / M.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Output, real ( kind = rk8 ) GASCON, the value of the specific gas
! constant, in J/(g degrees Kelvin).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) gascon
gascon = 0.461522D+00
return
end
subroutine ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
!*****************************************************************************80
!
!! ideal() computes ideal gas thermodynamic properties of water.
!
! Discussion:
!
! Values for thermodynamic properties of water in the ideal
! gas state were reported by Woolley. The formula for the ideal gas
! term of the Helmholtz function approximates a term by term summation of
! contributions from each of the rotation and vibration states.
! The formula, equation #6 in the reference, is:
!
! A(ideal)(T) = -R * T * ( 1 + ( C(1)/Tr + C(2) ) * ln(Tr)
! + Sum ( 3 <= I <= 18) C(I) * Tr**(I-6)
!
! where Tr=T/100 K. The C(i) are tabulated coefficients. Equation
! 6 can be used for temperatures below 3000 K, and is accurate to
! within the tolerance of the gas constant for 50<=T<=2000 K.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) AI, the Helmholtz function, in KJ/kg.
!
! Output, real ( kind = rk8 ) CPI, the heat capacity at constant pressure,
! in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) CVI, the heat capacity at constant volume,
! in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) GI, the Gibbs free energy, in KJ/kg.
!
! Output, real ( kind = rk8 ) HI, the enthalpy, in KJ/kg.
!
! Output, real ( kind = rk8 ) SI, the entropy, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) UI, the internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ai
real ( kind = rk8 ), parameter, dimension ( 18 ) :: c = (/ &
19.730271018D+00, 20.9662681977D+00, -0.483429455355D+00, &
6.05743189245D+00, 22.56023885D+00, -9.87532442D+00, &
-4.3135538513D+00, 0.458155781D+00, -0.047754901883D+00, &
0.0041238460633D+00, -0.00027929052852D+00, 0.14481695261D-04, &
-0.56473658748D-06, 0.16200446D-07, -0.3303822796D-09, &
0.451916067368D-11, -0.370734122708D-13, 0.137546068238D-15 /)
real ( kind = rk8 ) cpi
real ( kind = rk8 ) cvi
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gi
real ( kind = rk8 ) hi
integer i
real ( kind = rk8 ) si
real ( kind = rk8 ) t
real ( kind = rk8 ) tt
real ( kind = rk8 ) ui
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDEAL - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
tt = t / 100.0D+00
gi = - ( c(1) / tt + c(2) ) * log ( tt )
do i = 3, 18
gi = gi - c(i) * tt**(i-6)
end do
hi = c(2) + c(1) * ( 1.0D+00 - log ( tt ) ) / tt
do i = 3, 18
hi = hi + dble ( i - 6 ) * c(i) * tt**(i-6)
end do
cpi = c(2) - c(1) / tt
do i = 3, 18
cpi = cpi + dble ( ( i - 6 ) * ( i - 5 ) ) * c(i) * tt**(i-6)
end do
ai = gi - 1.0D+00
ui = hi - 1.0D+00
cvi = cpi - 1.0D+00
si = hi - gi
!
! Assign dimensions.
!
ai = gascon() * t * ai
cpi = gascon() * cpi
cvi = gascon() * cvi
gi = gascon() * t * gi
hi = gascon() * t * hi
si = gascon() * si
ui = gascon() * t * ui
return
end
subroutine prandtl ( t, p, pr )
!*****************************************************************************80
!
!! prandtl() computes the Prandtl number.
!
! Discussion:
!
! This routine was NOT working properly for large pressures,
! because the routine CORR was changing the input value of P.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 17 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) PR, the Prandtl number, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cv
logical, parameter :: debug = .false.
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) eta
real ( kind = rk8 ) g
real ( kind = rk8 ) h
real ( kind = rk8 ) lambda
real ( kind = rk8 ) p
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) rhov
real ( kind = rk8 ) s
real ( kind = rk8 ) t
real ( kind = rk8 ) t2
real ( kind = rk8 ) u
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PRANDTL - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PRANDTL - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was P = ', p
stop
end if
!
! Compute the density.
!
if ( debug ) then
write ( *, * ) 'PRANDTL - Call TSAT, with P = ', p
end if
rhol = 0.0D+00
rhov = 0.0D+00
call tsat ( p, t2, rhol, rhov )
if ( debug ) then
write ( *, * ) 'PRANDTL - T2 = ', t2
end if
if ( t < t2 ) then
rho_start = 1.9D+00
else
rho_start = 0.01D+00
end if
call dense ( p, t, rho_start, rho, dpdr )
if ( debug ) then
write ( *, * ) 'PRANDTL - RHO = ', rho
end if
!
! Now from T and RHO, compute CP, ETA and LAMBDA.
!
call therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
call viscosity ( t, rho, eta )
call thercon ( t, rho, lambda )
if ( debug ) then
write ( *, '(7f10.4)' ) t, p, rho, eta, cp, lambda, pr
end if
pr = eta * cp / lambda
return
end
subroutine prandtl_values ( n, tc, p, pr )
!*****************************************************************************80
!
!! prandtl_values() returns some values of the Prandtl number for testing.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 265.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) PR, the Prandtl number, dimensionless.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 35
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ) pr
real ( kind = rk8 ), save, dimension ( n_data ) :: prvec = (/ &
13.50D+00, 13.48D+00, 13.46D+00, 13.39D+00, 13.27D+00, &
13.15D+00, 13.04D+00, 12.93D+00, 12.83D+00, 12.73D+00, &
12.63D+00, 12.53D+00, 12.43D+00, 12.34D+00, 12.25D+00, &
12.08D+00, 11.92D+00, 11.77D+00, 11.62D+00, 11.48D+00, &
11.36D+00, 11.23D+00, 11.12D+00, 10.91D+00, 10.72D+00, &
10.55D+00, 6.137D+00, 3.555D+00, 2.378D+00, 1.000D+00, &
0.974D+00, 0.960D+00, 0.924D+00, 0.899D+00, 0.882D+00 /)
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 5.0D+00, 10.0D+00, 25.0D+00, 50.0D+00, &
75.0D+00, 100.0D+00, 125.0D+00, 150.0D+00, 175.0D+00, &
200.0D+00, 225.0D+00, 250.0D+00, 275.0D+00, 300.0D+00, &
350.0D+00, 400.0D+00, 450.0D+00, 500.0D+00, 550.0D+00, &
600.0D+00, 650.0D+00, 700.0D+00, 800.0D+00, 900.0D+00, &
1000.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, &
150.0D+00, 200.0D+00, 400.0D+00, 600.0D+00, 800.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
pr = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
pr = prvec(n)
end if
return
end
subroutine psat ( t, p, rhol, rhov )
!*****************************************************************************80
!
!! psat() calculates the vapor pressure, and the liquid and vapor densities.
!
! Discussion:
!
! These quantities correspond to the input temperature T, corrected
! so that the Gibbs functions for liquid and vapor phase are
! equal to within a tolerance.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the vapor temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) RHOL, the liquid density, in G/CM3.
!
! Output, real ( kind = rk8 ) RHOV, the vapor density, in G/CM3.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) bot
real ( kind = rk8 ) delg
real ( kind = rk8 ) dp
real ( kind = rk8 ) errtol
real ( kind = rk8 ) gascon
integer it
integer, parameter :: it_max = 100
real ( kind = rk8 ) p
real ( kind = rk8 ) p_consistent
real ( kind = rk8 ) p_old
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rhov
real ( kind = rk8 ) t
!
! Ensure that output quantities are initialized,, obliterating any
! input values.
!
p = 0.0D+00
rhol = 0.0D+00
rhov = 0.0D+00
!
! Set the error tolerance.
!
errtol = 100.0D+00 * sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Get an estimate for the saturation pressure.
!
call psat_est ( t, p )
dp = 0.0D+00
do it = 1, it_max
call corr ( t, p, p_consistent, rhol, rhov, delg )
bot = ( rhol - rhov ) / ( rhol * rhov )
if ( abs ( bot ) < errtol ) then
write ( *, * ) 'PSAT - Warning, what is this?'
bot = sign ( errtol, bot )
end if
dp = delg * gascon() * t / bot
p_old = p
p = p + dp
if ( abs ( dp ) <= errtol * ( abs ( p ) + 1.0D+00 ) ) then
return
end if
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Warning!'
write ( *, '(a)' ) ' The iterates have become nonpositive.'
write ( *, '(a,i6)' ) ' Iteration number = ', it
write ( *, '(a,g14.6)' ) ' Last iterate was ', p
write ( *, '(a,g14.6)' ) ' Previous iterate was ', p_old
write ( *, '(a,g14.6)' ) ' Last correction was ', dp
write ( *, '(a)' ) ' Trying to recover...'
p = 0.5D+00 * p_old
end if
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Fatal error!'
write ( *, '(a)' ) ' The iterates have become nonpositive.'
write ( *, '(a,i6)' ) ' Iteration number = ', it
write ( *, '(a,g14.6)' ) ' Last iterate was ', p
write ( *, '(a,g14.6)' ) ' Previous iterate was ', p_old
write ( *, '(a,g14.6)' ) ' Last correction was ', dp
stop
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' The number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Convergence tolerance was ', errtol
write ( *, '(a,g14.6)' ) ' Last iterate was ', p
write ( *, '(a,g14.6)' ) ' Last correction was ', dp
return
end
subroutine psat_est ( t, p )
!*****************************************************************************80
!
!! psat_est() makes a rough estimate of the vapor pressure.
!
! Discussion:
!
! The calculation agrees with tabulated data to within
! 0.02% for temperature to within a degree or so of the critical
! temperature. The approximate vapor pressure can be refined
! by imposing the condition that the Gibbs functions of the vapor
! and liquid phases be equal.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 21 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter, dimension ( 8 ) :: a = (/ &
-7.8889166D+00, 2.5514255D+00, -6.716169D+00, 33.239495D+00, &
-105.38479D+00, 174.35319D+00, -148.39348D+00, 48.631602D+00 /)
real ( kind = rk8 ) b
integer i
real ( kind = rk8 ) p
real ( kind = rk8 ) q
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.25D+00
real ( kind = rk8 ) v
real ( kind = rk8 ) w
real ( kind = rk8 ) z
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PSAT_EST - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
if ( t <= 314.0D+00 ) then
p = 0.1D+00 * exp ( 6.3573118D+00 - 8858.843D+00 / t &
+ 607.56335D+00 * t**( -0.6D+00 ) )
else
v = t / t_ref
w = abs ( 1.0D+00 - v )
b = 0.0D+00
do i = 1, 8
z = i
b = b + a(i) * w**( ( z + 1.0D+00 ) / 2.0D+00 )
end do
q = b / v
p = 22.093D+00 * exp ( q )
end if
return
end
subroutine psat_values ( n, tc, p )
!*****************************************************************************80
!
!! psat_values() returns some values of the saturation pressure.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 9-15.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the saturation pressure, in bar.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 12
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
0.0061173D+00, 0.0065716D+00, 0.0087260D+00, 0.12344D+00, 1.0132D+00, &
2.3201D+00, 4.7572D+00, 15.537D+00, 39.737D+00, 85.838D+00, &
165.21D+00, 220.55D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.01D+00, 1.0D+00, 5.0D+00, 50.0D+00, 100.0D+00, &
125.0D+00, 150.0D+00, 200.0D+00, 250.0D+00, 300.0D+00, &
350.0D+00, 373.976D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
end if
return
end
subroutine r8poly_val_horner ( n, c, x, cx )
!*****************************************************************************80
!
!! r8poly_val_horner() evaluates a polynomial using Horner's method.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 08 August 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer N, the dimension of C.
!
! Input, real ( kind = rk8 ) C(0:N), the polynomial coefficients.
! C(I) is the coefficient of X^I.
!
! Input, real ( kind = rk8 ) X, the point at which the polynomial
! is to be evaluated.
!
! Output, real ( kind = rk8 ) CX, the value of the polynomial at X.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer n
real ( kind = rk8 ) c(0:n)
real ( kind = rk8 ) cx
integer i
real ( kind = rk8 ) x
cx = c(n)
do i = n - 1, 0, -1
cx = cx * x + c(i)
end do
return
end
subroutine resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
!*****************************************************************************80
!
!! resid() calculates residual contributions to thermodynamic quantities.
!
! Discussion:
!
! The residual function consists of 40 terms. The first 36 are
! used in a global least squares fit to experimental data.
!
! Three terms were added that contribute only in the immediate
! neighborhood of the critical point
! (tk-5) <= T <= (tk+5) C
! 0.20 <= rho <= 0.44 g/cm3,
!
! A single term was added for the region of high pressure and
! low temperature: T < 75 C, 300 MPa < P.
!
! Except in these limited regions, the residual function is
! given by the first 36 terms. The equation is
!
! A(residual)(rho,T)=
! sum(i=1 to 36) (g(i)/k(i)) * (T0/T)**(l(i)) (1-exp(-rho))^(k(i))
! + sum(i=37 to 40) g(i)*delta(i)^(k(i))
! * exp(-alpha(i)*delta(i)^(k(i)) - beta(i)*tau(i)^2)
! (Equation 5)
!
! where
!
! g(i) are coefficients determined by fits to data,
! delta(i) are reduced densities (delta(i)=((rho-rho(i))/rho(i))
! tau(i) are reduced temperatures (tau(i)=((T-tau(i))/tau(i))
! rho(i) are specified densities.
! tau(i) are specified temperatures.
! The k(i) and l(i) are specified integers.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 22 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) AR, the residual contribution to the
! Helmholtz function, in KJ/kg.
!
! Output, real ( kind = rk8 ) CVR, the residual contribution to the
! isochoric (constant volume) heat capacity, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) DPDRR, the residual contribution to
! the partial derivative dP(T,RHO)/dRHO, with T held fixed, in
! (MegaPascals CM3)/G.
!
! Output, real ( kind = rk8 ) DPDTR, the residual contribution to
! the partial derivative dP(T,RHO)/dT, with RHO held fixed,
! in MegaPascals/degrees Kelvin.
!
! Output, real ( kind = rk8 ) GR, the residual contribution to the Gibbs
! function, in KJ/kg.
!
! Output, real ( kind = rk8 ) HR, the residual contribution to the
! enthalpy, in KJ/kg.
!
! Output, real ( kind = rk8 ) PR, the residual contribution to the pressure,
! in MegaPascals.
!
! Output, real ( kind = rk8 ) SR, the residual contribution to the entropy,
! in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) UR, the residual contribution to the
! internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter, dimension ( 4 ) :: aad = (/ &
34.0D+00, 40.0D+00, 30.0D+00, 1050.0D+00 /)
real ( kind = rk8 ), parameter, dimension ( 4 ) :: aat = (/ &
20000.0D+00, 20000.0D+00, 40000.0D+00, 25.0D+00 /)
real ( kind = rk8 ), parameter, dimension ( 4 ) :: adz = (/ &
0.319D+00, 0.319D+00, 0.319D+00, 1.55D+00 /)
real ( kind = rk8 ) ar
real ( kind = rk8 ) att
real ( kind = rk8 ), parameter, dimension ( 4 ) :: atz = (/ &
640.0D+00, 640.0D+00, 641.6D+00, 270.0D+00 /)
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dadt
real ( kind = rk8 ) ddz
real ( kind = rk8 ) del
real ( kind = rk8 ) dex
real ( kind = rk8 ) dfdt
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) e
real ( kind = rk8 ) errtol
real ( kind = rk8 ) ex0
real ( kind = rk8 ) ex1
real ( kind = rk8 ) ex2
real ( kind = rk8 ) fct
real ( kind = rk8 ), parameter, dimension ( 40 ) :: g = (/ &
-530.62968529023D+00, 0.22744901424408D+04, 0.78779333020687D+03, &
-69.830527374994D+00, 0.17863832875422D+05,-0.39514731563338D+05, &
0.33803884280753D+05, -0.13855050202703D+05,-0.25637436613260D+06, &
0.48212575981415D+06, -0.34183016969660D+06, 0.12223156417448D+06, &
0.11797433655832D+07, -0.21734810110373D+07, 0.10829952168620D+07, &
-0.25441998064049D+06, -0.31377774947767D+07, 0.52911910757704D+07, &
-0.13802577177877D+07, -0.25109914369001D+06, 0.46561826115608D+07, &
-0.72752773275387D+07, 0.41774246148294D+06, 0.14016358244614D+07, &
-0.31555231392127D+07, 0.47929666384584D+07, 0.40912664781209D+06, &
-0.13626369388386D+07, 0.69625220862664D+06,-0.10834900096447D+07, &
-0.22722827401688D+06, 0.38365486000660D+06, 0.68833257944332D+04, &
0.21757245522644D+05, -0.26627944829770D+04,-0.70730418082074D+05, &
-0.225D+00, -1.68D+00, 0.055D+00, -93.0D+00 /)
real ( kind = rk8 ) gascon
real ( kind = rk8 ) gr
real ( kind = rk8 ) hr
integer i
integer, parameter, dimension ( 40 ) :: ii = (/ &
0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6, &
8,8,8,8,2,2,0,4,2,2,2,4 /)
integer j
integer, parameter, dimension ( 40 ) :: jj = (/ &
2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,&
2,3,5,7,1,4,4,4,0,2,0,0 /)
integer k
integer l
integer nc
real ( kind = rk8 ) pr
real ( kind = rk8 ) q10
real ( kind = rk8 ) q20
real ( kind = rk8 ) q2a
real ( kind = rk8 ) q5t
real ( kind = rk8 ) qm
real ( kind = rk8 ) qp
real ( kind = rk8 ) qr(11)
real ( kind = rk8 ) qt(10)
real ( kind = rk8 ) rho
real ( kind = rk8 ) sr
real ( kind = rk8 ), parameter :: s_ref = 7.6180720166752D+00
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.073D+00
real ( kind = rk8 ) tau
real ( kind = rk8 ) tx
real ( kind = rk8 ), parameter :: u_ref = - 4328.4549774261D+00
real ( kind = rk8 ) ur
real ( kind = rk8 ) v
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RESID - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'RESID - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
nc = 36
dpdrr = 0.0D+00
pr = 0.0D+00
ar = 0.0D+00
dadt = 0.0D+00
cvr = 0.0D+00
dpdtr = 0.0D+00
ex0 = - rho
ex0 = max ( ex0, -225.0D+00 )
ex0 = min ( ex0, 225.0D+00 )
e = exp ( ex0 )
q10 = rho * rho * e
q20 = 1.0D+00 - e
qr(1) = 0.0D+00
qr(2) = q10
do i = 2, 10
qr(i+1) = qr(i) * q20
end do
v = t_ref / t
qt(1) = t / t_ref
do i = 2, 10
qt(i) = qt(i-1) * v
end do
do i = 1, nc
k = ii(i) + 1
l = jj(i)
qp = g(i) * qr(k+1) * qt(l+1)
pr = pr + qp
dpdrr = dpdrr + ( 2.0D+00 / rho - ( 1.0D+00 - e * dble ( k - 1 ) / &
( 1.0D+00 - e ) ) ) * qp
ar = ar + g(i) * qr(k+2) * qt(l+1) / ( rho**2 * e * dble ( k ) &
* gascon ( ) * t )
dfdt = ( 1.0D+00 - e )**k * dble ( 1 - l ) * qt(l+2) / t_ref / dble ( k )
dadt = dadt + g(i) * dfdt
dpdtr = dpdtr + g(i) * dfdt * rho**2 * e * dble ( k ) / ( 1.0D+00 - e )
cvr = cvr + g(i) * dble ( l ) * dfdt / gascon()
end do
qp = 0.0D+00
q2a = 0.0D+00
do j = 37, 40
k = ii(j)
ddz = adz(j-36)
del = rho / ddz - 1.0D+00
if ( abs ( del ) < errtol ) then
del = errtol
end if
ex1 = - aad(j-36) * del**k
ex1 = max ( ex1, -225.0D+00 )
ex1 = min ( ex1, 225.0D+00 )
dex = exp ( ex1 ) * del**jj(j)
att = aat(j-36)
tx = atz(j-36)
tau = ( t / tx ) - 1.0D+00
ex2 = - att * tau * tau
ex2 = max ( ex2, -225.0D+00 )
ex2 = min ( ex2, 225.0D+00 )
q10 = dex * exp ( ex2 )
qm = dble ( jj(j) ) / del - dble ( k ) * aad(j-36) * del**(k-1)
fct = qm * rho**2 * q10 / ddz
q5t = fct * ( 2.0D+00 / rho + qm / ddz ) - ( rho / ddz )**2 * q10 * &
( dble ( jj(j) ) / del**2 + dble ( k * ( k - 1 ) ) * aad(j-36) * &
del**(k-2) )
dpdrr = dpdrr + q5t * g(j)
qp = qp + g(j) * fct
dadt = dadt - 2.0D+00 * g(j) * att * tau * q10 / tx
dpdtr = dpdtr - 2.0D+00 * g(j) * att * tau * fct / tx
q2a = q2a + t * g(j) * att * ( 4.0D+00 * ex2 + 2.0D+00 ) * q10 / tx**2
ar = ar + q10 * g(j) / ( gascon() * t )
end do
cvr = cvr + q2a / gascon()
pr = pr + qp
sr = - dadt / gascon()
ur = ar + sr
!
! Assign dimensions.
!
ar = gascon() * t * ar
cvr = gascon() * cvr
sr = gascon() * sr
ur = gascon() * t * ur
!
! Adjust energies.
!
ar = ar + gascon ( ) * t * s_ref - gascon ( ) * u_ref
sr = sr - gascon ( ) * s_ref
ur = ur - gascon ( ) * u_ref
gr = ar + pr / rho
hr = ur + pr / rho
return
end
subroutine secvir ( t, vir )
!*****************************************************************************80
!
!! secvir() calculates the second virial coefficient at a given temperature.
!
! Discussion:
!
! The second virial coefficient VIR appears in the first correction term
! to the ideal gas equation of state:
!
! P = R * T / volume + VIR / volume^2 + ...
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 28 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) VIR, the second virial coefficient, in CM3/G.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) b1
real ( kind = rk8 ) b1t
real ( kind = rk8 ) b1tt
real ( kind = rk8 ) b2
real ( kind = rk8 ) b2t
real ( kind = rk8 ) b2tt
real ( kind = rk8 ), parameter, dimension ( 5 ) :: g = (/ &
-0.53062968529023D+03, 0.22744901424408D+04, -0.26627944829770D+04, &
0.78779333020687D+03, -0.69830527374994D+02 /)
real ( kind = rk8 ) gascon
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.073D+00
real ( kind = rk8 ) v
real ( kind = rk8 ) vir
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SECVIR - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
call bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
v = t_ref / t
vir = b2 + ( &
v * ( g(1) &
+ v * ( g(2) &
+ v * ( g(3) &
+ v * ( g(4) &
+ v * v * g(5) ))))) &
/ ( gascon ( ) * t )
return
end
subroutine secvir_values ( n, tc, vir )
!*****************************************************************************80
!
!! secvir_values() returns some values of the second virial coefficient.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 03 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 24-25.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) VIR, the second virial coefficient, in
! m^3/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 19
integer n
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 5.0D+00, 10.0D+00, 20.0D+00, 30.0D+00, &
40.0D+00, 60.0D+00, 90.0D+00, 120.0D+00, 150.0D+00, &
180.0D+00, 210.0D+00, 240.0D+00, 300.0D+00, 400.0D+00, &
500.0D+00, 700.0D+00, 1000.0D+00, 2000.0D+00 /)
real ( kind = rk8 ) vir
real ( kind = rk8 ), save, dimension ( n_data ) :: virvec = (/ &
-98.96D+00, -90.08D+00, -82.29D+00, -69.36D+00, -59.19D+00, &
-51.07D+00, -39.13D+00, -27.81D+00, -20.83D+00, -16.21D+00, &
-12.98D+00, -10.63D+00, -8.85D+00, -6.39D+00, -4.03D+00, &
-2.71D+00, -1.32D+00, -0.39D+00, 0.53D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
vir = 0.0D+00
else
n = n + 1
tc = tcvec(n)
vir = virvec(n)
end if
return
end
subroutine sound ( t, p, c )
!*****************************************************************************80
!
!! sound() computes the speed of sound given temperature and pressure.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 22 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) C, the speed of sound, in m/s.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ) ai
real ( kind = rk8 ) ar
real ( kind = rk8 ) c
real ( kind = rk8 ) cp
real ( kind = rk8 ) cpi
real ( kind = rk8 ) cv
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvi
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) gb
real ( kind = rk8 ) gi
real ( kind = rk8 ) gr
real ( kind = rk8 ) hb
real ( kind = rk8 ) hi
real ( kind = rk8 ) hr
real ( kind = rk8 ) p
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rho_start
real ( kind = rk8 ) rhov
real ( kind = rk8 ) sb
real ( kind = rk8 ) si
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) t2
real ( kind = rk8 ) ub
real ( kind = rk8 ) ui
real ( kind = rk8 ) ur
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SOUND - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative pressures.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SOUND - Fatal error!'
write ( *, '(a)' ) ' The input pressure P must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was P = ', p
stop
end if
!
! For the given pressure, compute the saturation temperature.
!
rhol = 0.0D+00
rhov = 0.0D+00
call tsat ( p, t2, rhol, rhov )
!
! Depending on whether the temperature is above or below the
! saturation temperature, we expect to compute the density of
! a liquid or vapor.
!
if ( t < t2 ) then
rho_start = 1.9D+00
else
rho_start = 0.01D+00
end if
call dense ( p, t, rho_start, rho, dpdr )
!
! From T and RHO, compute the thermodynamic properties.
!
call ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
cv = cvb + cvr + cvi
dpdr = dpdrb + dpdrr
dpdt = dpdtb + dpdtr
cp = cv + t * dpdt * dpdt / ( dpdr * rho**2 )
c = sqrt ( 1000.0D+00 * cp * dpdr / cv )
return
end
subroutine sound_values ( n, tc, p, c )
!*****************************************************************************80
!
!! sound_values() returns some values of the speed of sound.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 238-246.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) C, the speed of sound, in m/s.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 20
real ( kind = rk8 ) c
real ( kind = rk8 ), save, dimension ( n_data ) :: cvec = (/ &
1401.0D+00, 472.8D+00, 533.7D+00, 585.7D+00, 609.5D+00, &
632.2D+00, 674.6D+00, 713.9D+00, 802.0D+00, 880.1D+00, &
1017.8D+00, 1115.9D+00, 1401.7D+00,1402.6D+00, 1409.6D+00, &
1418.1D+00, 1443.1D+00, 1484.6D+00, 1577.1D+00, 1913.4D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 5.0D+00, 10.0D+00, 50.0D+00, &
100.0D+00, 250.0D+00, 500.0D+00, 1000.0D+00, 2500.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 100.0D+00, 200.0D+00, 300.0D+00, 350.0D+00, &
400.0D+00, 500.0D+00, 600.0D+00, 850.0D+00, 1100.0D+00, &
1600.0D+00, 2000.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
c = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
c = cvec(n)
end if
return
end
subroutine surten ( t, sigma )
!*****************************************************************************80
!
!! surten() returns the surface tension as a function of temperature.
!
! Discussion:
!
! SURTEN uses an equation that yields values of the surface tension to
! within the accuracy of measurements from the triple point to the
! critical point.
!
! Sigma = B * ( (TSTAR-T)/TSTAR)^Mu * (1+b*(TSTAR-T)/TSTAR)
!
! where:
!
! TSTAR = 647.15 Degrees Kelvin,
! B = 0.2358 Pascals * Meters
! b = -0.625,
! Mu = 1.256.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) SIGMA, the surface tension,
! in Pascal * m = Newton / m.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter :: b_cap = 0.2358D+00
real ( kind = rk8 ), parameter :: b_small = -0.625D+00
real ( kind = rk8 ), parameter :: mu = 1.256D+00
real ( kind = rk8 ) sigma
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_star = 647.15D+00
real ( kind = rk8 ) term
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'SURTEN - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
term = ( t_star - t ) / t_star
sigma = b_cap * term**mu * ( 1.0D+00 + b_small * term )
!
! Need this conversion to match the table, but justification is there none.
!
sigma = 1000.0D+00 * sigma
return
end
subroutine surten_values ( n, tc, sigma )
!*****************************************************************************80
!
!! surten_values() returns some values of the surface tension.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 267.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) SIGMA, the surface tension,
! in Pascal * m = Newton / m.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 14
integer n
real ( kind = rk8 ) sigma
real ( kind = rk8 ), save, dimension ( n_data ) :: sigmavec = (/ &
74.22D+00, 72.74D+00, 71.20D+00, 69.60D+00, 67.95D+00, &
58.92D+00, 48.75D+00, 37.68D+00, 26.05D+00, 14.37D+00, &
8.78D+00, 3.67D+00, 0.40D+00, 0.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
10.0D+00, 20.0D+00, 30.0D+00, 40.0D+00, 50.0D+00, &
100.0D+00, 150.0D+00, 200.0D+00, 250.0D+00, 300.0D+00, &
325.0D+00, 350.0D+00, 370.0D+00, 373.976D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
sigma = 0.0D+00
else
n = n + 1
tc = tcvec(n)
sigma = sigmavec(n)
end if
return
end
subroutine tdpsdt ( t, dp )
!*****************************************************************************80
!
!! tdpsdt() computes the quantity T * dP(Sat)/dT.
!
! Discussion:
!
! Here T is the temperature and P(Sat) is the vapor pressure.
! It is used by TSAT_EST and TSAT.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 13 January 2013
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) DP, the value T*(dP(Sat)/dT),
! in MegaPascals.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ), parameter, dimension ( 8 ) :: a = (/ &
-7.8889166D+00, 2.5514255D+00, -6.716169D+00, 33.239495D+00, &
-105.38479D+00, 174.35319D+00, -148.39348D+00, 48.631602D+00 /)
real ( kind = rk8 ) b
real ( kind = rk8 ) c
real ( kind = rk8 ) dp
integer i
real ( kind = rk8 ) q
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.25D+00
real ( kind = rk8 ) v
real ( kind = rk8 ) w
real ( kind = rk8 ) y
real ( kind = rk8 ) z
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TDPSDT - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
v = t / t_ref
w = 1.0D+00 - v
b = 0.0D+00
c = 0.0D+00
do i = 1, 8
z = dble ( i + 1 ) / 2.0D+00
y = a(i) * w**z
c = c + ( y / w ) * ( 0.5D+00 - 0.5D+00 * dble ( i ) - 1.0D+00 / v )
b = b + y
end do
q = b / v
dp = 22.093D+00 * exp ( q ) * c
return
end
subroutine thercon ( t, rho, lambda )
!*****************************************************************************80
!
!! thercon() calculates thermal conductivity for given temperature and density.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 20 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) LAMBDA, the thermal conductivity,
! in mW/(m degrees Kelvin).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ), parameter, dimension ( 0:3 ) :: acof = (/ &
2.02223D+00, 14.11166D+00, 5.25597D+00, -2.01870D+00 /)
real ( kind = rk8 ), parameter :: a_con = 18.66D+00
real ( kind = rk8 ) b(0:4,0:5)
real ( kind = rk8 ), parameter :: b_con = 1.00D+00
real ( kind = rk8 ), parameter :: c_con = 3.7711D-08
real ( kind = rk8 ) chi
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cv
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdr2
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) eta
real ( kind = rk8 ) g
real ( kind = rk8 ) h
integer i
integer j
real ( kind = rk8 ) lambda
real ( kind = rk8 ) lambda0
real ( kind = rk8 ) lambda_del
real ( kind = rk8 ), parameter :: omega = 0.4678D+00
real ( kind = rk8 ) p
real ( kind = rk8 ), parameter :: p_ref = 22.115D+00
real ( kind = rk8 ) power
real ( kind = rk8 ) rho
real ( kind = rk8 ), parameter :: rho_ref = 317.763D+00
real ( kind = rk8 ) rho2
real ( kind = rk8 ) s
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_ref = 647.27D+00
real ( kind = rk8 ) total
real ( kind = rk8 ) u
data b / &
1.3293046D+00, 1.7018363D+00, 5.2246158D+00, &
8.7127675D+00, -1.8525999D+00, &
-0.40452437D+00, -2.2156845D+00, -10.124111D+00, &
-9.5000611D+00, 0.93404690D+00, &
0.24409490D+00, 1.6511057D+00, 4.9874687D+00, &
4.3786606D+00, 0.0D+00, &
0.018660751D+00, -0.76736002D+00, -0.27297694D+00, &
-0.91783782D+00, 0.0D+00, &
-0.12961068D+00, 0.37283344D+00, -0.43083393D+00, &
0.0D+00, 0.0D+00, &
0.044809953D+00, -0.11203160D+00, 0.13333849D+00, &
0.0D+00, 0.0D+00 /
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERCON - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERCON - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
!
! Compute DPDR, DPDT, ETA.
!
call therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
call viscosity ( t, rho, eta )
!
! Convert RHO from G/CM3 to kg/M3,
! Convert DPDR from ? to ?.
!
rho2 = 1000.0D+00 * rho
dpdr2 = dpdr / 1000.0D+00
!
! Compute LAMBDA0.
!
total = 0.0D+00
do i = 0, 3
total = total + acof(i) * ( t_ref / t )**i
end do
lambda0 = sqrt ( t / t_ref ) / total
!
! Compute CHI.
!
chi = rho2 * p_ref / ( rho_ref**2 * dpdr2 )
!
! Compute delta_Lambda.
!
power = - a_con * ( ( t_ref - t ) / t )**2 - b_con * ( ( rho2 - rho_ref ) &
/ rho_ref )**4
lambda_del = ( c_con / eta ) * ( ( t * rho_ref ) / ( t_ref * rho ) )**2 &
* ( t_ref / p_ref )**2 * dpdt**2 * chi**omega * sqrt ( rho2 / rho_ref ) &
* exp ( power )
!
! Compute LAMBDA.
!
total = 0.0D+00
do i = 0, 4
do j = 0, 5
total = total + b(i,j) * ( ( t_ref - t ) / t )**i * &
( ( rho2 - rho_ref ) / rho_ref )**j
end do
end do
lambda = lambda0 * exp ( ( rho2 / rho_ref ) * total ) + lambda_del
!
! Temporary fix.
!
lambda = 1000.0D+00 * lambda
return
end
subroutine thercon_values ( n, tc, p, lambda )
!*****************************************************************************80
!
!! thercon_values() returns some values of the thermal conductivity.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 264.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) LAMBDA, the thermal conductivity, in
! mW/(m degrees Kelvin).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 35
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ) lambda
real ( kind = rk8 ), save, dimension ( n_data ) :: lambdavec = (/ &
561.0D+00, 561.3D+00, 561.5D+00, 562.4D+00, 563.7D+00, &
565.1D+00, 566.5D+00, 567.9D+00, 569.3D+00, 570.6D+00, &
572.0D+00, 573.4D+00, 574.8D+00, 576.1D+00, 577.5D+00, &
580.2D+00, 582.9D+00, 585.5D+00, 588.1D+00, 590.7D+00, &
593.3D+00, 595.8D+00, 598.3D+00, 603.1D+00, 607.8D+00, &
612.2D+00, 607.2D+00, 643.6D+00, 666.8D+00, 25.08D+00, &
28.85D+00, 33.28D+00, 54.76D+00, 79.89D+00, 107.3D+00 /)
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 5.0D+00, 10.0D+00, 25.0D+00, 50.0D+00, &
75.0D+00, 100.0D+00, 125.0D+00, 150.0D+00, 175.0D+00, &
200.0D+00, 225.0D+00, 250.0D+00, 275.0D+00, 300.0D+00, &
350.0D+00, 400.0D+00, 450.0D+00, 500.0D+00, 550.0D+00, &
600.0D+00, 650.0D+00, 700.0D+00, 800.0D+00, 900.0D+00, &
1000.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, &
150.0D+00, 200.0D+00, 400.0D+00, 600.0D+00, 800.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
lambda = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
lambda = lambdavec(n)
end if
return
end
subroutine therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
!*****************************************************************************80
!
!! therm() calculates thermodynamic functions given temperature and density.
!
! Discussion:
!
! Thermodynamic values were calculated from an analytic equation
! that approximates the Helmholtz function (specific Helmholtz
! energy) for ordinary water and steam, of the form A=A(rho,T)
! where A is the Helmholtz function, rho the density, and T
! the absolute (thermodynamic) temperature. Any thermodynamic
! value for any state, liquid, vapor or metastable, may be
! calculated by differentiation of this equation in accord with
! the first and second laws of thermodynamics.
!
! The International Association for the Properties of Steam
! has provisionally accepted this formulation for the range
! 273.15 <= T <= 1273.15 degrees Kelvin, where, for 423.15 <= T,
! the maximum pressure is Pmax = 1500 MPa = 15000 bar, and for
! 273.15 <= T < 423.15, the maximum pressure is
! Pmax = 100 * (5 + (T-273.15)/15) MPa.
!
! Close to the critical point, a small region is excluded:
! Abs(T-Tk) < 1, abs((rho-rhok)/rhok) < 0.3.
!
! The equation has a wider useful range, namely, fluid states
! of pure, undissociated water and steam defined by
! 260 <= T <= 2500 K and 0 <= P <= 3000 MPa.
!
! Thermodynamic property values for specific volume, density,
! specific internal energy, specific enthalpy, and specific
! entropy of water and steam were tabulated over the range
! 0 <= t <= 2000 C, 0 <= P <= 3000 MPa. The reference
! state is the liquid at the triple point, for which the
! internal energy and entropy have been assigned the value zero.
!
! Thermodynamic quantities are determined from the Helmholtz function
! A(rho,T), which is computed as the sum of three terms:
!
! A(rho,T) = A(base)(rho,T) + A(residual)(rho,T) + A(ideal)(T)
! (Equation 1)
!
! Because A(rho,T) is everywhere single valued and analytic,
! we can derive closed form relations for all other properties.
! In the following, unless otherwise indicated, the independent
! variables are temperature T and density RHO, and differentiation
! with respect to one variable is to imply that the other is fixed.
!
! Pressure: P = RHO^2 * dA/dRHO
! Density derivative: dP/dRHO = 2*P/RHO + RHO^2 * d2A/dRHO2
! Temperature derivative: dP/dT = RHO^2 * d2A/(dRHO dT)
! Specific entropy: S = - dA/dT
! Specific internal energy: U = A + T*S
! Specific enthalpy: H = U + P/RHO
! Specific heat capacity
! at constant volume: Cv = - T * d2A/dT2
! Specific Gibbs function: G = A + P/RHO
! Specific heat capacity
! at constant pressure: Cp = Cv + (T*(dP/dT)^2)/(RHO^2*dP/dRHO)
! Speed of sound: Omega = Sqrt ((Cp/Cv) * dP/dRHO)
! Second virial coefficient: B = 1/(2*R*T) * (d2P/dRHO2) (at RHO=0)
! Isothermal Joule-Thomson
! coefficient: DeltaT = (dH/dP) (fixed T) =
! (1/RHO)-(T*dP/dT)/(RHO^2*dP/dRHO)
! Joule-Thomson coefficient: Mu = (dT/dP) (fixed H) = DeltaT/Cp
! Isentropic temperature-
! pressure coefficient: BetaS = (dT/dP) (fixed S) =
! (DeltaT - 1/RHO)/Cp
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 19 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the fluid density, in G/CM3.
!
! Output, real ( kind = rk8 ) A, the Helmholtz function, in KJ/kg.
!
! Output, real ( kind = rk8 ) CJTH, the Joule-Thomson coefficient,
! in K/MegaPascals.
!
! Output, real ( kind = rk8 ) CJTT, the isothermal Joule-Thomson coefficient,
! in CM3/G.
!
! Output, real ( kind = rk8 ) CP, the isobaric (constant pressure) heat
! capacity, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) CV, the isochoric (constant volume) heat
! capacity, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) DPDR, the partial derivative
! dP(T,RHO)/dRHO, with T held fixed, in MegaPascals*CM3/G.
!
! Output, real ( kind = rk8 ) DPDT, the partial derivative
! dP(T,RHO)/dT, with RHO held fixed, in MegaPascals/degrees Kelvin.
!
! Output, real ( kind = rk8 ) G, the Gibbs free energy, in KJ/kg.
!
! Output, real ( kind = rk8 ) H, the enthalpy, in KJ/kg.
!
! Output, real ( kind = rk8 ) P, the pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) S, the entropy, in KJ/(kg degrees Kelvin).
!
! Output, real ( kind = rk8 ) U, the internal energy, in KJ/kg.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) a
real ( kind = rk8 ) ab
real ( kind = rk8 ) ai
real ( kind = rk8 ) ar
real ( kind = rk8 ) cjth
real ( kind = rk8 ) cjtt
real ( kind = rk8 ) cp
real ( kind = rk8 ) cpi
real ( kind = rk8 ) cv
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvi
real ( kind = rk8 ) cvr
logical, parameter :: debug = .false.
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) g
real ( kind = rk8 ) gb
real ( kind = rk8 ) gi
real ( kind = rk8 ) gr
real ( kind = rk8 ) h
real ( kind = rk8 ) hb
real ( kind = rk8 ) hi
real ( kind = rk8 ) hr
real ( kind = rk8 ) p
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) s
real ( kind = rk8 ) sb
real ( kind = rk8 ) si
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) u
real ( kind = rk8 ) ub
real ( kind = rk8 ) ui
real ( kind = rk8 ) ur
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERM - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'THERM - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' Input value was RHO = ', rho
stop
end if
call ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
a = ab + ar + ai
cv = cvb + cvr + cvi
if ( debug ) then
write ( *, * ) ' '
write ( *, * ) 'THERM:'
write ( *, * ) ' CVB = ', cvb
write ( *, * ) ' CVR = ', cvr
write ( *, * ) ' CVI = ', cvi
write ( *, * ) ' CV = ', cv
end if
dpdr = dpdrb + dpdrr
dpdt = dpdtb + dpdtr
p = pb + pr
s = sb + sr + si
u = ub + ur + ui
if ( debug ) then
write ( *, * ) ' '
write ( *, * ) 'THERM:'
write ( *, * ) ' UB = ', ub
write ( *, * ) ' UR = ', ur
write ( *, * ) ' UI = ', ui
end if
g = a + p / rho
h = u + p / rho
cp = cv + t * dpdt**2 / ( dpdr * rho**2 )
cjtt = 1.0D+00 / rho - t * dpdt / ( dpdr * rho**2 )
cjth = - cjtt / cp
return
end
subroutine tsat ( p, t, rhol, rhov )
!*****************************************************************************80
!
!! tsat() calculates the saturation temperature for a given pressure.
!
! Discussion:
!
! The corresponding liquid and vapor densities are also computed.
! The saturation temperature is also known as the "vapor temperature".
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) P, the vapor pressure, in MegaPascals.
!
! Output, real ( kind = rk8 ) T, the vapor temperature, in degrees Kelvin.
!
! Output, real ( kind = rk8 ) RHOL, the liquid density, in G/CM3.
!
! Output, real ( kind = rk8 ) RHOV, the vapor density, in G/CM3.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
logical, parameter :: debug = .false.
real ( kind = rk8 ) delg
real ( kind = rk8 ) dp
real ( kind = rk8 ) dp2
real ( kind = rk8 ) errtol
real ( kind = rk8 ) gascon
integer it
integer, parameter :: it_max = 50
real ( kind = rk8 ) p
real ( kind = rk8 ) p_consistent
real ( kind = rk8 ) rhol
real ( kind = rk8 ) rhov
real ( kind = rk8 ) t
!
! Initialize output quantities, obliterating any input value.
!
t = 0.0D+00
rhol = 0.0D+00
rhov = 0.0D+00
!
! Set the error tolerance.
!
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative pressure.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT - Fatal error!'
write ( *, '(a)' ) ' The input pressure must be positive!'
write ( *, '(a,g14.6)' ) ' Your value was P = ', p
stop
end if
!
! Estimate the saturation temperature.
!
call tsat_est ( p, t )
if ( debug ) then
write ( *, * ) ' '
write ( *, * ) 'TSAT:'
write ( *, '(2g14.6)' ) p, t, rhol, rhov
end if
do it = 1, it_max
call corr ( t, p, p_consistent, rhol, rhov, delg )
dp = delg * gascon ( ) * t * rhol * rhov / ( rhol - rhov )
call tdpsdt ( t, dp2 )
t = t * ( 1.0D+00 - dp / dp2 )
if ( debug ) then
write ( *, '(2g14.6)' ) p, t, rhol, rhov
end if
if ( abs ( delg ) < errtol ) then
return
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' Number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Last iterate was ', t
write ( *, '(a,g14.6)' ) ' Last DELG was ', delg
return
end
subroutine tsat_est ( p, t )
!*****************************************************************************80
!
!! tsat_est() makes a rough estimate of the saturation temperature.
!
! Discussion:
!
! The saturation temperature is also called the vapor temperature.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) P, the pressure, in MegaPascals. The tabulated
! range for P is
! 0.00061173 MegaPascals <= P <= P_CRIT = 22.055 MegaPascals.
! The input value of P must be positive.
!
! Output, real ( kind = rk8 ) T, the saturation temperature,
! in degrees Kelvin. This value will always be in the range
! [ 273.15, 647.126 ].
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: npol = 4
real ( kind = rk8 ), parameter, dimension ( 0:npol ) :: c = (/ &
372.83D+00, 27.7589D+00, 2.3819D+00, 0.24834D+00, 0.0193855D+00 /)
real ( kind = rk8 ) dp
real ( kind = rk8 ) dt
real ( kind = rk8 ) errtol
integer it
integer, parameter :: it_max = 8
real ( kind = rk8 ) p
real ( kind = rk8 ) pl
real ( kind = rk8 ), parameter :: p_crit = 22.055D+00
real ( kind = rk8 ) pp
real ( kind = rk8 ) t
real ( kind = rk8 ), parameter :: t_crit = 647.126D+00
real ( kind = rk8 ), parameter :: t_min = 273.15D+00
real ( kind = rk8 ) t_old
errtol = sqrt ( epsilon ( errtol ) )
!
! Refuse to handle zero or negative pressure.
!
if ( p <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT_EST - Fatal error!'
write ( *, '(a)' ) ' The input pressure must be positive!'
write ( *, '(a,g14.6)' ) ' Your value was P = ', p
stop
end if
if ( p_crit < p ) then
t = t_crit
return
end if
!
! The initial estimate for T uses a polyonmial in the logarithm of P.
!
pl = 2.302585D+00 + log ( p )
call r8poly_val_horner ( npol, c, pl, t )
t = min ( t, t_crit )
t = max ( t, t_min )
dt = 0.0D+00
do it = 1, it_max
call psat_est ( t, pp )
call tdpsdt ( t, dp )
if ( abs ( p - pp ) < errtol * p ) then
return
end if
dt = t * ( p - pp ) / dp
t_old = t
t = t * ( 1.0D+00 + ( p - pp ) / dp )
t = min ( t, t_crit )
t = max ( t, t_min )
if ( abs ( dt ) < errtol * ( abs ( t ) + 1.0D+00 ) ) then
return
else if ( abs ( t - t_old ) < errtol ) then
return
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TSAT_EST - Warning!'
write ( *, '(a)' ) ' The iteration did not converge.'
write ( *, '(a,i6)' ) ' Number of iterations was ', it_max
write ( *, '(a,g14.6)' ) ' Convergence tolerance was ', errtol
write ( *, '(a,g14.6)' ) ' Last iterate was ', t
write ( *, '(a,g14.6)' ) ' Last correction was ', dt
return
end
subroutine tsat_values ( n, p, tc )
!*****************************************************************************80
!
!! tsat_values() returns some values of the saturation temperature.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 05 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, pages 16-22.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) TC, the saturation temperature, in
! degrees Celsius.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 20
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
0.0061173D+00, 0.012D+00, 0.025D+00, 0.055D+00, 0.080D+00, &
0.11D+00, 0.16D+00, 0.25D+00, 0.50D+00, 0.75D+00, &
1.0D+00, 1.5D+00, 2.0D+00, 5.0D+00, 10.0D+00, &
20.0D+00, 50.0D+00, 100.0D+00, 200.0D+00, 220.55D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.010D+00, 9.655D+00, 21.080D+00, 34.589D+00, 41.518D+00, &
47.695D+00, 55.327D+00, 64.980D+00, 81.339D+00, 91.783D+00, &
99.632D+00, 111.378D+00, 120.443D+00, 151.866D+00, 179.916D+00, &
212.417D+00, 263.977D+00, 311.031D+00, 365.800D+00, 373.976D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
p = 0.0D+00
tc = 0.0D+00
else
n = n + 1
p = pvec(n)
tc = tcvec(n)
end if
return
end
subroutine viscosity ( t, rho, eta )
!*****************************************************************************80
!
!! viscosity() calculates the viscosity for given temperature and density.
!
! Discussion:
!
! On 02 February 2002, I discovered that the Haar/Gallagher/Kell
! reference apparently reversed the sign on the A3 coefficient.
! That made the results better, but still off.
!
! Apparently Haar/Gallagher/Kell had a transcription error in
! the value of B(4,1), which they list as -0.273093, but which
! should be -0.253093.
!
! These two corrections courtesy of Meyer/McClintock/Silvestri/Spencer.
!
! Now the results look proper! And just 12 years late...
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 02 February 2002
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! International Association for the Properties of Steam,
! Release on Dynamic Viscosity of Water Substance,
! National Bureau of Standards, Washington DC, 1975, revised 1983.
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the density, in G/CM3.
!
! Output, real ( kind = rk8 ) ETA, the viscosity, in MegaPascal seconds.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: npol_t = 3
real ( kind = rk8 ), parameter, dimension ( 0:npol_t ) :: a = (/ &
0.0181583D+00, 0.0177624D+00, 0.0105287D+00, -0.0036744D+00 /)
real ( kind = rk8 ) arg
real ( kind = rk8 ), dimension(0:5,0:4) :: b = reshape ( (/ &
0.501938D+00, 0.162888D+00, -0.130356D+00, &
0.907919D+00, -0.551119D+00, 0.146543D+00, &
0.235622D+00, 0.789393D+00, 0.673665D+00, &
1.207552D+00, 0.0670665D+00, -0.0843370D+00, &
-0.274637D+00, -0.743539D+00, -0.959456D+00, &
-0.687343D+00, -0.497089D+00, 0.195286D+00, &
0.145831D+00, 0.263129D+00, 0.347247D+00, &
0.213486D+00, 0.100754D+00, -0.032932D+00, &
-0.0270448D+00, -0.0253093D+00, -0.0267758D+00, &
-0.0822904D+00, 0.0602253D+00, -0.0202595D+00 /), &
(/ 6,5 /) )
logical, parameter :: debug = .false.
real ( kind = rk8 ) eta
real ( kind = rk8 ) eta0
integer i
integer j
real ( kind = rk8 ) rho
real ( kind = rk8 ), parameter :: rho_max = 1.050D+00
real ( kind = rk8 ), parameter :: rho_ref = 0.317763D+00
real ( kind = rk8 ) t
real ( kind = rk8 ) total
real ( kind = rk8 ), parameter :: t_max = 800.00D+00
real ( kind = rk8 ), parameter :: t_ref = 647.27D+00
logical, save :: warning1 = .false.
logical, save :: warning2 = .false.
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
if ( t_max < t .and. ( .not. warning1 ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Warning (once only)!'
write ( *, '(a,g14.6)' ) &
' The input temperature T should be no more than ', t_max
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
write ( *, '(a)' ) ' '
warning1 = .true.
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was RHO = ', rho
stop
end if
if ( rho_max < rho .and. ( .not. warning2 ) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VISCOSITY - Warning (once only)!'
write ( *, '(a,g14.6)' ) &
' The input density RHO should be no more than ', rho_max
write ( *, '(a,g14.6)' ) ' The input value was RHO = ', rho
write ( *, '(a)' ) ' '
warning2 = .true.
end if
!
! Compute ETA0.
!
arg = t_ref / t
call r8poly_val_horner ( npol_t, a, arg, total )
eta0 = sqrt ( t / t_ref ) / total
!
! Compute ETA.
!
total = 0.0D+00
do i = 0, 5
do j = 0, 4
total = total + b(i,j) * ( ( t_ref - t ) / t )**i * &
( ( rho - rho_ref ) / rho_ref )**j
end do
end do
eta = eta0 * exp ( ( rho / rho_ref ) * total )
return
end
subroutine viscosity_values ( n, tc, p, eta )
!*****************************************************************************80
!
!! viscosity_values() returns some values of the viscosity function for testing.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 04 February 2002
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3, page 263.
!
! Parameters:
!
! Input/output, integer N.
! On input, if N is 0, the first test data is returned, and N is set
! to the index of the test data. On each subsequent call, N is
! incremented and that test data is returned. When there is no more
! test data, N is set to 0.
!
! Output, real ( kind = rk8 ) TC, the temperature, in degrees Celsius.
!
! Output, real ( kind = rk8 ) P, the pressure, in bar.
!
! Output, real ( kind = rk8 ) ETA, the viscosity, in MegaPascal seconds.
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
integer, parameter :: n_data = 34
real ( kind = rk8 ) eta
real ( kind = rk8 ), save, dimension ( n_data ) :: etavec = (/ &
1792.0D+00, 1791.0D+00, 1790.0D+00, 1786.0D+00, 1780.0D+00, &
1775.0D+00, 1769.0D+00, 1764.0D+00, 1759.0D+00, 1754.0D+00, &
1749.0D+00, 1744.0D+00, 1739.0D+00, 1735.0D+00, 1731.0D+00, &
1722.0D+00, 1714.0D+00, 1707.0D+00, 1700.0D+00, 1694.0D+00, &
1687.0D+00, 1682.0D+00, 1676.0D+00, 1667.0D+00, 1659.0D+00, &
1653.0D+00, 890.8D+00, 547.1D+00, 378.4D+00, 12.28D+00, &
16.18D+00, 24.45D+00, 32.61D+00, 40.38D+00 /)
integer n
real ( kind = rk8 ) p
real ( kind = rk8 ), save, dimension ( n_data ) :: pvec = (/ &
1.0D+00, 5.0D+00, 10.0D+00, 25.0D+00, 50.0D+00, &
75.0D+00, 100.0D+00, 125.0D+00, 150.0D+00, 175.0D+00, &
200.0D+00, 225.0D+00, 250.0D+00, 275.0D+00, 300.0D+00, &
350.0D+00, 400.0D+00, 450.0D+00, 500.0D+00, 550.0D+00, &
600.0D+00, 650.0D+00, 700.0D+00, 800.0D+00, 900.0D+00, &
1000.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00, &
1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /)
real ( kind = rk8 ) tc
real ( kind = rk8 ), save, dimension ( n_data ) :: tcvec = (/ &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00, &
0.0D+00, 25.0D+00, 50.0D+00, 75.0D+00, 100.0D+00, &
200.0D+00, 400.0D+00, 600.0D+00, 800.0D+00 /)
if ( n < 0 ) then
n = 0
end if
if ( n_data <= n ) then
n = 0
tc = 0.0D+00
p = 0.0D+00
eta = 0.0D+00
else
n = n + 1
tc = tcvec(n)
p = pvec(n)
eta = etavec(n)
end if
return
end
subroutine volume ( t, rho, v, dvdt, dvdr )
!*****************************************************************************80
!
!! volume() computes specific volume derivatives given temperature and density.
!
! Discussion:
!
! Because A(rho,T) is everywhere single valued and analytic,
! we can derive closed form relations for all other properties.
!
! The independent variables are temperature T and density RHO,
! and differentiation with respect to one variable is to imply that
! the other is held at a fixed value.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 28 November 1998
!
! Author:
!
! Original Fortran77 version by Lester Haar, John Gallagher, George Kell.
! This version by John Burkardt.
!
! Reference:
!
! Lester Haar, John Gallagher, George Kell,
! NBS/NRC Steam Tables:
! Thermodynamic and Transport Properties and Computer Programs
! for Vapor and Liquid States of Water in SI Units,
! Hemisphere Publishing Corporation, Washington, 1984,
! LC: TJ270.H3
!
! C A Meyer, R B McClintock, G J Silvestri, R C Spencer,
! ASME Steam Tables: Thermodynamic and Transport Properties of Steam,
! American Society of Mechanical Engineers,
! Fifth Edition, 1983,
! LC: TJ270.A75.
!
! Parameters:
!
! Input, real ( kind = rk8 ) T, the temperature, in degrees Kelvin.
!
! Input, real ( kind = rk8 ) RHO, the fluid density, in G/CM3.
!
! Output, real ( kind = rk8 ) V, the specific volume, in CM3/G.
!
! Output, real ( kind = rk8 ) DVDT, the partial derivative dV(T,RHO)/dT,
! where V is the specific volume, in CM3 / (G * degrees Kelvin).
!
! Output, real ( kind = rk8 ) DVDR, the partial derivative dV(T,RHO)/dRHO,
! where V is the specific volume, in CM3^2 / ( G^2 ).
!
implicit none
integer, parameter :: rk8 = kind ( 1.0D+00 )
real ( kind = rk8 ) ab
real ( kind = rk8 ) ar
real ( kind = rk8 ) cvb
real ( kind = rk8 ) cvr
real ( kind = rk8 ) dpdr
real ( kind = rk8 ) dpdrb
real ( kind = rk8 ) dpdrr
real ( kind = rk8 ) dpdt
real ( kind = rk8 ) dpdtb
real ( kind = rk8 ) dpdtr
real ( kind = rk8 ) dvdr
real ( kind = rk8 ) dvdt
real ( kind = rk8 ) gb
real ( kind = rk8 ) gr
real ( kind = rk8 ) hb
real ( kind = rk8 ) hr
real ( kind = rk8 ) pb
real ( kind = rk8 ) pr
real ( kind = rk8 ) rho
real ( kind = rk8 ) sb
real ( kind = rk8 ) sr
real ( kind = rk8 ) t
real ( kind = rk8 ) ub
real ( kind = rk8 ) ur
real ( kind = rk8 ) v
!
! Refuse to handle zero or negative temperatures.
!
if ( t <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VOLUME - Fatal error!'
write ( *, '(a)' ) ' The input temperature T must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was T = ', t
stop
end if
!
! Refuse to handle zero or negative density.
!
if ( rho <= 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'VOLUME - Fatal error!'
write ( *, '(a)' ) ' The input density RHO must be positive.'
write ( *, '(a,g14.6)' ) ' The input value was RHO = ', rho
stop
end if
call resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
call base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
dpdr = dpdrb + dpdrr
dpdt = dpdtb + dpdtr
dvdt = dpdt / ( dpdr * rho**2 )
!
! Because V = 1/Rho, dV/dRho = -1/Rho**2
!
dvdr = - 1.0D+00 / rho**2
v = 1.0D+00 / rho
return
end
! =========================================================================
! Source File: steam_calc.f90
! =========================================================================
program steam_calc
implicit none
integer, parameter :: rk8 = kind(1.0D+00)
! Input variables
integer :: input_type
double precision :: val1, val2
! Conversions
double precision :: t_c, t_k, p_kpa, p_mpa
double precision :: x_qual, spec_vol, spec_u, spec_h, spec_s, density
double precision :: cp_val, cv_val, speed_sound, visc_eta, therm_lambda
! Saturated values
double precision :: t_sat_k, p_sat_mpa, rhol, rhov
double precision :: h_liq, h_vap, s_liq, s_vap, u_liq, u_vap
double precision :: cp_liq, cp_vap, cv_liq, cv_vap, sound_liq, sound_vap
double precision :: visc_liq, visc_vap, therm_liq, therm_vap
! Calculation values
character(len=30) :: state_name
double precision :: rho, rho_start, dpdr, dpdt
double precision :: ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub
double precision :: ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur
double precision :: a_term, cjth, cjtt, cp_term, cv_term, g_term, h_term, p_term, s_term, u_term
! Solver variables
double precision :: t_min, t_max, t_mid, h_mid, s_mid, err
integer :: iter
logical :: is_mix
! Functions from steam.f90
double precision :: gascon
! Read from stdin
read(*,*) input_type
read(*,*) val1
read(*,*) val2
! Default quality set to -1 (meaning not in mixture region, single phase)
x_qual = -1.0d0
is_mix = .false.
select case(input_type)
! -------------------------------------------------------------
! 1. Temperature and Pressure (val1 = T [C], val2 = P [kPa])
! -------------------------------------------------------------
case(1)
t_c = val1
t_k = t_c + 273.15d0
p_kpa = val2
p_mpa = p_kpa / 1000.0d0
! Check boundaries
if (t_k < 273.16d0 .or. t_k > 1273.15d0) then
write(*,*) 'ERROR: Temperature must be between 0.01 C and 1000 C'
stop
end if
if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
stop
end if
if (p_mpa < 22.055d0) then
call tsat(p_mpa, t_sat_k, rhol, rhov)
if (abs(t_k - t_sat_k) < 1.0d-3) then
! On the saturation boundary, assume superheated or mixture
! But T & P is not enough to define mixture state, so default to saturated liquid
state_name = 'Saturated boundary'
rho_start = rhol
else if (t_k < t_sat_k) then
state_name = 'Subcooled Liquid'
rho_start = 1.11d0 - 0.0004d0 * t_k
else
state_name = 'Superheated Vapor'
rho_start = p_mpa / (gascon() * t_k)
end if
else
state_name = 'Supercritical Fluid'
if (t_k < 647.13d0) then
rho_start = 1.11d0 - 0.0004d0 * t_k
else
rho_start = p_mpa / (gascon() * t_k)
end if
end if
call dense(p_mpa, t_k, rho_start, rho, dpdr)
call therm(t_k, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, spec_h, p_term, spec_s, spec_u)
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
! -------------------------------------------------------------
! 2. Temperature and Quality (val1 = T [C], val2 = x)
! -------------------------------------------------------------
case(2)
t_c = val1
t_k = t_c + 273.15d0
x_qual = val2
if (t_k < 273.16d0 .or. t_k >= 647.13d0) then
write(*,*) 'ERROR: Saturated state requires T between 0.01 C and 373.9 C'
stop
end if
if (x_qual < 0.0d0 .or. x_qual > 1.0d0) then
write(*,*) 'ERROR: Quality x must be between 0 and 1'
stop
end if
call psat(t_k, p_sat_mpa, rhol, rhov)
p_mpa = p_sat_mpa
p_kpa = p_mpa * 1000.0d0
is_mix = .true.
if (x_qual == 0.0d0) then
state_name = 'Saturated Liquid'
else if (x_qual == 1.0d0) then
state_name = 'Saturated Vapor'
else
state_name = 'Saturated Mixture'
end if
! Calculate properties of saturated liquid
call therm(t_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
! Calculate properties of saturated vapor
call therm(t_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
! Blend properties
spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
! -------------------------------------------------------------
! 3. Pressure and Quality (val1 = P [kPa], val2 = x)
! -------------------------------------------------------------
case(3)
p_kpa = val1
p_mpa = p_kpa / 1000.0d0
x_qual = val2
if (p_mpa < 0.000611d0 .or. p_mpa >= 22.055d0) then
write(*,*) 'ERROR: Saturated state requires P between 0.61 kPa and 22055 kPa'
stop
end if
if (x_qual < 0.0d0 .or. x_qual > 1.0d0) then
write(*,*) 'ERROR: Quality x must be between 0 and 1'
stop
end if
call tsat(p_mpa, t_sat_k, rhol, rhov)
t_k = t_sat_k
t_c = t_k - 273.15d0
is_mix = .true.
if (x_qual == 0.0d0) then
state_name = 'Saturated Liquid'
else if (x_qual == 1.0d0) then
state_name = 'Saturated Vapor'
else
state_name = 'Saturated Mixture'
end if
! Calculate properties of saturated liquid
call therm(t_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
! Calculate properties of saturated vapor
call therm(t_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
! Blend properties
spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
! -------------------------------------------------------------
! 4. Pressure and Enthalpy (val1 = P [kPa], val2 = h [kJ/kg])
! -------------------------------------------------------------
case(4)
p_kpa = val1
p_mpa = p_kpa / 1000.0d0
spec_h = val2
if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
stop
end if
if (p_mpa < 22.055d0) then
! Check if it falls in saturation boundary
call tsat(p_mpa, t_sat_k, rhol, rhov)
call therm(t_sat_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call therm(t_sat_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
if (spec_h >= h_liq .and. spec_h <= h_vap) then
! Saturated mixture!
x_qual = (spec_h - h_liq) / (h_vap - h_liq)
t_k = t_sat_k
t_c = t_k - 273.15d0
is_mix = .true.
state_name = 'Saturated Mixture'
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
spec_s = (1.0d0 - x_qual) * s_liq + x_qual * s_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
else if (spec_h < h_liq) then
! Subcooled liquid
state_name = 'Subcooled Liquid'
t_min = 273.16d0
t_max = t_sat_k
rho_start = 1.0d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (h_mid < spec_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
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
else
! Superheated vapor
state_name = 'Superheated Vapor'
t_min = t_sat_k
t_max = 1273.15d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
rho_start = p_mpa / (gascon() * t_mid)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (h_mid < spec_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
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
else
! Supercritical region
state_name = 'Supercritical Fluid'
t_min = 273.16d0
t_max = 1273.15d0
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
if (t_mid < 647.13d0) then
rho_start = 1.0d0
else
rho_start = p_mpa / (gascon() * t_mid)
end if
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (h_mid < spec_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
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
! -------------------------------------------------------------
! 5. Pressure and Entropy (val1 = P [kPa], val2 = s [kJ/(kg K)])
! -------------------------------------------------------------
case(5)
p_kpa = val1
p_mpa = p_kpa / 1000.0d0
spec_s = val2
if (p_mpa < 0.000611d0 .or. p_mpa > 100.0d0) then
write(*,*) 'ERROR: Pressure must be between 0.61 kPa and 100000 kPa'
stop
end if
if (p_mpa < 22.055d0) then
! Check if it falls in saturation boundary
call tsat(p_mpa, t_sat_k, rhol, rhov)
call therm(t_sat_k, rhol, a_term, cjth, cjtt, cp_liq, cv_liq, dpdr, dpdt, &
g_term, h_liq, p_term, s_liq, u_liq)
call therm(t_sat_k, rhov, a_term, cjth, cjtt, cp_vap, cv_vap, dpdr, dpdt, &
g_term, h_vap, p_term, s_vap, u_vap)
if (spec_s >= s_liq .and. spec_s <= s_vap) then
! Saturated mixture!
x_qual = (spec_s - s_liq) / (s_vap - s_liq)
t_k = t_sat_k
t_c = t_k - 273.15d0
is_mix = .true.
state_name = 'Saturated Mixture'
call viscosity(t_k, rhol, visc_liq)
call thercon(t_k, rhol, therm_liq)
call sound(t_k, p_mpa, sound_liq)
call viscosity(t_k, rhov, visc_vap)
call thercon(t_k, rhov, therm_vap)
call sound(t_k, p_mpa, sound_vap)
spec_h = (1.0d0 - x_qual) * h_liq + x_qual * h_vap
spec_u = (1.0d0 - x_qual) * u_liq + x_qual * u_vap
spec_vol = (1.0d0 - x_qual) / (rhol * 1000.0d0) + x_qual / (rhov * 1000.0d0)
density = 1.0d0 / (spec_vol * 1000.0d0)
cp_val = (1.0d0 - x_qual) * cp_liq + x_qual * cp_vap
cv_val = (1.0d0 - x_qual) * cv_liq + x_qual * cv_vap
speed_sound = (1.0d0 - x_qual) * sound_liq + x_qual * sound_vap
visc_eta = (1.0d0 - x_qual) * visc_liq + x_qual * visc_vap
therm_lambda = (1.0d0 - x_qual) * therm_liq + x_qual * therm_vap
rho = density
else if (spec_s < s_liq) then
! Subcooled liquid
state_name = 'Subcooled Liquid'
t_min = 273.16d0
t_max = t_sat_k
rho_start = 1.0d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (s_mid < spec_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
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
else
! Superheated vapor
state_name = 'Superheated Vapor'
t_min = t_sat_k
t_max = 1273.15d0
! Bisection solver for T
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
rho_start = p_mpa / (gascon() * t_mid)
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (s_mid < spec_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
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
else
! Supercritical region
state_name = 'Supercritical Fluid'
t_min = 273.16d0
t_max = 1273.15d0
do iter = 1, 40
t_mid = 0.5d0 * (t_min + t_max)
if (t_mid < 647.13d0) then
rho_start = 1.0d0
else
rho_start = p_mpa / (gascon() * t_mid)
end if
call dense(p_mpa, t_mid, rho_start, rho, dpdr)
call therm(t_mid, rho, a_term, cjth, cjtt, cp_val, cv_val, dpdr, dpdt, &
g_term, h_mid, p_term, s_mid, u_term)
if (s_mid < spec_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
spec_h = h_mid
spec_s = s_mid
spec_u = u_term
call viscosity(t_k, rho, visc_eta)
call thercon(t_k, rho, therm_lambda)
call sound(t_k, p_mpa, speed_sound)
end if
end select
! Density conversion: rho is in g/cm3. Output specific volume v in m3/kg.
if (.not. is_mix) then
density = rho * 1000.0d0 ! kg/m3
spec_vol = 1.0d0 / density ! m3/kg
end if
! If it's single phase, compute dynamic viscosity in Pa s
! visc_eta is in microPascal-seconds. Let's convert to Pa s: * 1.0E-6
visc_eta = visc_eta * 1.0d-6
! Convert thermal conductivity to W/(m K): therm_lambda is in mW/(m K) => / 1000.d0
therm_lambda = therm_lambda / 1000.d0
! Convert quality
if (.not. is_mix) then
if (state_name == 'Subcooled Liquid') then
x_qual = 0.0d0
else if (state_name == 'Superheated Vapor') then
x_qual = 1.0d0
else
x_qual = -1.0d0 ! Supercritical has no quality defined
end if
end if
! Print output in parsable format
write(*,'(A,A)') 'STATE=', trim(state_name)
write(*,'(A,F14.4)') 'T=', t_c
write(*,'(A,F14.4)') 'P=', p_kpa
write(*,'(A,F14.6)') 'V=', spec_vol
write(*,'(A,F14.4)') 'H=', spec_h
write(*,'(A,F14.6)') 'S=', spec_s
write(*,'(A,F14.4)') 'U=', spec_u
write(*,'(A,F14.4)') 'CP=', cp_val
write(*,'(A,F14.4)') 'CV=', cv_val
write(*,'(A,F14.2)') 'SOUND=', speed_sound
write(*,'(A,E14.6)') 'VISC=', visc_eta
write(*,'(A,F14.6)') 'THERM=', therm_lambda
write(*,'(A,F14.4)') 'QUALITY=', x_qual
end program steam_calc
Solver Description
Thermodynamic properties of water and steam are modeled using the IAPWS-IF97 industrial formulation, which divides the state space into 5 regions, mapping Helmholtz free energy equations $a(\rho, T)$ and Gibbs free energy equations $g(P, T)$ to obtain pressure, volume, temperature, enthalpy, and entropy.
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 steam.f90 steam_calc.f90 -o steam_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
steam_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Input combination mode (1=T/P, 2=T/x, 3=P/x, 4=P/h, 5=P/s)
1
! Parameter 1 (Temperature [°C] or Pressure [kPa])
150.0
! Parameter 2 (Pressure [kPa] or Quality [0-1] or Enthalpy [kJ/kg] or Entropy [kJ/kg-K])
101.325
1
! Parameter 1 (Temperature [°C] or Pressure [kPa])
150.0
! Parameter 2 (Pressure [kPa] or Quality [0-1] or Enthalpy [kJ/kg] or Entropy [kJ/kg-K])
101.325