C VAGIR.f

C VIS-5D analysis function to compute the V-Component of ageostrophic 
C irrotational wind (in m/s)

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:
c      include "/usr/local/src/vis5d/vis5d-4.2/src/v5df.h"
      integer ir,ic,il
      integer var,ig,iu,ug,vg,iv,ivagnd,vag,uag,ilat,imap
      real g,lat,Pi,Rt,Per,cor,uderiv,vderiv,ugeo,vgeo,par,ds
      real und(200,200,30),vage(200,200,30)
      real vnd(200,200,30),dtini(200,200,30),med(30)
      real vortini(200,200,30),max(30),max2(30)
      real dt,medini(30),vort,urd,vrd,urr,vrr
      integer control(30)
      parameter (Pi=3.141596,Rt=6371227.,Per=86400.,g=9.807,par=0.05)
      
C     Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0

C     Find the variables
      ig = -1
      iu = -1
      iv = -1
      ug = -1
      vg = -1
      ivagnd=-1
      vag=-1
      uag=-1
      imap=-1
      ilat=-1     
      do var=1,nvars
	 if (names(var) .eq. 'UAG') uag = var
	 if (names(var) .eq. 'VAG') vag = var
	 if (names(var) .eq. 'VAGND') ivagnd = var
	 if (names(var) .eq. 'UG') ug = var
         if (names(var) .eq. 'VG') vg = var
         if (names(var) .eq. 'H') ig = 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 LATITCRS or MAPFACCR not found, return error 1
      if (ilat.eq.-1 .or. imap.eq.-1) then
         print *, "Couldn't find 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

C If ugeo or vgeo variables are not found, we calculate them
	if ((vag.eq.-1).or.(uag.eq.-1)) then
	if ((ug.eq.-1).or.(vg.eq.-1)) then
	print*,'Calculating geostrophic wind internally'
C     If Geo not found, return error 1
      if ((ig.eq.-1).or.(iu.eq.-1).or.(iv.eq.-1)) then
         print *, "Couldn't find GEO or U variables!"
         userfunc = 1
         return
      endif

C     Compute ageostrophic wind
       do 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)
	    uderiv=(ingrid(ir-1,ic,il,ig)-ingrid(ir+1,ic,il,ig))/
     & (2.*ds)
            uderiv=uderiv*ingrid(ir,ic,1,imap)         
	    vderiv=(ingrid(ir,ic+1,il,ig)-ingrid(ir,ic-1,il,ig))/
     & (2.*ds)
            vderiv=vderiv*ingrid(ir,ic,1,imap)      
	    ugeo=-uderiv*g/cor
	    vgeo=vderiv*g/cor		
               und(ir,ic,il) = ingrid(ir,ic,il,iu)-ugeo
               vnd(ir,ic,il) = ingrid(ir,ic,il,iv)-vgeo
            end do
         end do
        do ir=2,nr-1
	   und(ir,1,il)=2.*und(ir,2,il)-und(ir,3,il)
	   und(ir,nc,il)=2.*und(ir,nc-1,il)-und(ir,nc-2,il)
	   vnd(ir,1,il)=2.*vnd(ir,2,il)-vnd(ir,3,il)
	   vnd(ir,nc,il)=2.*vnd(ir,nc-1,il)-vnd(ir,nc-2,il)
	enddo

	do ic=2,nc-1
	  und(1,ic,il)=2.*und(2,ic,il)-und(3,ic,il)
	  und(nr,ic,il)=2.*und(nr-1,ic,il)-und(nr-2,ic,il)
	  vnd(1,ic,il)=2.*vnd(2,ic,il)-vnd(3,ic,il)
	  vnd(nr,ic,il)=2.*vnd(nr-1,ic,il)-vnd(nr-2,ic,il)
	enddo
	und(1,1,il)=2.*und(2,2,il)-und(3,3,il)
	und(1,nc,il)=2.*und(2,nc-1,il)-und(3,nc-2,il)
       und(nr,1,il)=2.*und(nr-1,2,il)-und(nr-2,3,il)
       und(nr,nc,il)=2.*und(nr-1,nc-1,il)-und(nr-2,nc-2,il)
	vnd(1,1,il)=2.*vnd(2,2,il)-vnd(3,3,il)
	vnd(1,nc,il)=2.*vnd(2,nc-1,il)-vnd(3,nc-2,il)
       vnd(nr,1,il)=2.*vnd(nr-1,2,il)-vnd(nr-2,3,il)
       vnd(nr,nc,il)=2.*vnd(nr-1,nc-1,il)-vnd(nr-2,nc-2,il)
       end do
