real(8) function entropy(TK,p,qtot)
!-----------------------------------------------------------------------------------------
!
! TK(K),p(mb),qtot(kg/kg)
! from Raymond and Blyth 1992
!
     integer , parameter :: r8     = selected_real_kind(12)
     real(r8), parameter :: cpliq  = 4.188e3_r8    ! specific heat of fresh h2o ~ J/kg 
     real(r8), parameter :: cpres  = 1.00464e3_r8  ! specific heat of dry air ~ J/kg/K
     real(r8), parameter :: cpwv   = 1.810e3_r8    ! specific heat of water vap ~ J/kg
     real(r8), parameter :: eps1   = 18.016_r8/28.966_r8                     ! ratio of h2o to dry air molecular weights
     real(r8), parameter :: rgas   = 6.02214e26_r8*1.38065e-23_r8/28.966_r8  ! Dry air gas constant ~ J/K/kg
     real(r8), parameter :: rh2o   = 6.02214e26_r8*1.38065e-23_r8/18.016_r8  ! Water vapor gas constant ~ J/K/kg
     real(r8), parameter :: rl     = 2.501e6_r8    ! latent heat of evaporation ~ J/kg
     real(r8), parameter :: tfreez = 273.15_r8     ! freezing T of fresh water ~ K
!
     real(r8), intent(in) :: p,qtot,TK
     real(r8) :: qv,qst,e,est,L,eref,pref

pref = 1000.0_r8           ! mb
eref = 6.106_r8            ! sat p at tfreez (mb)

L = rl - (cpliq - cpwv)*(TK-tfreez)         ! T IN CENTIGRADE

! Replace call to satmixutils.

call qmmr_hPa(TK, p, est, qst)

qv = min(qtot,qst)                         ! Partition qtot into vapor part only.
e = qv*p / (eps1 +qv)

entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + &
        L*qv/TK - qv*rh2o*log(qv/qst)
! 
return
end FUNCTION entropy
