C v_avg.f
C
C VIS-5D analysis function to compute VERTICAL AVERAGE of FIELD (units of
C FIELD) between specified pressure LEVELS

C See below ROMU: For selection of FIELD and pressure levels

C ==============================================================================

      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
      integer ivar,numil
      real p(22),pbot,ptop

C     Specify number of levels in OUTGRID
      OUTNL = 1
      OUTLOWLEV = 0
  
C ROMU: Set the bottom and top pressure levels (hPa)
      pbot=1000.
      ptop=700. 
  
C ROMU: Set the FIELD variable of interest
      ivar = -1
      do var=1,nvars
         if (names(var) .eq. 'FORZQfil') ivar = var
      end do

C If not found, return error 1
      if (ivar.eq.-1) then
         print *, "Couldn't find VAR variable!"
         userfunc = 1
         return
      endif

C Fill p(hPa)       
      do il=1,maxnl
        p(il)=1012.5*exp(vert_args(il)/-7.2)	
      enddo

C Vertical integration between pbot and ptop
      Do 1000 ir=1,nr
      Do 1000 ic=1,nc
        numil=0
        outgrid(ir,ic,OUTNL)=0.
        Do il=1,maxnl
          If (pbot.Ge.p(il).And.p(il).Ge.ptop) Then
            numil=numil+1
            outgrid(ir,ic,OUTNL)=outgrid(ir,ic,OUTNL)+
     $                           ingrid(ir,ic,il,ivar)
          Endif
        Enddo
        outgrid(ir,ic,OUTNL)=outgrid(ir,ic,OUTNL)/Real(numil)  
1000  Continue
      

      USERFUNC = 0
      RETURN
     
      END