C If uageo or vageo has been found we use it and dont recalculate
	else
	print*,'Using geostrophic calculated winds'
	do il=1,maxnl
	 do ic=1,nc
	    do ir=1,nr
          und(ir,ic,il) =ingrid(ir,ic,il,iu)-ingrid(ir,ic,il,ug)
          vnd(ir,ic,il) =ingrid(ir,ic,il,iv)-ingrid(ir,ic,il,vg)
	  end do
	end do
	end do
	endif
	else
	print*,'Using ageostrophic calculated winds'
	do il=1,maxnl
	 do ic=1,nc
	    do ir=1,nr
          und(ir,ic,il) =ingrid(ir,ic,il,uag)
          vnd(ir,ic,il) =ingrid(ir,ic,il,vag)
	  end do
	end do
	end do
	endif
	do il=1,maxnl
	do ic=1,nc
	do ir=1,nr
	vage(ir,ic,il)=vnd(ir,ic,il)
	enddo
	enddo
	enddo
	if (ivagnd.eq.-1) then
	print*,'Calculating non-divergent ageostrophic wind internally'
C Calcul de la vorticitat inicial necessaria en tot el programa
	call vortic(nr,nc,maxnl,und,vnd,vortini,ds)
	
        Do il=1,maxnl
        Do ir=1,nr
        Do ic=1,nc 
         vortini(ir,ic,il)=vortini(ir,ic,il)*ingrid(ir,ic,1,imap)  
        Enddo
        Enddo
        Enddo	

C Calcul de la mitja del camp inicial a calcular
	call medlevels(vnd,medini,nr,nc,maxnl)
C Calcul de la divergencia inicial
	do il=1,maxnl
	call div(nr,nc,maxnl,und,vnd,dtini,ds,max,control,il)
           Do ir=1,nr
           Do ic=1,nc 
            dtini(ir,ic,il)=dtini(ir,ic,il)*ingrid(ir,ic,1,imap)  
           Enddo
           Enddo		
	enddo

C Calcul del 5% de control per a aturar la iteracio
	call maxlevls(dtini,max,nr,nc,maxnl)
	do il=1,maxnl
	max2(il)=max(il)
	max(il)=max(il)*par
	control(il)=-1
	enddo

C Inici de les iteracions

	do 36 il=1,maxnl

100         do ic=2,nc-1
          do ir=2,nr-1
C Calcul de la divergencia en el punt (ir,ic,il)
	dt=(und(ir,ic+1,il)-und(ir,ic-1,il))/
     & (2.*ds)+(vnd(ir-1,ic,il)-
     & vnd(ir+1,ic,il))/(2.*ds)
C Ajust de camp per fer-lo no divergent
	urd=-0.5*ds*dt
	vrd=-0.5*ds*dt
	und(ir,ic+1,il)=und(ir,ic+1,il)+urd
	und(ir,ic-1,il)=und(ir,ic-1,il)-urd
	vnd(ir+1,ic,il)=vnd(ir+1,ic,il)-vrd
	vnd(ir-1,ic,il)=vnd(ir-1,ic,il)+vrd
	enddo 
	enddo

        do ic=2,nc-1
          do ir=2,nr-1 
	vort=(vnd(ir,ic+1,il)-vnd(ir,ic-1,il))/
     & (2.*ds)-(und(ir-1,ic,il)-
     & und(ir+1,ic,il))/(2.*ds)
	urr=-0.5*ds*(vort-vortini(ir,ic,il))
	vrr=-0.5*ds*
     & (vort-vortini(ir,ic,il))
	und(ir-1,ic,il)=und(ir-1,ic,il)-urr
	und(ir+1,ic,il)=und(ir+1,ic,il)+urr
	vnd(ir,ic+1,il)=vnd(ir,ic+1,il)+vrr
	vnd(ir,ic-1,il)=vnd(ir,ic-1,il)-vrr
	enddo
	enddo
	call div(nr,nc,maxnl,und,vnd,dtini,ds,max,control,il)
 
        Do ir=1,nr
        Do ic=1,nc 
         dtini(ir,ic,il)=dtini(ir,ic,il)*ingrid(ir,ic,1,imap)  
        Enddo
        Enddo
	
	call maxlevls(dtini,max2,nr,nc,maxnl)
	if (control(il).eq.-1) goto 100

36	enddo

C Resituam els camps amb la mateixa mitja original
	call medlevels(vnd,med,nr,nc,maxnl)
	do il=1,maxnl
	 do ir=1,nr
	 do ic=1,nc
	    outgrid(ir,ic,il)=vage(ir,ic,il)-(vnd(ir,ic,il)-
     & med(il)+medini(il))
	 enddo
	enddo
	enddo
	else
	print*,'Using ageostrophic non-divergent calculated winds'
	do il=1,maxnl
	 do ir=1,nr
	 do ic=1,nc
	    outgrid(ir,ic,il)=vage(ir,ic,il)-ingrid(ir,ic,il,ivagnd)
	 enddo
	enddo
	enddo
	endif
      USERFUNC = 0
      RETURN
      END


