C
C CAPEN.f

C Function to compute the inhibition energy (negative area) (in J/kg)


      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, il
      integer it,ihum,igeo
      real p(22)

c VARIABLES INTERPOLADAS PARA EL SONDEO EN INCREMENTOS DE 10 hPa
	DIMENSION PP(105),HP(105),TP(105),TDP(105),QP(105),TETAP(105)
     $            ,TETAWP(105),TETAPART(105),TEPOTEQP(105),
     $		  TCP(105),PCP(105)
c VARIABLES CALCULADAS SOBRE LOS NIVELES DE PRESION DE REFERENCIA
	DIMENSION TETA(22),TD(22),Q(22)
     	DIMENSION TEPOTEQ(300,300,22)


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

C     Find the T, RH and H 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
         if (names(var) .eq. 'H') igeo = var
      end do

C     If T, RH or H not found, return error 1
      if (it .eq. -1 .or. ihum .eq. -1 .or. igeo .eq. -1) then
         print *, "Couldn't find T, RH, and/or H variables!"
         userfunc = 1
         return
      endif
c Cercam p(i)=1000. mb i p(i)=700 mb

	do i=1,maxnl
	p(i)=1012.5*exp(vert_args(i)/-7.2)
	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

		DO 1101 K=1,maxnl
	TETA(K)=((273.16+ingrid(i,j,k,it))*(1000./P(K))**.28571482)
     &	-273.16
	EWS=VAPFW(ingrid(i,j,k,it))
	TD(K)=DPTOF(EWS*ingrid(i,j,k,ihum)/100.)
	Q(K)=.622*(VAPFW(TD(K))/(P(K)-VAPFW(TD(K))))
	TL=TCONOF(ingrid(i,j,k,iT),TD(K))+273.16
	TEPOTEQ(I,J,K)=((ingrid(I,J,K,iT)+273.16)*((1000./P(K))**(.2854*
     $(1.-.28*Q(K)))*EXP(((3.376/TL)-.00254)*1E3*Q(K)*(1.+.81*Q(K)))))
     $-273.16
 1101		CONTINUE

C	-------------------------------------------------------
C	INT ES LA INSTRUCCION FORTRAN PARA OBTENER EL ENTERO INFERIOR
C	NINT ES LA INSTRUCCION PARA OBTENER EL ENTERO MAS PROXIMO
C	--------------------------------------------------------
C	...............................................................
C	INTERPOLACION A INTERVALOS DE 10 MB

			DO 733 K=2,maxnl
	KK1=NINT(P(K-1)/10.)
	KK2=NINT(P(K)/10.)
		DO 732 KK=KK1,KK2,-1
	PP(KK)=FLOAT(KK)*10.
	TP(KK)=ingrid(I,J,K-1,it)+(ingrid(I,J,K,it)-ingrid(I,J,K-1,it))*
     &  LOG(P(K-1)/10./FLOAT(KK))/LOG(P(K-1)/P(K))
	TDP(KK)=TD(K-1)+(TD(K)-TD(K-1))*LOG(P(K-1)/10./FLOAT(KK))/
     #         LOG(P(K-1)/P(K))
	QP(KK)=Q(K-1)+(Q(K)-Q(K-1))*LOG(P(K-1)/10./FLOAT(KK))/
     #         LOG(P(K-1)/P(K))
	TEPOTEQP(KK)=TEPOTEQ(I,J,K-1)+(TEPOTEQ(I,J,K)-TEPOTEQ(I,J,K-1))
     #         *LOG(P(K-1)/10./FLOAT(KK))/LOG(P(K-1)/P(K))
	TETAP(KK)=((273.16+TP(KK))*(1000./PP(KK))**.28571482)-273.16
	TCP(KK)=TCONOF(TP(KK),TDP(KK))
	PCP(KK)=PCONOF(PP(KK),TP(KK),TCP(KK))
	TETAWP(KK)=TETAP(KK)-WOBF(TETAP(KK))+WOBF(TCP)

