C UQ.f
C
C VIS-5D analysis function to compute U Component of Q Vector
C (in  1.E-10 m2 kg-1 s-1)

      INTEGER FUNCTION USERFUNC( OUTGRID, OUTNL, OUTLOWLEV,
     *                           INGRID, NR, NC, NL, LOWLEV, MAXNL,
     *                           NVARS, NAMES,
     *                           DATE, TIME,
     *                           PROJECTION, PROJ_ARGS,
     *                           VERTICAL, VERT_ARGS )
C     ARGUMENTS:
	implicit none
      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,it,ig,i
      integer ilat,imap
      real lat,Pi,Rt,Per
      real q1
      real p(22),r
      real g,cor,xincx,xincy,b20,b21,ds
      parameter (Pi=3.141596,Rt=6371227.,Per=86400.,g=9.807,r=287.04)

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

C   Find the T, H and LATITCRS and MAPFACCR variables
      it = -1
      ig = -1
      ilat = -1
      imap = -1
      do var=1,nvars
         if (names(var) .eq. 'T') it = var    
         if (names(var) .eq. 'LATITCRS') ilat = var
         if (names(var) .eq. 'MAPFACCR') imap = var            
         if (names(var) .eq. 'H') ig = var         
      end do
C  If T, H or LATITCRS and MAPFACC not found, return error 1
      if (it .eq. -1 .or. ig .eq. -1
     $    .or. ilat .eq. -1 .or. imap. eq. -1) then
         print *, "Couldn't find T, H  or LATITCRS and 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

      
	do i=1,maxnl
	p(i)=1012.5*exp(vert_args(i)/-7.2)*100.
	enddo

        xincy=2.*ds

C  Bucle sobre todos los niveles de presion
      do 100 il=1,maxnl

C CALCULO DE LA COMPONENTE Q1 DEL VECTOR Q (en unidades de 1E-10 SI):

         do ic=3,nc-2
           do ir=3,nr-2
		lat=ingrid(ir,ic,1,ilat)*Pi/180.
		cor=4.*Pi/Per*Sin(lat)
		xincx=2.*ds
	
          B20=(-ingrid(ir+1,ic+1,il,ig)+ingrid(ir-1,ic+1,il,ig)+
     & ingrid(ir+1,ic-1,il,ig)-ingrid(ir-1,ic-1,il,ig))/(xincx*xincy)
        B20=B20*(ingrid(ir,ic+1,il,it)-ingrid(ir,ic-1,il,it))/xincx
        B20=B20*(ingrid(ir,ic,1,imap)**3.)        
          B21=(ingrid(ir,ic+2,il,ig)+ingrid(ir,ic-2,il,ig)-
     &    2.*ingrid(ir,ic,il,ig))/xincx**2
	 B21=B21*(-ingrid(ir+1,ic,il,it)+ingrid(ir-1,ic,il,it))/xincy
	 B21=B21*(ingrid(ir,ic,1,imap)**3.)	 
	  Q1=g*(B20-B21)*1.e10/cor
	 outgrid(ir,ic,il)=Q1
         enddo
        enddo 

        Do ir=3,nr-2
        outgrid(ir,2,il)=2*outgrid(ir,3,il)-outgrid(ir,4,il)
         outgrid(ir,1,il)=2*outgrid(ir,2,il)-outgrid(ir,3,il)
        outgrid(ir,nc-1,il)=2*outgrid(ir,nc-2,il)-outgrid(ir,nc-3,il)
         outgrid(ir,nc,il)=2*outgrid(ir,nc-1,il)-outgrid(ir,nc-2,il)
        Enddo
        Do ic=1,nc
         outgrid(2,ic,il)=2*outgrid(3,ic,il)-outgrid(4,ic,il)
         outgrid(1,ic,il)=2*outgrid(2,ic,il)-outgrid(3,ic,il)
         outgrid(nr-1,ic,il)=2*outgrid(nr-2,ic,il)-outgrid(nr-3,ic,il)
         outgrid(nr,ic,il)=2*outgrid(nr-1,ic,il)-outgrid(nr-2,ic,il)
        Enddo

100      end do
      USERFUNC = 0
      RETURN
      END
