source: trunk/LMDZ.GENERIC/libf/phystd/phasechange_kcm.F90 @ 848

!  ----------------------------------------------------------------
!  Purpose: Thermodynamic data on H2O, NH3
!  Authour: Adapted from various sources by R. Wordsworth (2011)


! ---------------------
! NH3
! ---------------------
subroutine psat_NH3 ( T, p )
!
!  PSAT_NH3 makes a rough estimate of the vapor pressure.
!
!  Interpolated from www.engineeringtoolbox.com data by RDW 21/09/11 for
!  temperatures between 223 and 323 K.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Output, double precision P, the vapor pressure, in MegaPascals.
!

  implicit none

  double precision p,T

  p=exp(-1.5609d-004*T**2 + 0.1236*T - 9.1530)/1d6

end subroutine psat_NH3

subroutine latheat_NH3 ( T, sc, sv )
!
!  PSAT_NH3 makes a rough estimate of the entropies of condensation /
!  vapourisation.
! 
!  Interpolated from www.engineeringtoolbox.com data by RDW 21/09/11 for
!  temperatures between 223 and 323 K.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Output, double precision sc, the entropy of condensate, in J kg^-1 K^-1.
!            double precision sv, the entropy of gas, in J kg^-1 K^-1.
!

  implicit none

  double precision T,sc,sv

  sv = 0.0492*T**2 - 40.4199*T + 1.2708e+004
  sc = -0.0215*T**2 + 28.7138*T - 5.5267e+003

  !L = -11.2373*T**2 + 2.5326d+03*T + 1.4099d+06

end subroutine latheat_NH3


! ---------------------
! H2O
! ---------------------

subroutine base ( t, rho, ab, cvb, dpdrb, dpdtb, gb, hb, pb, sb, ub )
!
!*******************************************************************************
!
!! 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
!
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Modified:
!
!    03 February 2002
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Input, double precision RHO, the density, in G/CM3.
!
!    Output, double precision AB, the base value of the Helmholtz function,
!    in KJ/kg.
!
!    Output, double precision CVB, the base value of the isochoric (constant 
!    volume) heat capacity, in KJ/(kg degrees Kelvin).
!
!    Output, double precision DPDRB, the base value of the partial
!    derivative dP(T,RHO)/dRHO, with T held fixed, in (MegaPascals CM3)/G.
!
!    Output, double precision DPDTB, the base value of the partial
!    derivative dP(T,RHO)/dT, with RHO held fixed, in 
!    MegaPascals/degrees Kelvin.
!
!    Output, double precision GB, the base value of the Gibbs free energy,
!    in KJ/kg.
!
!    Output, double precision HB, the base value of enthalpy, in KJ/kg.
!
!    Output, double precision PB, the base pressure, in MegaPascals.
!
!    Output, double precision SB, the base value of entropy, 
!    in KJ/(kg degrees Kelvin).
!
!    Output, double precision UB, the base value of internal energy, 
!    in KJ/kg.
!
  implicit none
!
  double precision ab
  double precision, parameter :: alpha = 11.0D+00
  double precision b1
  double precision b1t
  double precision b1tt
  double precision b2
  double precision b2t
  double precision b2tt
  double precision, parameter :: beta = 44.333333333333D+00
  double precision cvb
  double precision dpdrb
  double precision dpdtb
  double precision dz
  double precision dz0
  double precision, parameter :: gamma = 3.5D+00
  double precision gascon
  double precision gb
  double precision hb
  double precision, parameter :: p_zero = 0.101325D+00
  double precision pb
  double precision rho
  double precision sb
  double precision t
  double precision ub
  double precision x
  double precision y
  double precision z
  double precision z0
!
!  Refuse to handle zero or negative temperatures.
!
  if ( t <= 0.0 ) 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 base

subroutine bb ( t, b1, b2, b1t, b2t, b1tt, b2tt )
!
!*******************************************************************************
!
!! 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
!
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Output, double precision B1, the coefficient b from equation 3, 
!    in CM3/G.
!
!    Output, double precision B2, the coefficient Bbar from equation 4, 
!    in CM3/G.
!
!    Output, double precision B1T, the derivative dB1/dT, 
!    in (CM3)/(G Degrees Kelvin).
!
!    Output, double precision B2T, the derivative dB2/dT, 
!    in (CM3)/(G Degrees Kelvin).
!
!    Output, double precision B1TT, the second derivative of B1 with 
!    respect to T, in (CM3)/(G (Degrees Kelvin)**2 ).
!
!    Output, double precision B2TT, the second derivative of B2 with 
!    respect to T, in (CM3)/(G (Degrees Kelvin)**2 ).
!
  implicit none