C ---------------------------------------------------------------------------
C -- Calcula la divergencia del nivel il (2D)  de un campo vectorial 3D    --
C -- (campinu,campinv) y la pone en diver. Llena los contornos linealmente.--               --

      SUBROUTINE div(nr,nc,maxnl,campinu,campinv,diver,ds,epsil,
     &control,il)
	implicit none
      integer nr,nc,maxnl,il,ic,ir,control(il)
      real campinu(200,200,30),campinv(200,200,30),diver(200,200,30)
      real Pi,Rt,epsil(maxnl),ds
      parameter (Pi=3.141596,Rt=6371227.)
c      do 35 il=1,maxnl
	control(il)=1
         do ic=2,nc-1
            do ir=2,nr-1
	diver(ir,ic,il)=(campinu(ir,ic+1,il)-campinu(ir,ic-1,il))/
     & (2.*ds)+(campinv(ir-1,ic,il)-
     & campinv(ir+1,ic,il))/(2.*ds)
	    if (diver(ir,ic,il).gt.epsil(il)) control(il)=-1
           end do
         end do
C OMPLIM ELS CONTORNS LINEALMENT
	        do ir=2,nr-1
	   diver(ir,1,il)=2.*diver(ir,2,il)-diver(ir,3,il)
	   diver(ir,nc,il)=2.*diver(ir,nc-1,il)-diver(ir,nc-2,il)
	enddo

	do ic=2,nc-1
	  diver(1,ic,il)=2.*diver(2,ic,il)-diver(3,ic,il)
	  diver(nr,ic,il)=2.*diver(nr-1,ic,il)-diver(nr-2,ic,il)
	enddo
	diver(1,1,il)=2.*diver(2,2,il)-diver(3,3,il)
	diver(1,nc,il)=2.*diver(2,nc-1,il)-diver(3,nc-2,il)
       diver(nr,1,il)=2.*diver(nr-1,2,il)-diver(nr-2,3,il)
       diver(nr,nc,il)=2.*diver(nr-1,nc-1,il)-diver(nr-2,nc-2,il)
c35	end do
	end

C ----------------------------------------------------------------------------

	SUBROUTINE vortic(nr,nc,maxnl,campinu,campinv,vort,ds)
	implicit none
      integer nr,nc,maxnl,il,ic,ir
      real campinu(200,200,30),campinv(200,200,30),vort(200,200,30)
      real Pi,Rt,ds
      parameter (Pi=3.141596,Rt=6371227.)
      do 37 il=1,maxnl
         do ic=2,nc-1
            do ir=2,nr-1
	vort(ir,ic,il)=(campinv(ir,ic+1,il)-campinv(ir,ic-1,il))/
     & (2.*ds)-(campinu(ir-1,ic,il)-
     & campinu(ir+1,ic,il))/(2.*ds)
           end do
         end do
C OMPLIM ELS CONTORNS LINEALMENT
	        do ir=2,nr-1
	   vort(ir,1,il)=2.*vort(ir,2,il)-vort(ir,3,il)
	   vort(ir,nc,il)=2.*vort(ir,nc-1,il)-vort(ir,nc-2,il)
	enddo

	do ic=2,nc-1
	  vort(1,ic,il)=2.*vort(2,ic,il)-vort(3,ic,il)
	  vort(nr,ic,il)=2.*vort(nr-1,ic,il)-vort(nr-2,ic,il)
	enddo
	vort(1,1,il)=2.*vort(2,2,il)-vort(3,3,il)
	vort(1,nc,il)=2.*vort(2,nc-1,il)-vort(3,nc-2,il)
       vort(nr,1,il)=2.*vort(nr-1,2,il)-vort(nr-2,3,il)
       vort(nr,nc,il)=2.*vort(nr-1,nc-1,il)-vort(nr-2,nc-2,il)
37	end do
	end

C ----------------------------------------------------------------------------
C DONAT UN CAMP NRxNCxNL ens dona la mitja per nivells
	SUBROUTINE medlevels (campin,med,nr,nc,nl)
	implicit none
	integer nr,nc,nl,il,ic,ir
       real campin(200,200,30),med(nl)
	do il=1,nl
	   do ir=1,nr
	   do ic=1,nc
	     med(il)=med(il)+campin(ir,ic,il)
	   enddo
	   enddo
	   med(il)=med(il)/(nr*nc)
	enddo
	end
	
C ----------------------------------------------------------------------------
C DONAT UN CAMP NRxNCxNL ens dona el maxim per nivells
	SUBROUTINE maxlevls (campin,max,nr,nc,nl)
	implicit none
	integer nr,nc,nl,il,ic,ir
       real campin(200,200,30),max(nl)
	do il=1,nl
	   max(il)=campin(1,1,il)
	   do ir=1,nr
	   do ic=1,nc
c	print*,campin(ir,ic,il),max(il),il
	     if (campin(ir,ic,il).gt.max(il)) max(il)=campin(ir,ic,il)
	   enddo
	   enddo
	enddo
	end
	
