Index: doninip.ftn
===================================================================
--- doninip.ftn	2004-11-15 17:15:45.000000000 +0000
+++ doninip.ftn	2004-11-15 15:58:00.000000000 +0000
@@ -62,6 +62,12 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
          call rsdata ( refip1(1,2),lnk,unf,'ES  ',datestp,ip3,nv(3,1),
      $                 donut )
       endif
+      if ((nv(3,1).eq.'!@#$%^&*').and.(gngalsig.eq.0)
+     $                           .and.(nv(2,1)(1:2).eq.'TT')) then
+	 nv(3,2) = 'HR'
+	 call rsdata ( refip1(1,2),lnk,unf,'HR  ',datestp,ip3,nv(3,1),
+     $		       donut )
+      endif
 *
       inv = 4
       if (gngalsig.ne.1) then
Index: rsdata.ftn
===================================================================
--- rsdata.ftn	2004-11-15 17:15:45.000000000 +0000
+++ rsdata.ftn	2004-11-15 16:22:09.000000000 +0000
@@ -18,6 +18,7 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
       integer ip2,i,ii,j,jj,k,err_rdh,cnt,s_rdhint,iproc,i0,in,j0,jn,
      $        donutx
       external s_rdhint
+      logical haveVariable
       real wk1(Grd_ni,Grd_nj),wk2(Grd_ni*Grd_nj)
 *----------------------------------------------------------------------
 *  
@@ -25,13 +26,24 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
       typvar = ' '
       ip2    = -1
       donutx = -1
+      haveVariable = .false.
       if (donut.gt.0) donutx = donut + 1
 *      
       do 110 k=1,lnk
 *
          err_rdh = s_rdhint (wk1,xpx,ypx,Grd_ni,Grd_nj,nomvar,idatev,
      $          levels(k),ip2,ip3,etiket,typvar,.false.,hint_ntr,iun,6)
-         if (err_rdh.lt.0) goto 500
+*  Allow for missing humidity fields above the lowest level (set to 0.)
+	 if (err_rdh < 0 .and. haveVariable .and. (
+     $        trim(nomvar) .eq. 'HR' .or.
+     $        trim(nomvar) .eq. 'HU' .or.
+     $        trim(nomvar) .eq. 'ES')) then
+	    write(6,1000) levels(k)
+	    wk1 = 0.
+         else if (err_rdh.lt.0) then
+	    goto 500
+	 endif
+	 haveVariable = .true.
 *
          if (donut.le.0) then
 *
@@ -105,6 +117,7 @@ c     $                          1,1,1,l
 *
       nv = nomvar
 *
+ 1000 format(' (RSDATA) NO DATA FOR HUMIDITY FIELD FOR ',i4,': SET TO 0.')
 *----------------------------------------------------------------------
  500  return
       end
Index: s_rdvint.ftn
===================================================================
--- s_rdvint.ftn	2004-11-15 17:15:45.000000000 +0000
+++ s_rdvint.ftn	2004-11-15 17:14:09.000000000 +0000
@@ -37,7 +37,7 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
       external bmf_get,get_bmf
       integer i,j,k,n,ng,lka,err,i0,in,j0,jn,nest,
      $        l_i0,l_in,l_j0,l_jn,suv,km1,kp1
-      real td,trpl,eps1,eps2,delta,knams
+      real td,trpl,eps1,eps2,delta,knams,es
       real , dimension (:  ), allocatable :: p_ref
       real , dimension (:,:), allocatable :: topo_anal,gotsr,ortsr,
      $     qtsr,ntsr2,hgeow,hgeot,hgeom,t_anal,hgeow_anal,hgeot_anal,
@@ -167,18 +167,45 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
       err = get_bmf ('HU  ',0,0,q_anal,i0,in,j0,jn,nka,nia,nja)
       if (err.ne.0) then
          allocate (p_ref(nka))
-         err = get_bmf ('ES  ',0,0, q_anal,i0,in,j0,jn,nka,nia,nja)
-         if (err.ne.0) call mc2stop(-1)
-         err = get_bmf ('TT  ',0,0, t_anal,i0,in,j0,jn,nka,nia,nja)
+	 err = get_bmf ('TT  ',0,0, t_anal,i0,in,j0,jn,nka,nia,nja)
          if (err.ne.0) call mc2stop(-1)
          err = get_bmf ('PREF',0,0, p_ref ,1,1,1,1,nka,1,1)
          if (err.ne.0) call mc2stop(-1)
-         do k=1,nka
-         do i=1,ng
-            td = (t_anal(i,k) - q_anal(i,k)) + tcdk_8
-            q_anal(i,k)  = foqst(td,p_ref(k)*100.)
-         end do
-         end do
+         err = get_bmf ('ES  ',0,0, q_anal,i0,in,j0,jn,nka,nia,nja)  ! Try dewpoint depression
+	 if (err.ne.0) then  ! No? ... try relative humidity
+	    err = get_bmf ('HR  ',0.,0., q_anal,i0,in,j0,jn,nka,nia,nja) ! Now convert HR
+	    if (err.ne.0) call mc2stop(-1)
+	    if (maxval(q_anal) > 10.) q_anal = q_anal/100.  ! convert to frac if required
+*
+***  Use equations from Gill (1982) to convert relative humidity (HR) to specific
+***  humidity (HU).  Note that these equations are only strictly valid from -40C to 40C.
+*
+	    do k=1,nka
+	    do i=1,ng
+	       es = min(q_anal(i,k),1.)*10.**((.7859+.03477*(t_anal(i,k))) / 
+     +            (1.+.00412*(t_anal(i,k)))+2.)
+	       q_anal(i,k) = max(.622*es / (p_ref(k)*100.-0.378*es), 0.)
+	    enddo
+	    enddo
+*
+***  Use equations from Rogers and Yau to convert relative humidity (HR) to specific
+***  humidity (HU).  These equations are technically valid over a larger range.
+*
+c	    do k=1,nka
+c	    do i=1,ng	
+c	       es = (3.8/p_ref(k))*exp(17.27*(t_anal(i,k) / 
+c     +            (t_anal(i,k)+tcdk_8-35.86)))
+c	       q_anal(i,k) = max(es * min(q_anal(i,k),1.), 0.)
+c	    enddo
+c	    enddo
+	 else  ! Yes ... now convert dewpoint depression
+            do k=1,nka
+            do i=1,ng
+               td = (t_anal(i,k) - q_anal(i,k)) + tcdk_8
+               q_anal(i,k)  = foqst(td,p_ref(k)*100.)
+            end do
+            end do
+         endif
          deallocate (p_ref)
       endif
       call vertint3 (qv,q_anal,posit(ng*gnk*2+1),htt,ng,gnk,nka)