!
  double precision b1
  double precision b1t
  double precision b1tt
  double precision b2
  double precision b2t
  double precision b2tt
  double precision, 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 /)
  double precision, 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
  double precision t
  double precision, parameter :: t_ref = 647.073D+00
  double precision 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 bb


subroutine ideal ( t, ai, cpi, cvi, gi, hi, si, ui )
!
!*******************************************************************************
!
!! 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.
!     
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    TJ270.H3
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Output, double precision AI, the Helmholtz function, in KJ/kg.
!
!    Output, double precision CPI, the heat capacity at constant pressure,
!    in KJ/(kg degrees Kelvin).
!
!    Output, double precision CVI, the heat capacity at constant volume,
!    in KJ/(kg degrees Kelvin).
!
!    Output, double precision GI, the Gibbs free energy, in KJ/kg.
!
!    Output, double precision HI, the enthalpy, in KJ/kg.
!
!    Output, double precision SI, the entropy, in KJ/(kg degrees Kelvin).
!
!    Output, double precision UI, the internal energy, in KJ/kg.
!
  implicit none
!
  double precision ai
  double precision, 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 /)
  double precision cpi
  double precision cvi
  double precision gascon
  double precision gi
  double precision hi
  integer i
  double precision si
  double precision t
  double precision temp
  double precision tt
  double precision 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 ideal


subroutine resid ( t, rho, ar, cvr, dpdrr, dpdtr, gr, hr, pr, sr, ur )
!
!*******************************************************************************
!
!! 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, P > 300 MPa.
!
!    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.
!   
!  Modified:
!
!    22 November 1998
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Input, double precision RHO, the density, in G/CM3.
!
!    Output, double precision AR, the residual contribution to the 
!    Helmholtz function, in KJ/kg.
!
!    Output, double precision CVR, the residual contribution to the 
!    isochoric (constant volume) heat capacity, in KJ/(kg degrees Kelvin).
!
!    Output, double precision DPDRR, the residual contribution to 
!    the partial derivative dP(T,RHO)/dRHO, with T held fixed, in 
!    (MegaPascals CM3)/G.
!
!    Output, double precision DPDTR, the residual contribution to 
!    the partial derivative dP(T,RHO)/dT, with RHO held fixed, 
!    in MegaPascals/degrees Kelvin.
!
!    Output, double precision GR, the residual contribution to the Gibbs 
!    function, in KJ/kg.
!
!    Output, double precision HR, the residual contribution to the 
!    enthalpy, in KJ/kg.
!
!    Output, double precision PR, the residual contribution to the pressure, 
!    in MegaPascals.
!
!    Output, double precision SR, the residual contribution to the entropy, 
!    in KJ/(kg degrees Kelvin).
!
!    Output, double precision UR, the residual contribution to the 
!    internal energy, in KJ/kg.
!
  implicit none
!
  double precision, parameter, dimension ( 4 ) :: aad = (/ &
    34.0D+00, 40.0D+00, 30.0D+00, 1050.0D+00 /)
  double precision, parameter, dimension ( 4 ) :: aat = (/ &
    20000.0D+00, 20000.0D+00, 40000.0D+00, 25.0D+00 /)
  double precision, parameter, dimension ( 4 ) :: adz = (/ &
    0.319D+00, 0.319D+00, 0.319D+00, 1.55D+00 /)
  double precision ar
  double precision att
  double precision, parameter, dimension ( 4 ) ::  atz = (/ &
    640.0D+00, 640.0D+00, 641.6D+00, 270.0D+00 /)
  double precision cvr
  double precision dadt
  double precision ddz
  double precision del
  double precision dex
  double precision dfdt
  double precision dpdrr
  double precision dpdtr
  double precision e
  double precision errtol
  double precision ex0
  double precision ex1
  double precision ex2
  double precision fct
  double precision, 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 /)
  double precision gascon
  double precision gr
  double precision 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
  double precision pr
  double precision q10
  double precision q20
  double precision q2a
  double precision q5t
  double precision qm
  double precision qp
  double precision qr(11)
  double precision qt(10)
  double precision rho
  double precision sr
  double precision, parameter :: s_ref = 7.6180720166752D+00
  double precision t
  double precision, parameter :: t_ref = 647.073D+00
  double precision tau
  double precision tx
  double precision, parameter :: u_ref = - 4328.4549774261D+00
  double precision ur
  double precision 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**2
