C EPV.f
C
C VIS-5D analysis function to compute Equivalent Potential Vorticity
C (in PVU)

      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, it,ipt,ihur,ilat,imap
      real pv,lat,Pi,Rt,Per,cor,vx,thetax,uy,thetay,thetap,up,vp
      real p(22),rocp,g,xincx,xincy,tpot(300,300,22),exnf
      real t,hur,pl,qv,td,tlden,the1,the2,tl,dewpoint,hum,ds
      parameter (Pi=3.141596,Rt=6371227.,Per=86400.,rocp=2./7.,g=9.807)

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

C Find the U, V, T, RH, LATITCRS and MAPFACCR variables
      iu = -1
      iv = -1
      it = -1
      ihur=-1
      ilat = -1
      imap = -1     
      ipt= -1
      do var=1,nvars
         if (names(var) .eq. 'U') iu = var
         if (names(var) .eq. 'V') iv = var
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'RH') ihur = var
         if (names(var) .eq. 'LATITCRS') ilat = var
         if (names(var) .eq. 'MAPFACCR') imap = var                 
         if (names(var) .eq. 'EPT') ipt = var
      end do

C  If U, V, T, RH, or LATITCRS and MAPFACCR not found, return error 1
      if (iu.eq.-1.or.iv.eq.-1.or.it.eq.-1.or.ihur.eq.-1
     $    .or. ilat.eq.-1 .or. imap.eq.-1) then
         print *, "Couldn't find U, V, T, RH, LATITCRS and MAPFACCR !"
         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
     
      
       do il=1,maxnl
	p(il)=1012.5*exp(vert_args(il)/-7.2)*100.
	enddo
	if (ipt.eq.-1) then
	print*,'Calculating EPT internally'
	do 123 il=1,maxnl
	do 123 ic=1,nc
	do 123 ir=1,nr
	t=ingrid(ir,ic,il,it)+273.16
	hur=ingrid(ir,ic,il,ihur)
	pl=p(il)
           qv=hum(t,hur)/pl
           td=dewpoint(qv,pl)
           tlden=1./(td-56.)+log(t/td)/800.
           tl=1./tlden+56.
           the1=t*((100000./pl)**(rocp*(1.-0.28*qv)))
           the2=exp(1000.*qv*(1.+0.81*qv)*(3.376/tl-0.00254)) 
           tpot(ir,ic,il)=the1*the2
123	enddo
	else
	print*,'Using calculated EPT'
	do 124 il=1,maxnl
	do 124 ic=1,nc
	do 124 ir=1,nr
          tpot(ir,ic,il)=ingrid(ir,ic,il,ipt)
124	enddo
	endif

C Empieza el calculo de PV
 
        xincy=2.*ds


c Calculo del PV en pv units (10E-6 SI)
        do 105 il=1,maxnl
         exnf=(100000./p(il))**rocp   
         do ic=2,nc-1
           do ir=2,nr-1
		lat=ingrid(ir,ic,1,ilat)*Pi/180.
		xincx=2.*ds
		cor=4.*Pi/Per*Sin(lat) 
             vx=(ingrid(ir,ic+1,il,iv)-ingrid(ir,ic-1,il,iv))/xincx
             vx=vx*ingrid(ir,ic,1,imap)             
	     thetax=(tpot(ir,ic+1,il)-tpot(ir,ic-1,il))/xincx
             thetax=thetax*ingrid(ir,ic,1,imap)	     
	     uy=(ingrid(ir-1,ic,il,iu)-ingrid(ir+1,ic,il,iu))/xincy
             uy=uy*ingrid(ir,ic,1,imap)	     
	     thetay=(tpot(ir-1,ic,il)-tpot(ir+1,ic,il))/xincy
             thetay=thetay*ingrid(ir,ic,1,imap)	     
	if (il.eq.1) then
	     up=(ingrid(ir,ic,1,iu)-ingrid(ir,ic,2,iu))/(p(1)-p(2))
	     vp=(ingrid(ir,ic,1,iv)-ingrid(ir,ic,2,iv))/(p(1)-p(2))
	     thetap=(tpot(ir,ic,1)-tpot(ir,ic,2))/(p(1)-p(2))
	    else
	    if (il.eq.maxnl) then
	     up=(ingrid(ir,ic,maxnl-1,iu)-ingrid(ir,ic,maxnl,iu))/
     &       (p(maxnl-1)-p(maxnl))
	     vp=(ingrid(ir,ic,maxnl-1,iv)-ingrid(ir,ic,maxnl,iv))/
     &       (p(maxnl-1)-p(maxnl))
	     thetap=(tpot(ir,ic,maxnl-1)-tpot(ir,ic,maxnl))/
     &       (p(maxnl-1)-p(maxnl))
	    else
	     up=((p(il)-p(il+1))*(ingrid(ir,ic,il-1,iu)-
     & ingrid(ir,ic,il,iu))/(p(il-1)-p(il))+(p(il-1)-p(il))*
     & (ingrid(ir,ic,il,iu)-ingrid(ir,ic,il+1,iu))/(p(il)-p(il+1)))/
     & (p(il-1)-p(il+1))
	     vp=((p(il)-p(il+1))*(ingrid(ir,ic,il-1,iv)-
     & ingrid(ir,ic,il,iv))/(p(il-1)-p(il))+(p(il-1)-p(il))*
     & (ingrid(ir,ic,il,iv)-ingrid(ir,ic,il+1,iv))/(p(il)-p(il+1)))/
     & (p(il-1)-p(il+1))
	     thetap=((p(il)-p(il+1))*(tpot(ir,ic,il-1)-
     & tpot(ir,ic,il))/(p(il-1)-p(il))+(p(il-1)-p(il))*
     & (tpot(ir,ic,il)-tpot(ir,ic,il+1))/(p(il)-p(il+1)))/
     & (p(il-1)-p(il+1))
		endif
		endif
           pv=-g*(cor+vx-uy)*thetap+g*(thetax*vp-thetay*up)
	  outgrid(ir,ic,il) = pv*1.E6
          enddo
        enddo 
     

        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(nr,nc,il)=2.*outgrid(nr-1,nc-1,il)-outgrid(nr-2,nc-2,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)

105      enddo
      USERFUNC = 0
      RETURN
      END


      function hum(T,H)
       real t,h,hum          
       hum=380.042*exp(17.27*(T-273.16)/(T-35.86))*H/100.
       Return
      end

      function dewpoint(qv,P)
        real qv,p,bb,cc,dd,qvpln,dewpoint           
       parameter(bb=17.27,cc=17.27*273.16,dd=35.86)
       qvpln=log(qv*P/380.042)
       dewpoint=(qvpln*dd-cc)/(qvpln-bb) 
       Return
      end




