undef ("relhum_water")
function relhum_water(tk:numeric,qw:numeric,p:numeric)
; ---
; computes rh (relative humidity) from p, t and qw (mixing ratio)
; ---
;  Calculate rel hum with respect to water 
;  Same constants as mixhum_ptrh reversible calculation
; ---
; NCL version of f77 subroutines donated by: Alan Brammer, U. ALbany
; ---
; definition of mixing ratio
;    es  - is the saturation mixing ratio
;    ep  - is the ratio of the molecular weights of water vapor to dry air
;    p   - is the atmospheric pressure
;    rh  - is the relative humidity (given as a percent)

;    rh =  100*  q / ( (ep*es)/(p-es) )

; input-
;    P      - pressure (Pa)  ... this will be converted in the code to hPa (mb)
;             p      - pressure (hPa or mb)
;    tk     - temperature (k)
;    qw     - mixing ratio (kg/kg)

; output-
;    rh     - relative humidity as %

local t0,ep,onemep,es0w,aw,bw,es0i,ai,bi,est,qst,pa2mb
begin
; local: constants for water and ice

      t0     = 273.15d
      ep     = 0.622d
      onemep = 1-ep
      pa2mb  = 0.01d
          
; local: constants for water (W)

      es0w   = 6.11d
      aw     = 17.269d
      bw     = 35.86d
      

      est    = es0w*exp((aw*(tk-t0)) / (tk-bw)) 
      
      qst = (ep*est)/ ((p*pa2mb)-onemep*est)

      if (all(qst.gt.0)) then
          rh  =  100* (qw / qst)
      else                                           ; ? numerical issues ?
          if (isatt(tk,"_FillValue")) then
              qst@_FillValue = tk@_FillValue 
          else
              qst@_FillValue = default_fillvalue(typeof(qst))
          end if
          rh  =  100* (qw / where(qst.le.0, qst@_FillValue, qst) )
      end if

      if (typeof(tk).eq."float" .or. typeof(qw).eq."float" \
                                .or. typeof(p) .eq."float") then
          rh := tofloat(rh)
      end if

      rh@long_name = "relative humidity"
      rh@units     = "%"
      rh@tag       = "function: relhum_water"
      copy_VarCoords(tk, rh)

      return (rh)
end