!   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 resid

subroutine psat_H2O ( t, p )
!
!*******************************************************************************
!
!! PSAT_H2O 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.
!
!  Modified:
!
!    21 November 1998
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Output, double precision P, the vapor pressure, in MegaPascals.
!
  implicit none
!
  double precision, 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 /)
  double precision b
  integer i
  double precision p
  double precision q
  double precision t
  double precision, parameter :: t_ref = 647.25D+00
  double precision v
  double precision w
  double precision z
!
!  Refuse to handle zero or negative temperatures.
!
  if ( t <= 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PSAT_H2O - 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_H2O


subroutine tdpsdt ( t, dp )
!
!*******************************************************************************
!
!! 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.
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Output, double precision DP, the value T*(dP(Sat)/dT), 
!    in MegaPascals.
!
  implicit none
!
  double precision, 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 /)
  double precision b
  double precision c
  double precision dp
  integer i
  double precision q
  double precision t
  double precision, parameter :: t_ref = 647.25D+00
  double precision v
  double precision w
  double precision y
  double precision 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 tdpsdt



subroutine therm ( t, rho, a, cjth, cjtt, cp, cv, dpdr, dpdt, g, h, p, s, u )
!
!*******************************************************************************
!
!! 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
!   
!  Modified:
!
!    19 November 1998
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Parameters:
!
!    Input, double precision T, the temperature, in degrees Kelvin.
!
!    Input, double precision RHO, the fluid density, in G/CM3.
!
!    Output, double precision A, the Helmholtz function, in KJ/kg.
!
!    Output, double precision CJTH, the Joule-Thomson coefficient,
!    in K/MegaPascals.
!
!    Output, double precision CJTT, the isothermal Joule-Thomson coefficient,
!    in CM3/G.
!
!    Output, double precision CP, the isobaric (constant pressure) heat
!    capacity, in KJ/(kg degrees Kelvin).
!
!    Output, double precision CV, the isochoric (constant volume) heat capacity,
!    in KJ/(kg degrees Kelvin).
!
!    Output, double precision DPDR, the partial derivative 
!    dP(T,RHO)/dRHO, with T held fixed, in MegaPascals*CM3/G.
!
!    Output, double precision DPDT, the partial derivative 
!    dP(T,RHO)/dT, with RHO held fixed, in MegaPascals/degrees Kelvin.
!
!    Output, double precision G, the Gibbs free energy, in KJ/kg.
!
!    Output, double precision H, the enthalpy, in KJ/kg.
!
!    Output, double precision P, the pressure, in MegaPascals.
!
!    Output, double precision S, the entropy, in KJ/(kg degrees Kelvin).
!
!    Output, double precision U, the internal energy, in KJ/kg.
!
  implicit none
!
  double precision a
  double precision ab
  double precision ai
  double precision ar
  double precision cjth
  double precision cjtt
  double precision cp
  double precision cpi
  double precision cv
  double precision cvb
  double precision cvi
  double precision cvr
  logical, parameter :: debug = .false.
!  logical, parameter :: debug = .true.
  double precision dpdr
  double precision dpdrb
  double precision dpdrr
  double precision dpdt
  double precision dpdtb
  double precision dpdtr
  double precision g
  double precision gb
  double precision gi
  double precision gr
  double precision h
  double precision hb
  double precision hi
  double precision hr
  double precision p
  double precision pb
  double precision pr
  double precision rho
  double precision s
  double precision sb
  double precision si
  double precision sr
  double precision t
  double precision u
  double precision ub
  double precision ui
  double precision 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 therm

function gascon ( )
!
!*******************************************************************************
!
!! GASCON returns the value of the specific gas constant.
!
!
!  Note:
!
!    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.
!
!  Reference:
!
!    Lester Haar, John Gallagher and 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,
!    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, 1967.
!
!  Parameters:
!
!    Output, double precision GASCON, the value of the specific gas 
!    constant, in J/(g degrees Kelvin).
!
  implicit none
!
  double precision gascon
!
  gascon = 0.461522D+00
 
  return
end function gascon