732		CONTINUE
733			CONTINUE

C	------------------------------------------------------
C	CALCULO ALTURA DE LAS SUPERFICIES ISOBARICAS CADA 10 MB

	IS=INT(P(1)/10.)
	HP(IS)=ingrid(I,J,1,igeo)

		DO 823 KK=INT(P(1)/10.),NINT(P(maxnl)/10.),-1
	RM=(461.51/287.05)*((QP(KK-1)+QP(KK))/2.)+
     $     (1.-(QP(KK-1)+QP(KK))/2.)
	TVM=((TP(KK-1)+TP(KK))/2.+273.15)*RM
	HP(KK-1)=HP(KK)+(287.05/9.8)*TVM*LOG(FLOAT(KK)/FLOAT(KK-1))
823		CONTINUE


C	----------------------------------------------------------
C	------------------------------------------------------
C	TEMPERATURA POTENCIAL TERMOM HUMEDO PARTICULA SUPERFICIE
	TCS=TCONOF(ingrid(I,J,1,it),TD(1))
     	PCS=PCONOF(P(1),ingrid(I,J,1,it),TCS)
	TETAWS=TETA(1)-WOBF(TETA(1))+WOBF(TCS)
C	--------------------------------------------------------
C	-----------------------------------------------------
C	------------------------------------------------------	
C	CALCULO DE LA ENERGIA POTENCIAL CONVECTIVA
C	DISPONIBLE (CAPE) O CONVECTIVE BOUYANCY ENERGY (CBE).

	KK1=NINT(P(1)/10.)
	KK2=NINT(PCS/10.)
		DO 567 KK=KK1,KK2,-1
	TETAPART(KK)=((273.16+ingrid(I,J,1,it))*(1000./P(1))**.28571482)-273.16
567		CONTINUE


	KK2=KK2-1		
		DO WHILE (KK2.GE.10)
	PPART=FLOAT(KK2)*10.
	TPART=SATLFT(TETAWS,PPART)
	TETAPART(KK2)=((273.16+TPART)*(1000./PPART)**.28571482)-273.16
	KK2=KK2-1
		END DO


	KK=NINT(P(1)/10.)-1
	ENERG_POS1=0.
		DO WHILE(TETAP(KK).GE.TETAPART(KK))
	TETAPARTM=(TETAPART(KK)+TETAPART(KK+1))/2.
	TETAPM=(TETAP(KK)+TETAP(KK+1))/2.
        HPM=(HP(KK)-HP(KK+1))
	XINTEGRAL=((TETAPM-TETAPARTM)/(TETAPM+273.16))*HPM
	IF(XINTEGRAL.LT.0.) XINTEGRAL=0.
	ENERG_POS1=ENERG_POS1+XINTEGRAL
	IF(KK.EQ.10) GOTO 679
	KK=KK-1	
		END DO
679	ENERG_POS1=ENERG_POS1*9.8

	ENERG_NEG1=0.
	IF(KK.EQ.10) GOTO 673

		DO WHILE(TETAP(KK).LE.TETAPART(KK))	
	TETAPARTM=(TETAPART(KK)+TETAPART(KK+1))/2.
	TETAPM=(TETAP(KK)+TETAP(KK+1))/2.
        HPM=(HP(KK)-HP(KK+1))
	XINTEGRAL=(TETAPARTM-TETAPM)/(TETAPM+273.16)*HPM
	IF(XINTEGRAL.LT.0.) XINTEGRAL=0.
	ENERG_NEG1=ENERG_NEG1+XINTEGRAL
	IF(KK.EQ.10) GOTO 681
	KK=KK-1	
		END DO


681	ENERG_NEG1=ENERG_NEG1*9.8
	IWMAX1=NINT(SQRT(2.*ENERG_NEG1*9.8))
	PCN=P(1)
	HCN=HP(NINT(P(1)/10.))


