SUBROUTINE ientropy (rcall,icol,lchnk,s,p,qt,T,qst,Tfg)
!-----------------------------------------------------------------------------------------
!
! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg). 
! Inverts entropy, pressure and total water qt 
! for T and saturated vapor mixing ratio
! 

!     use phys_grid, only: get_rlon_p, get_rlat_p
     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
!
     integer, intent(in) :: icol, lchnk, rcall
     real(r8), intent(in)  :: s, p, Tfg, qt
     real(r8), intent(out) :: qst, T
     real(r8) :: qv,Ts,dTs,fs1,fs2,est
     real(r8) :: pref,eref,L,e
     !!real(r8) :: this_lat,this_lon
     integer :: LOOPMAX,i

LOOPMAX = 100                   !* max number of iteration loops 

! Values for entropy
pref = 1000.0_r8           ! mb ref pressure.
eref = 6.106_r8           ! sat p at tfreez (mb)

! Invert the entropy equation -- use Newton's method

Ts = Tfg                  ! Better first guess based on Tprofile from conv.

converge: do i=0, LOOPMAX

   L = rl - (cpliq - cpwv)*(Ts-tfreez) 

   call qmmr_hPa(Ts, p, est, qst)
   qv = min(qt,qst) 
   e = qv*p / (eps1 +qv)  ! Bolton (eq. 16)
   fs1 = (cpres + qt*cpliq)*log( Ts/tfreez ) - rgas*log( (p-e)/pref ) + &
        L*qv/Ts - qv*rh2o*log(qv/qst) - s
   
   L = rl - (cpliq - cpwv)*(Ts-1._r8-tfreez)         

   call qmmr_hPa(Ts-1._r8, p, est, qst)
   qv = min(qt,qst) 
   e = qv*p / (eps1 +qv)
   fs2 = (cpres + qt*cpliq)*log( (Ts-1._r8)/tfreez ) - rgas*log( (p-e)/pref ) + &
        L*qv/(Ts-1._r8) - qv*rh2o*log(qv/qst) - s 
   
   dTs = fs1/(fs2 - fs1)
   Ts  = Ts+dTs
   if (abs(dTs).lt.0.001_r8) exit converge
   if (i .eq. LOOPMAX - 1) then
      !this_lat = get_rlat_p(lchnk, icol)*57.296_r8
      !this_lon = get_rlon_p(lchnk, icol)*57.296_r8
      !write(iulog,*) '*** ZM_CONV: IENTROPY: Failed and about to exit, info follows ****'
      !write(iulog,100) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, &
      ! ' lat: ',this_lat,' lon: ',this_lon, &
      ! ' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._r8*qt, &
      ! ' qst(g/kg) = ', 1000._r8*qst,', s(J/kg) = ',s
      !call endrun('**** ZM_CONV IENTROPY: Tmix did not converge ****')
      print *, '**** ZM_CONV IENTROPY: Tmix did not converge ****'
      print *,"lchnk=",lchnk
      print *,"rcall=",rcall
      stop
   end if
end do converge

! Replace call to satmixutils.

call qmmr_hPa(Ts, p, est, qst)

qv = min(qt,qst)                             !       /* check for saturation */
T = Ts 

 100    format (A,I1,I4,I4,7(A,F6.2))

return
end SUBROUTINE ientropy
