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

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 )
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 ivg,iug,igavor,igrvor,ig
      real xincx,xincy,gradx,grady,Pi,lat,Rt,escal,ds
      real gavor(300,300,22),per,relv,ugeo(300,300,22)
      real vgeo(300,300,22),cor,lapla,g
      PARAMETER(Pi=3.141596,Rt=6371227.,Per=86400.,escal=1.e9,
     & g=9.81)   
C     Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0

	print*,'Reescalat en 1.e9'
C     Specify number of levels in OUTGRID
      OUTNL = MAXNL

C  Find the H variable
      ivg = -1
      iug = -1
      igavor = -1
      igrvor = -1
      ig = -1
      ilat = -1
      imap = -1

      do var=1,nvars
         if (NAMES(var).eq."GRVOR") igrvor = var
         if (NAMES(var).eq."GAVOR") igavor = var
         if (NAMES(var).eq."UG") iug = var
         if (NAMES(var).eq."VG") ivg = var
	 if (NAMES(var).eq."H") ig = var
         if (names(var) .eq. 'LATITCRS') ilat = var
         if (names(var) .eq. 'MAPFACCR') imap = var  	 
      end do
      
C  If H or LATITCRS or MAPFACCR not found, return error 1
      if (ig.eq.-1 .or. ilat.eq.-1 .or. imap.eq.-1) then
         print *, "Couldn't find H, 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
    
      	xincy=2.*ds
	if (igavor.eq.-1) then
	if (igrvor.eq.-1) then
	if (iug.eq.-1.or.ivg.eq.-1) then
C Calculo del viento geostrofico y GAVOR y adveccion
	print*,'Calculating UG, VG and GAVOR internally'
      do 100 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	    lat=ingrid(ir,ic,1,ilat)*Pi/180.
	    cor=4.*Pi/Per*Sin(lat)
	ugeo(ir,ic,il)=-g/cor*(ingrid(ir-1,ic,il,ig)-
     &ingrid(ir+1,ic,il,ig))/(2.*ds)
        ugeo(ir,ic,il)=ugeo(ir,ic,il)*ingrid(ir,ic,1,imap) 
	vgeo(ir,ic,il)=g/cor*(ingrid(ir,ic+1,il,ig)-
     & ingrid(ir,ic-1,il,ig))/(2.*ds)
        vgeo(ir,ic,il)=vgeo(ir,ic,il)*ingrid(ir,ic,1,imap)      	
            end do
         end do
         do ic=2,nc-1
            do ir=2,nr-1
	    lat=ingrid(ir,ic,1,ilat)*Pi/180.
	    cor=4.*Pi/Per*Sin(lat)
	    lapla=(ingrid(ir,ic+1,il,ig)+ingrid(ir,ic-1,il,ig)-2.*
     & ingrid(ir,ic,il,ig))/(ds**2)            
     & +(ingrid(ir-1,ic,il,ig)+ingrid(ir+1,ic,il,ig)-2.*
     & ingrid(ir,ic,il,ig))/(ds**2)
           lapla=lapla*(ingrid(ir,ic,1,imap)**2)
	    relv=lapla*g/cor
	  gavor(ir,ic,il)=relv+cor
	enddo
	enddo

         do ic=3,nc-2
            do ir=3,nr-2
	xincx=2.*ds
  	gradx=(gavor(ir,ic+1,il)-gavor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)
	grady=(gavor(ir-1,ic,il)-gavor(ir+1,ic,il))/xincy
  	grady=grady*ingrid(ir,ic,1,imap)	
        outgrid(ir,ic,il) =-(ugeo(ir,ic,il)*gradx+
     &  vgeo(ir,ic,il)*grady)*escal
            end do
         end do
100	enddo
	else
C Lectura del viento geostrofico,calculo de GAVOR i adveccion
	print*,'Using calculated UG, VG and calculating GAVOR 
     & internally'
      do 101 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	ugeo(ir,ic,il)=ingrid(ir,ic,il,iug)
	vgeo(ir,ic,il)=ingrid(ir,ic,il,ivg)	
            end do
         end do
         do ic=2,nc-1
            do ir=2,nr-1
	    lat=ingrid(ir,ic,1,ilat)*Pi/180.
	    cor=4.*Pi/Per*Sin(lat)
	    lapla=(ingrid(ir,ic+1,il,ig)+ingrid(ir,ic-1,il,ig)-2.*
     & ingrid(ir,ic,il,ig))/(ds**2)
     & +(ingrid(ir-1,ic,il,ig)+ingrid(ir+1,ic,il,ig)-2.*
     & ingrid(ir,ic,il,ig))/(ds**2)
           lapla=lapla*(ingrid(ir,ic,1,imap)**2)
	    relv=lapla*g/cor
	  gavor(ir,ic,il)=relv+cor
	enddo
	enddo

         do ic=3,nc-2
            do ir=3,nr-2
	xincx=2.*ds
  	gradx=(gavor(ir,ic+1,il)-gavor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)
	grady=(gavor(ir-1,ic,il)-gavor(ir+1,ic,il))/xincy
	grady=grady*ingrid(ir,ic,1,imap)
        outgrid(ir,ic,il) =-(ugeo(ir,ic,il)*gradx+
     &  vgeo(ir,ic,il)*grady)*escal
            end do
         end do
101	enddo
	endif
c Else del if de igrvor
	else
c Tenemos GRVOR pero hay que ver si tenemos UG y VG
	if (iug.eq.-1.or.ivg.eq.-1) then