673		CONTINUE


	IF(ENERG_NEG1.gt.0.) goto 1000


	IF(ENERG_NEG1.LE.0.) THEN

	NN1=NINT(P(1)/10.)-1
		DO 568 NN=NN1,50,-1
	KK1=NN
	KK2=NINT(PCP(NN)/10.)
		DO 368 KK=KK1,KK2,-1
	TETAPART(KK)=((273.16+TP(NN))*(1000./PP(NN))**.28571482)-273.16
368		CONTINUE


	KK2=KK2-1		
		DO WHILE (KK2.GE.10)
	PPART=FLOAT(KK2)*10.
	TPART=SATLFT(TETAWP(NN),PPART)
	TETAPART(KK2)=((273.16+TPART)*(1000./PPART)**.28571482)-273.16
	KK2=KK2-1
		END DO

	KK=NN-1

	ENERG_POS1=0.
		DO WHILE(TETAP(KK).GE.TETAPART(KK))
	TETAPARTM=(TETAPART(KK)+TETAPART(KK+1))/2.
	TETAPM=(TETAP(KK)+TETAP(KK+1))/2.
        HPM=(HP(KK)-HP(KK+1))
	XINTEGRAL=((TETAPM-TETAPARTM)/(TETAPM+273.16))*HPM
	IF(XINTEGRAL.LT.0.) XINTEGRAL=0.
	ENERG_POS1=ENERG_POS1+XINTEGRAL
C	PRINT *,PP(KK),TETAP(KK),TETAPM,TETAPART(KK),TETAPARTM
	IF(KK.EQ.10) GOTO 779
	KK=KK-1	
		END DO
779	ENERG_POS1=ENERG_POS1*9.8

	ENERG_NEG1=0.
	IF(KK.EQ.10) GOTO 781

		DO WHILE(TETAP(KK).LE.TETAPART(KK))	
	TETAPARTM=(TETAPART(KK)+TETAPART(KK+1))/2.
	TETAPM=(TETAP(KK)+TETAP(KK+1))/2.
        HPM=(HP(KK)-HP(KK+1))
	XINTEGRAL=(TETAPARTM-TETAPM)/(TETAPM+273.16)*HPM
	IF(XINTEGRAL.LT.0.) XINTEGRAL=0.
	ENERG_NEG1=ENERG_NEG1+XINTEGRAL
C	PRINT *,PP(KK),TETAP(KK),TETAPM,TETAPART(KK),TETAPARTM
	IF(KK.EQ.10) GOTO 781
	KK=KK-1	
		END DO


781	ENERG_NEG1=ENERG_NEG1*9.8
	IWMAX1=NINT(SQRT(2.*ENERG_NEG1*9.8))

	PCN=PP(NN)
	HCN=HP(NN)

	IF(ENERG_NEG1.EQ.0.) THEN
		HCN=0.
		PCN=P(1)
	ENDIF
	
	IF(ENERG_NEG1.GT.0.) GOTO 1000


568	CONTINUE
		END IF


1000		CONTINUE

	do il=1,outnl
		outgrid(i,j,il)=ENERG_POS1
	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 DETERMINAR LA PRESION DEL NIVEL DE CONDENSACION
C 	POR ELEVACION. ENTRADAS: PRESION INICIAL, TEMPERATURA INICIAL
C	Y TEMPERATURA DEL NIVEL DE CONDENSACION (OBTENIDA CON TCONOF)
C	PRES EN HPA. TEMP Y TCONOF EN C o K.


	FUNCTION PCONOF(PRES,TEMP,TCONOF)

	IF (TEMP-100.)4,4,5
4	TEMP=TEMP+273.16
	TCONOF=TCONOF+273.16

5	TETA=TEMP*(1000./PRES)**.28571428
	PCONOF=1000.*(TCONOF/TETA)**(1/.28571428)

	TEMP=TEMP-273.16
	TCONOF=TCONOF-273.16

	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 ===========================================================================
