*************************************************************************

function qvtd(P,QV)

*-----------------------------------------------------------------------
*
*  Calculates the dewpoint temperature (Td, in K) given the air
*  pressure (Pa) and water vapor mixing ratio (kg[H2O]/kg[dry air])
*
*  Solve Teten's formula
*
*         QVsat = (380./P)*exp(17.2693882*(Td-Tfr)/(Td-35.86))
*
*  for the dewpoint temperature, Td:
*
*         Td = (Tfr - (c/b)*ln(qvsat/a))/(1 - ln(qvsat/a)/b)
*
*     where
*
*         Tfr = freezing temperature (K)
*         P = air pressure (Pa)
*         QVsat = actual water vapor mixing ratio (kg[H2O]/kg[dry air])
*                 (= saturation mixing ratio at Td)
*
*         a = 380./P
*         b = 17.2693882
*         c = 35.86
*
*  Compare to Equation 13.1.4a, page 551 in
*  "An Introduction to Boundary Layer Meteorology" by Roland B. Stull (1988).
*
*  Note, uses function "math_log(x)" to calculate the natural logarithm
*  from:  http://grads.iges.org/grads/gadoc/mathfunctions.html
*
*-----------------------------------------------------------------------
*
say 'Function "qvtd":'
*
* constants required for function
*
Tfr = 273.15
b = 17.2693882
c = 35.86
a = 380./P

say 'P(hPa) = 'P*0.01
say 'QV(kg/kg) = 'QV
say 'QV(g/kg)  = 'QV*1000.

qvsat = QV
qvsat_ratio = qvsat/a
if( qvsat_ratio > 0. )
  Td = (Tfr-(c/b)*math_log(qvsat/a))/(1-math_log(qvsat/a)/b)
else
  say
  say '"qvsat/a" is not greater than zero;  qvsat/a = 'qvsat/a
  say 'exit script...'
  say
  exit
endif

say 'Td(C) = 'Td-Tfr

return(Td)

*************************************************************************
