C ILIFTPA.f
C
C VIS-5D analysis function to compute Lifted Parcel Index (in C)

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	implicit none
C     ARGUMENTS:
      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, i, j, i500,it,ihum
      real p(22),ews,td,tcs,tetaws,tetaws_5
		
C     Specify number of levels in OUTGRID
      OUTNL = 1
      OUTLOWLEV = 0
      
	print*,'En unidades de °C'
	
C  Find the T, and RH variables
      it = -1
      ihum = -1
      igeo = -1
      do var=1,nvars
         if (names(var) .eq. 'T') it = var
         if (names(var) .eq. 'RH') ihum = var
      end do

C  If T, or RH not found, return error 1
      if (it .eq. -1 .or. ihum .eq. -1) then
         print *, "Couldn't find T, and/or RH variables!"
         userfunc = 1
         return
      endif
      
c Cercam p(i)=500. mb 
	do i=1,maxnl
	p(i)=1012.5*exp(vert_args(i)/-7.2)
	if (p(i).lt.501.and.p(i).gt.499) i500=i
	enddo
	
		DO 5770 J=1,nc
		DO 5770 I=1,nr

c	Entran T,RH,p y H

C	CALCULO VARIABLES EN CADA NIVEL PARA CADA GRID
	EWS=VAPFW(ingrid(i,j,1,it))
	TD=DPTOF(EWS*ingrid(i,j,1,ihum)/100.)

C	----------------------------------------------------------
C	------------------------------------------------------
C	TEMPERATURA POTENCIAL TERMOM HUMEDO PARTICULA SUPERFICIE
	TCS=TCONOF(ingrid(I,J,1,it),TD)
	TETAWS=ingrid(i,j,1,it)-WOBF(ingrid(i,j,1,it))+WOBF(TCS)
	tetaws_5=satlft(tetaws,500.)
	
	do il=1,outnl
		outgrid(i,j,il)=ingrid(i,j,i500,it)-tetaws_5
	enddo

5770 	continue
      USERFUNC = 0
      RETURN
      END



c ===========================================================================
C	ALGORITMO PARA CALCULAR LA PRESION DE VAPOR SATURANTE
C	RESPECTO AL AGUA EN hPa.
c	ENTRADA TEMPERATURA EN C.

	FUNCTION VAPFW(T)

	X=T
	IF(100.-X)3,4,4

3	X=X-273.16

C	AJUSTE DE LA CURVA ENTRE -50 < T < 100
4	POL=.99999683E-0+X*(-.90826951E-2+X*(.78736169E-4+
     1	X*(-.61117958E-6+X*(.43884187E-8+X*(-.29883885E-10+
     2	X*(.21874425E-12+X*(-.17892321E-14+X*(.11112018E-16+
     3	X*(-.30994571E-19)))))))))
	POL=POL*POL
	POL=POL*POL
	VAPFW=6.1078/(POL*POL)
	RETURN
	END
c ===========================================================================


c ===========================================================================
C	ALGORITMO PARA INVERTIR LA FORMULA DE TETENS
C	CALCULA LA TEMPERATURA EN GRADOS C DADA LA PRESION DE VAPOR
C	SATURANTE EN hPa

	FUNCTION DPTOF(EW)

	TOL=.0001
	IF (EW-.21382876E-9)20,20,30

20	DPTOF=-10000.
	RETURN

30	IF (1013.-EW)20,100,100

100	X=ALOG(EW/6.1078)
	BOT=17.269388-X
	DPTOF=(273.3*X)/BOT
	BOT=BOT*EW
200	EDP=VAPFW(DPTOF)
	DTDE=(DPTOF+237.7)/BOT
	DELT=DTDE*(EW-EDP)
	DPTOF=DPTOF+DELT

	IF(ABS(DELT)-TOL)300,300,200

300	DPTOF=DPTOF-TOL
	RETURN
	END
c ===========================================================================


c ===========================================================================
C	-------------------------------------------------------------
C	ALGORITMO PARA CALCULAR LA TEMPERATURA DE CONDENSACION
C	PO ELEVACION. ENTRADA TEMPERATURA Y PUNTO DE ROCIO AMBOS
C	EN C o K. LA TEMPERATURA DE CONDENSACION SE DEVUELVE EN C.


	FUNCTION TCONOF(TEMP,DEWPT)
	S=TEMP-DEWPT
	T=TEMP
	IF(100.0-TEMP)4,5,5
4	T=TEMP-273.16
C	AJUSTE DE UNA CURVA
5	DLT=S*(1.2185+.001278*T+S*(-.00219+11.73E-6*S-5.2E-6*T))
	TCONOF=T-DLT
	RETURN
	END
c ===========================================================================


c ===========================================================================
C	ALGORITMO PARA ELEVAR UNA PARTICULA A LO LARGO DE UNA 
C	PSEUDOADIABATICA.
C	CALCULA LA TEMPERATURA EN C A LA CUAL UNA DETERMINADA
C	PSEUDOADIABAICA EN C (THM) CRUZA UNA ISOBARA (P) EN hPa.

	FUNCTION SATLFT(THM,P)

	ROVCP=.28571428
	IF(ABS(P-1000.)-.001)100,100,200
100	SATLFT=THM
	RETURN

200	PWRP=(P/1000.)**ROVCP
	TONE=(THM+273.16)*PWRP-273.16

	EONE=WOBF(TONE)-WOBF(THM)
	RATE=1.
	GO TO 330
300	CONTINUE

	RATE=(TTWO-TONE)/(ETWO-EONE)
	TONE=TTWO
	EONE=ETWO
330	CONTINUE

	TTWO=TONE-EONE*RATE

	ETWO=(TTWO+273.16)/PWRP-273.16
	ETWO=ETWO+WOBF(TTWO)-WOBF(ETWO)-THM

	EOR=ETWO*RATE
	IF(ABS(EOR)-.1)400,400,300

400	SATLFT=TTWO-EOR
	RETURN
	END
c ===========================================================================


c ===========================================================================
C	----------------------------------------------------------
C	ALGORITMO PARA CALCULAR LA FUNCION DE WOBUS
C	ENTRADA TEMPERATURA EN C

	FUNCTION WOBF(T)
	
	X=T-20.
	IF(X) 10,10,20
10	CONTINUE

	POL=1.+X*(-8.8416605E-3+X*(1.4714143E-4+X*(-9.671989E-7
     #  +X*(-3.2607217E-8+X*(-3.8598073E-10)))))
	POL=POL*POL
	WOBF=15.130/(POL*POL)
	RETURN

20	POL=1.+X*(3.6182989E-3+X*(-1.3603273E-5+X*(4.9618922E-7
     #  +X*(-6.1059365E-9+X*(3.9401551E-11+X*(-1.2588129E-13
     #  +X*(1.668828E-16)))))))
	POL=POL*POL
	WOBF=29.930/(POL*POL)+.96*X-14.8
	RETURN
	END
c ===========================================================================
