C AVAD.f
C
C VIS-5D analysis function to compute Absolute vorticity advection
C (in 1.E-9 s-2) 

      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,ilat,imap
      integer iavor,iu,iv,irvor
      real xincx,xincy,gradx,grady,Pi,lat,Rt,escal,ds
      real avor(300,300,22),per,relv
      PARAMETER(Pi=3.141596,Rt=6371227.,Per=86400.,escal=1.e9)   

	print*,'Reescalat en 1.e9'
C     Specify number of levels in OUTGRID
      OUTLOWLEV = 0
      OUTNL = MAXNL
      
C     Find the U, V, LATITCRS and MAPFACCR variables
      iu = -1
      iv = -1
      ilat = -1
      imap = -1
      iavor = -1
      irvor = -1

      do var=1,nvars
         if (NAMES(var).eq."RVOR") irvor = var
         if (NAMES(var).eq."AVOR") iavor = var
         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            
      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 (iavor.eq.-1) then
	if (irvor.eq.-1) then
	
c Calculo interno de AVOR
	print*,'Calculating RVOR&AVOR internally'
	do 365 il=1,maxnl
	do 365 ir=2,nr-1
	lat=ingrid(ir,ic,1,ilat)*Pi/180.
	do 365 ic=2,nc-1

	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)
	avor(ir,ic,il)=relv+4.*Pi/Per*Sin(lat)
365	enddo


	else
	print*,'Using calculated RVOR'
	do 368 ir=2,nr-1
	lat=ingrid(ir,ic,1,ilat)*Pi/180.
	do 368 ic=2,nc-1
	do 368 il=1,maxnl
	avor(ir,ic,il)=ingrid(ir,ic,il,irvor)/1.e5+4.*Pi/Per*Sin(lat)
368	enddo
	endif
	
	else
c Lectura de AVOR ya calculado
	print*,'Using calculated AVOR'
	do 345 ir=2,nr-1
	do 345 ic=2,nc-1
	do 345 il=1,maxnl
	avor(ir,ic,il)=ingrid(ir,ic,il,iavor)/1.e5
345	enddo
	endif
	
C Compute Absolute vorticity Advection
	xincy=2.*ds
      do il=1,maxnl
         do ic=3,nc-2
            do ir=3,nr-2
	lat=ingrid(ir,ic,1,ilat)*Pi/180.
	xincx=2.*ds
  	gradx=(avor(ir,ic+1,il)-avor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)
	grady=(avor(ir-1,ic,il)-avor(ir+1,ic,il))/xincy
  	grady=grady*ingrid(ir,ic,1,imap)	
        outgrid(ir,ic,il) =-(ingrid(ir,ic,il,iu)*gradx+
     &  ingrid(ir,ic,il,iv)*grady)*escal
            end do
         end do
        do ir=2,nr-1
	 outgrid(ir,2,il)=3.*outgrid(ir,3,il)-3.*outgrid(ir,4,il)+
     &  outgrid(ir,5,il)
	 outgrid(ir,1,il)=6.*outgrid(ir,3,il)-8.*outgrid(ir,4,il)+
     &  3.*outgrid(ir,5,il)
	 outgrid(ir,ic-1,il)=3.*outgrid(ir,ic-2,il)-3.*
     &  outgrid(ir,ic-3,il)+outgrid(ir,ic-4,il)
	 outgrid(ir,ic,il)=6.*outgrid(ir,ic-2,il)-8.*
     &  outgrid(ir,ic-3,il)+3.*outgrid(ir,ic-4,il)
	enddo
	do ic=2,nc-1
	outgrid(2,ic,il)=3.*outgrid(3,ic,il)-3.*outgrid(4,ic,il)+
     &  outgrid(5,ic,il)
	outgrid(1,ic,il)=6.*outgrid(3,ic,il)-8.*outgrid(4,ic,il)+
     &  3.*outgrid(5,ic,il)
	outgrid(nr-1,ic,il)=3.*outgrid(nr-2,ic,il)-3.*
     &  outgrid(nr-3,ic,il)+outgrid(nr-4,ic,il)
	outgrid(nr,ic,il)=6.*outgrid(nr-2,ic,il)-8.*
     &  outgrid(nr-3,ic,il)+3.*outgrid(nr-4,ic,il)
	enddo
	outgrid(1,1,il)=3.*outgrid(2,2,il)-3.*outgrid(3,3,il)+
     &  outgrid(4,4,il)
	outgrid(nr,nc,il)=3.*outgrid(nr-1,nc-1,il)-3.*
     &  outgrid(nr-2,nc-2,il)+outgrid(nr-3,nc-3,il)
	outgrid(1,nc,il)=3.*outgrid(2,nc-1,il)-3.*outgrid(3,nc-2,il)+
     &  outgrid(4,nc-3,il)
	outgrid(nr,1,il)=3.*outgrid(nr-1,2,il)-3.*outgrid(nr-2,3,il)+
     &  outgrid(nr-3,4,il)
      end do

      USERFUNC = 0
      RETURN
      END