C Calculo del viento geostrofico y GAVOR(a partir de RGVOR) y adveccion
	print*,'Using calculated RGVOR and calculating UG, VG & GAVOR
     & internally'
      do 102 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	    lat=ingrid(ir,ic,1,ilat)*Pi/180.
	    cor=4.*Pi/Per*Sin(lat)
	ugeo(ir,ic,il)=-g/cor*(ingrid(ir-1,ic,il,ig)-
     &ingrid(ir+1,ic,il,ig))/(2.*ds)
        ugeo(ir,ic,il)=ugeo(ir,ic,il)*ingrid(ir,ic,1,imap)      
     	vgeo(ir,ic,il)=g/cor*(ingrid(ir,ic+1,il,ig)-
     & ingrid(ir,ic-1,il,ig))/(2.*ds)
        vgeo(ir,ic,il)=vgeo(ir,ic,il)*ingrid(ir,ic,1,imap) 	
	  gavor(ir,ic,il)=ingrid(ir,ic,il,igrvor)/1.e5+cor
	enddo
	enddo
         do ic=3,nc-2
            do ir=3,nr-2
	xincx=2.*ds
  	gradx=(gavor(ir,ic+1,il)-gavor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)  	
	grady=(gavor(ir-1,ic,il)-gavor(ir+1,ic,il))/xincy
  	grady=grady*ingrid(ir,ic,1,imap)	
        outgrid(ir,ic,il) =-(ugeo(ir,ic,il)*gradx+
     &  vgeo(ir,ic,il)*grady)*escal
            end do
         end do
102	enddo
	else
C Lectura del viento geostrofico,calculo de AVOR(a partir de RGVOR) y adveccion
	print*,'Using calculated UG, VG & RGVOR and calculating GAVOR
     &internally'
      do 103 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	    lat=ingrid(ir,ic,1,ilat)*Pi/180.
	    cor=4.*Pi/Per*Sin(lat)
	ugeo(ir,ic,il)=ingrid(ir,ic,il,iug)
	vgeo(ir,ic,il)=ingrid(ir,ic,il,ivg)
	gavor(ir,ic,il)=ingrid(ir,ic,il,igrvor)/1.e5+cor
            end do
         end do
         do ic=3,nc-2
            do ir=3,nr-2
	xincx=2.*ds
  	gradx=(gavor(ir,ic+1,il)-gavor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)  	  	
	grady=(gavor(ir-1,ic,il)-gavor(ir+1,ic,il))/xincy
	grady=grady*ingrid(ir,ic,1,imap)  	
        outgrid(ir,ic,il) =-(ugeo(ir,ic,il)*gradx+
     &  vgeo(ir,ic,il)*grady)*escal
            end do
         end do
103	enddo
	endif
	endif
c ELSE del if de igavor
	else
c Tenemos GAVOR pero hay que ver si tenemos UG y VG
	if (iug.eq.-1.or.ivg.eq.-1) then
C Lectura de GAVOR i calculo del viento geostrofico y adveccion
	print*,'Using calculated GAVOR and calculating UG & VG'
      do 104 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	    lat=ingrid(ir,ic,1,ilat)*Pi/180.
	    cor=4.*Pi/Per*Sin(lat)
	ugeo(ir,ic,il)=-g/cor*(ingrid(ir-1,ic,il,ig)-
     &ingrid(ir+1,ic,il,ig))/(2.*ds)     
        ugeo(ir,ic,il)=ugeo(ir,ic,il)*ingrid(ir,ic,1,imap)      
	vgeo(ir,ic,il)=g/cor*(ingrid(ir,ic+1,il,ig)-
     & ingrid(ir,ic-1,il,ig))/(2.*ds)
        vgeo(ir,ic,il)=vgeo(ir,ic,il)*ingrid(ir,ic,1,imap)         
	gavor(ir,ic,il)=ingrid(ir,ic,il,igavor)/1.e5	
            end do
         end do
         do ic=3,nc-2
            do ir=3,nr-2
	xincx=2.*ds
  	gradx=(gavor(ir,ic+1,il)-gavor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)  	  	
	grady=(gavor(ir-1,ic,il)-gavor(ir+1,ic,il))/xincy
  	grady=grady*ingrid(ir,ic,1,imap)  		
        outgrid(ir,ic,il) =-(ugeo(ir,ic,il)*gradx+
     &  vgeo(ir,ic,il)*grady)*escal
            end do
         end do
104	enddo
	else
C Lectura del viento geostrofico, de AVOR y calculo de la adveccion
	print*,'Using calculated GAVOR, UG and VG'
      do 105 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	ugeo(ir,ic,il)=ingrid(ir,ic,il,iug)
	vgeo(ir,ic,il)=ingrid(ir,ic,il,ivg)
	gavor(ir,ic,il)=ingrid(ir,ic,il,igavor)/1.e5	
            end do
         end do
         do ic=3,nc-2
            do ir=3,nr-2
	xincx=2.*ds
  	gradx=(gavor(ir,ic+1,il)-gavor(ir,ic-1,il))/xincx
  	gradx=gradx*ingrid(ir,ic,1,imap)  	  	
	grady=(gavor(ir-1,ic,il)-gavor(ir+1,ic,il))/xincy
  	grady=grady*ingrid(ir,ic,1,imap)  		
        outgrid(ir,ic,il) =-(ugeo(ir,ic,il)*gradx+
     &  vgeo(ir,ic,il)*grady)*escal
            end do
         end do
105	enddo
	endif	
	endif
	do il=1,maxnl
        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
