subroutine genqsat1(sph,qsat,ges_prsl,ges_tv,ice,npts,nlevs)
! this subroutine was extracted from the GSI version operational
! at NCEP in Dec. 2007. Only difference is that this version takes
! single precision 2d arrays instead of default real 3d arrays.
!
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    genqsat
!   prgmmr: derber           org: np23                date: 1998-01-14
!
! abstract: obtain saturation specific humidity for given temperature.
!
! program history log:
!   1998-01-14  derber
!   1998-04-05  weiyu yang
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   1903-10-07  Wei Gu, bug fixes,if qs<0,then set qs=0; merge w/ GSI by R Todling
!   2003-12-23  kleist, use guess pressure, adapt module framework
!   2004-05-13  kleist, documentation
!   2004-06-03  treadon, replace ggrid_g3 array with ges_* arrays
!   2005-02-23  wu, output dlnesdtv
!   2005-11-21  kleist, derber  add dmax array to decouple moisture
!               from temp and pressure for questionable qsat
!   2006-02-02  treadon - rename prsl as ges_prsl
!   2006-09-18  derber - modify to limit saturated values near top
!   2006-11-22  derber - correct bug:  es<esmax should be es<=esmax
!
!   input argument list:
!     ges_prsl  - guess pressure (mb)
!     ges_tv    - guess virtual temp (K)
!     qsat      - guess specific humidity (in), saturation value (out)
!     ice       - logical flag:  T=include ice and ice-water effects,
!                 depending on t, in qsat calcuations.
!                 otherwise, compute qsat with respect to water surface
!
!   output argument list:
!     qsat      - saturation specific humidity (output)
!
! remarks: see modules used
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_double,i_kind,r_single,r_kind
  use constants, only: xai,tmix,xb,omeps,eps,xbi,one,zero,&
       xa,psat,ttp,fv
  implicit none

  logical,intent(in):: ice
  integer(i_kind), intent(in) :: npts,nlevs
  real(r_double),dimension(npts,nlevs),intent(out):: qsat
  real(r_single),intent(in),dimension(npts,nlevs) :: ges_prsl,ges_tv,sph

  integer(i_kind) k,i
  real(r_double) pw,tdry,tr,es
  real(r_double) w,onep2,esmax
  real(r_double),dimension(npts):: mint,estmax
  real(r_double),dimension(npts,nlevs):: ges_tsen
  integer(i_kind),dimension(npts):: lmint

  onep2 = 1.e2_r_double

  ! compute sensible from virtual temp.
  ges_tsen=ges_tv/(one+fv*sph)
  !print &
  !*,'ges_tsen',minval(ges_tsen),maxval(ges_tsen),minval(ges_prsl),maxval(ges_prsl),ice,psat

  lmint=1
  do i=1,npts
    mint(i)=340._r_double
  end do
  do k=1,nlevs
      do i=1,npts
        if((ges_prsl(i,k) < 30._r_double .and.  &
            ges_prsl(i,k) > 2._r_double) .and.  &
            ges_tsen(i,k) < mint(i))then
           lmint(i)=k
           mint(i)=ges_tsen(i,k)
         end if
      end do
  end do
  do i=1,npts
      tdry = mint(i)
      tr = ttp/tdry
      if (tdry >= ttp .or. .not. ice) then
        estmax(i) = psat * (tr**xa) * exp(xb*(one-tr))
      elseif (tdry < tmix) then
        estmax(i) = psat * (tr**xai) * exp(xbi*(one-tr))
      else
        w  = (tdry - tmix) / (ttp - tmix)
        estmax(i) =  w * psat * (tr**xa) * exp(xb*(one-tr)) &
                + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr))
      endif
  end do
  if (ice) then
    do k = 1,nlevs
        do i = 1,npts

          pw = onep2*ges_prsl(i,k) ! convert to Pa from mb.
              
          tdry = ges_tsen(i,k)
          tr = ttp/tdry
          if (tdry >= ttp) then
            es = psat * (tr**xa) * exp(xb*(one-tr))
          elseif (tdry < tmix) then
            es = psat * (tr**xai) * exp(xbi*(one-tr))
          else
            w  = (tdry - tmix) / (ttp - tmix)
            es =  w * psat * (tr**xa) * exp(xb*(one-tr)) &
                  + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr))
          endif

          esmax = es
          if(lmint(i) < k)then
             esmax=0.1*pw
             esmax=min(esmax,estmax(i))
          end if

          if(es > esmax)then
            es = esmax
          end if
          qsat(i,k) = eps * es / (pw - omeps * es)
          qsat(i,k) = max(tiny(qsat(i,k)),qsat(i,k))

        end do
    end do

! Compute saturation values with respect to water surface
  else
    do k = 1,nlevs
        do i = 1,npts
              
          pw = onep2*ges_prsl(i,k)  ! convert to Pa from mb

          tdry = ges_tsen(i,k)
          tr = ttp/tdry
          es = psat * (tr**xa) * exp(xb*(one-tr))
          esmax = es
          if(lmint(i) < k)then
             esmax=0.1*pw
             esmax=min(esmax,estmax(i))
          end if
          if(es > esmax)then
            es = esmax
          end if
          qsat(i,k) = eps * es / (pw - omeps * es)
          qsat(i,k) = max(tiny(qsat(i,k)),qsat(i,k))

        end do
    end do
     
  endif   ! end if ice

  ! enforce min value.
  where (qsat < 1.e-6) qsat = 1.e-6
  
  return
end subroutine genqsat1