C AVOR.f
C
C VIS-5D analysis function to compute Absolute vorticity (in 1.E-5 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
      integer iu, iv,irv, ilat, imap
      real relv,lat,Pi,Rt,Per,absrel,ds
      parameter (Pi=3.141596,Rt=6371227.,Per=86400.)

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

C     Find the U, V, LATITCRS and MAPFACCR variables
      iu = -1
      iv = -1
      ilat = -1
      imap = -1
      irv= -1
      do var=1,nvars
         if (names(var) .eq. 'U') iu = var
         if (names(var) .eq. 'V') iv = var
         if (names(var) .eq. 'LATITCRS') ilat = var
         if (names(var) .eq. 'MAPFACCR') imap = var                  
	 if (names(var) .eq. 'RVOR') irv=var
      end do

C     If U, V, or LATITCRS or MAPFACCR not found, return error 1
      if (iu .eq. -1 .or. iv .eq. -1 .or. ilat .eq. -1
     $    .or. imap .eq. -1) then
         print *, "Couldn't find U, V, or LATITCRS, 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

	If (irv.eq.-1) then
	print*,'Calculating Relative vorticity internally'
C     Compute relative vorticity
      do il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
		lat=ingrid(ir,ic,1,ilat)*Pi/180.
		relv=(ingrid(ir,ic+1,il,iv)-ingrid(ir,ic-1,il,iv))/
     & (2.*ds)-(ingrid(ir-1,ic,il,iu)-ingrid(ir+1,ic,il,iu))/(2*ds)
                relv=relv*ingrid(ir,ic,1,imap)
		absrel=relv+4.*Pi/Per*Sin(lat)              		
               outgrid(ir,ic,il) = absrel*1.e+5
            end do
         end do
        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)
      end do
	
	else
	print*,'Reading calculated RVOR'
	 do il=1,maxnl
          do ic=1,nc
           do ir=1,nr
	lat=ingrid(ir,ic,1,ilat)*Pi/180.
        outgrid(ir,ic,il)=ingrid(ir,ic,il,irv)+4.*Pi/Per*Sin(lat)*1.e+5
            end do
         end do
	end do
	endif 
      USERFUNC = 0
      RETURN
      END
