C QV.f

C VIS-5D analysis function to compute the water vapor mixing ratio (in g/kg)

      INTEGER FUNCTION USERFUNC( OUTGRID, OUTNL, OUTLOWLEV,
     *                           INGRID, NR, NC, NL, LOWLEV, MAXNL,
     *                           NVARS, NAMES,
     *                           DATE, TIME,
     *                           PROJECTION, PROJ_ARGS,
     *                           VERTICAL, VERT_ARGS )
	implicit none
C     ARGUMENTS:
      INTEGER NVARS
      INTEGER NR, NC, NL(NVARS), LOWLEV(NVARS), MAXNL
      REAL OUTGRID(NR,NC,MAXNL)
      INTEGER OUTNL, OUTLOWLEV
      REAL INGRID(NR,NC,MAXNL,NVARS)
      CHARACTER*8 NAMES(NVARS)
      INTEGER DATE, TIME
      INTEGER PROJECTION
      REAL PROJ_ARGS(*)
      INTEGER VERTICAL
      REAL VERT_ARGS(*)

C     LOCAL VARS:
      integer var, il,i,j
      integer it,ihum
      real p,grav,rd,cp,pi,rt,a1,a2,a4
      real q,pres,escal

c VARIABLES CALCULADAS SOBRE LOS NIVELES DE PRESION DE REFERENCIA
	DIMENSION Q(100,100,22)

      PARAMETER(rd=287.04,cp=1004.64,Pi=3.141596,Rt=6371227.,grav=9.8
     & ,escal=1.e3,pres=700.)

C   Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0

C   Find the T and RH variables
      it = -1
      ihum = -1
      do var=1,nvars
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'RH') ihum = var
      end do

C   If T and RH not found, return error 1
      if (it .eq. -1 .or. ihum .eq. -1) then
         print *, "Couldn't find T, RH variables!"
         userfunc = 1
         return
      endif
C ---------------------------------------------------------------
c Calcula los niveles de presion (hPA)
	do il=1,maxnl
	p=1012.5*exp(vert_args(il)/-7.2)
C ---------------------------------------------------------------
C Calcula la humedad especifica en g/Kg
		DO 677 J=1,nc
		DO 677 I=1,nr
	a1=ingrid(i,j,il,it)
	a2=ingrid(i,j,il,ihum)
	call procio(a1,a2,a4)
	call humesp(p,a4,q(i,j,il))
	outgrid(i,j,il)=Q(i,j,il)*1000.

677 	continue
	enddo
      USERFUNC = 0
      RETURN
      END

c -----------------------------------------------------------------
	subroutine procio(t,u,e)

	es=6.112*EXP((17.67*t)/(t+243.5))
	e=es*u/100.
	return
	end

c ------------------------------------------------------------------
	subroutine humesp (p,e,qq)

	r=.622*e/(p-e)
	qq=r/(1+r)
	return
	end
