C EPT.f
C
C VIS-5D analysis function to compute the equivalent 
c potencial temperature (in C)

      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 it, ihur
      real t, hur, the,dewpoint,hum
      real qv,p,td,tlden,the1,the2,rd,cp,tl
	      PARAMETER(rd=287.04,cp=1004.64)

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

C     Find the T and RH variables
      it = -1
      ihur = -1

      do var=1,nvars
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'RH') ihur = var
      end do

C     If T or RH not found, return error 1
      if (it .eq. -1 .or. ihur .eq. -1) then
         print *, "Couldn't find T or RH variables!"
         userfunc = 1
         return
      endif

C     Compute tpoteq
      do il=1,maxnl
         do ic=1,nc
            do ir=1,nr
               t = ingrid(ir,ic,il,it)
               hur = max(ingrid(ir,ic,il,ihur),0.1)
		t=t+273.16
c	Pressio en mil.libars
	       p=1012.5*exp(vert_args(il)/-7.2)
C              Check for missing data
               if (t.ge.1.0e30.or.hur.ge.1.0e30) then
		  print*,'Missing values added'
                  the = 1.0e35
               else
                 

c Calculo de la temperatura potencial equivalente 
           qv=hum(t,hur)/(p*100.)
           td=dewpoint(qv,p*100.)
           tlden=1./(td-56.)+log(t/td)/800.
           tl=1./tlden+56.
           the1=t*((1000./p)**(rd/cp*(1.-0.28*qv)))
           the2=exp(1000.*qv*(1.+0.81*qv)*(3.376/tl-0.00254)) 
           the=the1*the2
      
               end if
               outgrid(ir,ic,il) = the-273.16
            end do
         end do
      end do

      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
