!
      module m_relhum
!
! Module for transforming between specific and relative humidities
!
! For the saturation vapor pressure v this uses a linear interpolation
! between values from a look-up table constructed in the routine 
! relhum_setup.   These are starurations with repect to liquid water. 
! Then, given ta=t-tmin it 
! determines the largest value of t1<ta such that t1 is an integer 
! multiple of qsatvp(3).  tx-ip is the weight w=(ta-t1)/qsatvp(3)
!         qs=ep2*v/(p-v)
!
! Initial code by Ronald Errico December 2007
!
      use m_kinds
!
      implicit none
!
      private
      public :: relhum_setup
      public :: relhum_rh2q
      public :: relhum_q2rh
      public :: relhum_td2q
!
      real(rkind1), parameter :: ep2=0.622
      integer, parameter :: nsatvp=10000  ! size of look-up table for sat v.p
      real(rkind1), save :: satvp(nsatvp) ! table of e(T)=saturation vapor p.
!
      contains
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine relhum_setup
!
!  Construct tables for values for saturation vapor pressure for 
!  selected values of T.  Values of T in units of (tmax-tmin)/(nsatvp-4). 
!  For the interpolation algorithm, see routines qsat.
! local 
      implicit none
      integer :: n
      real(rkind1), parameter :: svpt0=273.15, svp1=611.2
      real(rkind1), parameter :: svp2=17.67, svp3=29.65
      real(rkind1), parameter :: tmin=100., tmax=340.
      real(rkind1) :: tr,t,a
!     
      satvp(1)=tmax
      satvp(2)=tmin
      tr=(tmax-tmin)/(nsatvp-4)
      satvp(3)=tr
      do n=4,nsatvp 
        t=tmin+(n-4)*tr
        a=svp2*(t-svpt0)/(t-svp3)
        satvp(n)=svp1*exp(a)
      enddo 
!
        print *,' '
        print *,' Table for satvp filled in module m_relhum setup'
!
      end subroutine relhum_setup
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine relhum_q2rh (t,p,q,ier)
!
! Replace specific humidity by relative humidity
!
      implicit none
      integer      :: ier
      real(rkind2) :: t, p, q
! local
      integer      :: ip
      real(rkind1) :: tx, v
!
      tx=(t-satvp(2))/satvp(3)  ! t-tmin in units of satvp(2)
      ip=int(tx)
      v=satvp(ip+4)+(tx-ip)*(satvp(ip+5)-satvp(ip+4))
      if (p > v) then
        q=q*(p-v)/(v*ep2)
      else
        ier=ier+1
      endif
!
      end subroutine relhum_q2rh 
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine relhum_rh2q (t,p,q,ier)
!
! Replace specific humidity by relative humidity
!
      implicit none
      integer      :: ier
      real(rkind1) :: t, p, q
! local
      integer      :: ip
      real(rkind1) :: tx, v
!
      tx=(t-satvp(2))/satvp(3)  ! t-tmin in units of satvp(2)
      ip=int(tx)
      v=satvp(ip+4)+(tx-ip)*(satvp(ip+5)-satvp(ip+4))
      if (p > v) then
        q=ep2*q*v/(p-v)
      else
        ier=ier+1
      endif
!
      end subroutine relhum_rh2q 
!
! 
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine relhum_td2q (tq,p)
!
! Replace dew point temperature (input in array tq) by specific humidity 
! (output array tq).
!
      implicit none
      integer      :: ier
      real(rkind2) :: tq, p
! local
      integer      :: ip
      real(rkind1) :: tx, v, w
!
      tx=(tq-satvp(2))/satvp(3)  ! t-tmin in units of satvp(2)
      ip=int(tx)
      v=satvp(ip+4)+(tx-ip)*(satvp(ip+5)-satvp(ip+4))
      w=ep2*v/(p-v)
      tq=w/(1.+w)
!
      end subroutine relhum_td2q 
!
!    
      end module m_relhum