C WVFD85.f

C VIS-5D analysis function to compute the Water Vapor Flux Divergence
C averaged over 1000-850 hPa (in 1.E-3 g m-2 s-1)

C See example.f for more information about writing functions.


      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, ir, ic, il,iu,iv,i,ipres,k,j
      integer it,ihum,imap
      real p(22),grav,fw,rd,cp,pi,rt,a1,a2,a4,ds
      real xincx,xincy,divk,divk2,escal,pres

c VARIABLES CALCULADAS SOBRE LOS NIVELES DE PRESION DE REFERENCIA
	REAL Q(200,200,22)

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

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

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

C   If T, RH or U, V or MAPFACCR not found, return error 1
      if (it .eq. -1 .or. ihum .eq. -1.or.iu.eq.-1.or.iv.eq.-1
     $    .or. imap .eq. -1) then
         print *, "Couldn't find T,RH, or U, V, MAPFACCR vars!"
         userfunc = 1
         return
      endif
      
      if (projection.Eq.2) then
        ds=proj_args(6)*1000.                     
      else if (projection.Eq.3) then  
        ds=proj_args(5)*1000. 
      else
        print *, "Uncomprensible projection data!"
        userfunc = 1
        return
      endif   

C ---------------------------------------------------------------
c Calcula los niveles de presion (hPA)
	do i=1,maxnl
	p(i)=1012.5*exp(vert_args(i)/-7.2)
	if (p(i).gt.(pres-1.).and.p(i).lt.(pres+1.)) ipres=i
	enddo

C ---------------------------------------------------------------
C Calcula la humedad especifica en g/Kg
		DO 1101 K=1,ipres
		DO 677 J=1,nc
		DO 677 I=1,nr

	a1=ingrid(i,j,k,it)
	a2=ingrid(i,j,k,ihum)
	call procio(a1,a2,a4)
	call humesp(p(k),a4,q(i,j,k))
	Q(i,j,K)=Q(i,j,K)*1000.
677 	continue
1101	CONTINUE

C Calculo de divergencia y integral
C ----------------------------------------------------------------

c Calculo de los diferenciales en x
	xincy=2.*ds
		DO 70 I=2,nr-1
	xincx=2.*ds
		DO 70 J=2,nc-1

C	------------------------------------------------------
C	CALCULO DIVERGENCIA DEL FLUJO DE VAPOR DE AGUA
	fw=0.
	do 345 k=2,ipres
	divk=(ingrid(i,j+1,k,iu)*q(i,j+1,k)-ingrid(i,j-1,k,iu)*
     & q(i,j-1,k))/xincx+(ingrid(i-1,j,k,iv)*q(i-1,j,k)-ingrid
     & (i+1,j,k,iv)*q(i+1,j,k))/xincy
        divk=divk*ingrid(i,j,1,imap)
	divk2=(ingrid(i,j+1,k-1,iu)*q(i,j+1,k-1)-ingrid(i,j-1,k-1,iu)*
     & q(i,j-1,k-1))/xincx+(ingrid(i-1,j,k-1,iv)*q(i-1,j,k-1)-ingrid
     & (i+1,j,k-1,iv)*q(i+1,j,k-1))/xincy
        divk2=divk2*ingrid(i,j,1,imap)
C Calculo de la integral usando el punto medio y pasando el dP a Pa
	fw=fw+((divk+divk2)/2.)*(p(k-1)-p(k))*100.
345	enddo
C Division final por la gravedad (Unidades finales g m-2 s-1)
	fw=fw/grav

	do il=1,outnl
		outgrid(i,j,il)=fw*escal
	enddo

70 	continue

	do il=1,outnl
	do ir=2,nr-1
	   outgrid(ir,1,il)=2.*outgrid(ir,2,il)-outgrid(ir,3,il)
	   outgrid(ir,nc,il)=2.*outgrid(ir,nc-1,il)-outgrid(ir,nc-2,il)
	enddo
	do ic=2,nc-1
	  outgrid(1,ic,il)=2.*outgrid(2,ic,il)-outgrid(3,ic,il)
	  outgrid(nr,ic,il)=2.*outgrid(nr-1,ic,il)-outgrid(nr-2,ic,il)
	enddo
	outgrid(1,1,il)=2.*outgrid(2,2,il)-outgrid(3,3,il)
	outgrid(1,nc,il)=2.*outgrid(2,nc-1,il)-outgrid(3,nc-2,il)
       outgrid(nr,1,il)=2.*outgrid(nr-1,2,il)-outgrid(nr-2,3,il)
       outgrid(nr,nc,il)=2.*outgrid(nr-1,nc-1,il)-outgrid(nr-2,nc-2,il)
	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
