C CI1050.f
C
C VIS-5D analysis function to compute the convective inestability
C (EPT500 - EPT1000, in degrees)

      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,il1000,il500
      integer it, ihur,itpe
      real t, hur,dewpoint,hum
      real qv,p,td,tlden,the1,the2,rd,cp,tl,p1000,p500,the1000,the500
      PARAMETER(rd=287.04,cp=1004.64)

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

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

      do var=1,nvars
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'RH') ihur = var
         if (names(var) .eq. 'EPT') itpe = 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   Find pressure levels
      do il=1,maxnl
      
       p=1012.5*exp(vert_args(il)/-7.2)
       If ((p.Gt.999.).And.(p.Lt.1001.)) Then
        il1000=il
        p1000=p
       Endif
       If ((p.Gt.499.).And.(p.Lt.501.)) Then
        il500=il
        p500=p
       Endif
      enddo

	if (itpe.eq.-1) then
      print*,'Calculating Potencial equivalent temperature internally'
         do ic=1,nc
            do ir=1,nr
c Calculs per 1000 mb	    
               t = ingrid(ir,ic,il1000,it)
               hur = ingrid(ir,ic,il1000,ihur)
               t=t+273.16

c Calculo de la temperatura potencial equivalente 
           qv=hum(t,hur)/(p1000*100.)
           td=dewpoint(qv,p1000*100.)
           tlden=1./(td-56.)+log(t/td)/800.
           tl=1./tlden+56.
           the1=t*((1000./p1000)**(rd/cp*(1.-0.28*qv)))
           the2=exp(1000.*qv*(1.+0.81*qv)*(3.376/tl-0.00254)) 
           the1000=the1*the2

c Calculs per 500 mb	    
               t = ingrid(ir,ic,il500,it)
               hur = ingrid(ir,ic,il500,ihur)
               t=t+273.16

c Calculo de la temperatura potencial equivalente 
           qv=hum(t,hur)/(p500*100.)
           td=dewpoint(qv,p500*100.)
           tlden=1./(td-56.)+log(t/td)/800.
           tl=1./tlden+56.
           the1=t*((1000./p500)**(rd/cp*(1.-0.28*qv)))
           the2=exp(1000.*qv*(1.+0.81*qv)*(3.376/tl-0.00254)) 
           the500=the1*the2      
	

           do il=1,OUTNL  
               outgrid(ir,ic,il) = the500-the1000
           enddo        

         end do
      end do

	else
	print*,'Reading calculated EPT'
        do il=1,OUTNL  
	 do ic=1,nc
          do ir=1,nr
      outgrid(ir,ic,il)=ingrid(ir,ic,il500,itpe)-
     &ingrid(ir,ic,il1000,itpe)
	enddo
	enddo
	enddo
	endif

      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
