       SUBROUTINE INDEX(SST,PSL,P,T,R,U,V,ALAT,NA,N,PMIN,VMAX,GINDX,IFL)
C
C   ***   This subroutine calculates the mimimum central pressure  ***
C   ***            and maximum 10 m 1 minute wind speed            ***
C   ***    achievable in tropical cyclones, given a sounding       ***
C   ***             and a sea surface temperature.                 ***
C   ***           A genesis index is also calculated               ***
C
C  INPUT:   SST: Sea surface temperature in C
C
C           PSL: Sea level pressure (mb)
C
C           P,T,R: One-dimensional arrays of dimension NA
C             containing pressure (mb), temperature (C),
C             and mixing ratio (g/kg). The arrays MUST be
C             arranged so that the lowest index corresponds
C             to the lowest model level, with increasing index
C             corresponding to decreasing pressure. The sounding
C             should extend to at least the tropopause and 
C             preferably to the lower stratosphere.
C
C           U,V: One dimensional arrays of dimension NA containing
C             zonal and meridional wind components (m/s), ordered
C             the same way as P, T and R. Any missing winds should
C             be filled in by interpolation.
C
C           ALAT: The latitude of the column, in degrees. 
C
C           NA: The dimension of P,T and R
C
C           N:  The actual number of points in the sounding
C                (N is less than or equal to NA)
C
C  OUTPUT:  PMIN is the minimum central pressure, in mb
C
C           VMAX is the maximum surface wind speed, in m/s
C                  (reduced to reflect surface drag)
C
C           GINDX is a genesis index, varying from 0 to 100,
C              with higher values indicating higher probabilities
C              of genesis.
C
C           IFL is a flag: A value of 1 means OK; a value of 0
C              indicates no convergence (hypercane); a value of 2
C              means that the CAPE routine failed.
C
C-----------------------------------------------------------------------------
	REAL T(NA), P(NA), R(NA), U(NA), V(NA)
C
C   ***   Adjustable constant: Ratio of C_k to C_D    ***
C
	CKCD=1.0
C
C   ***   Adjustable constant for buoyancy of displaced parcels:  ***
C   ***    0=Reversible ascent;  1=Pseudo-adiabatic ascent        ***
C
        SIG=0.0
C
C   ***  Adjustable switch: if IDISS = 0, no dissipative heating is   ***
C   ***     allowed; otherwise, it is                                 ***
C
	IDISS=1
C
C   ***  Adjustable constants used in genesis index    ***
C
	C1=0.2
	C2=8.0
C
C   ***   Normalize certain quantities   ***
C
	SSTK=SST+273.15
	ES0=6.112*EXP(17.67*SST/(243.5+SST))
	DO 40 I=1,N
	 R(I)=R(I)*0.001
	 T(I)=T(I)+273.15
   40	CONTINUE
C
C   ***   Default values   ***
C
        VMAX=0.0
	PMIN=PSL
	GINDX=0.0
	F1=0.0
	F2=0.0
	IFL=1
C
	NP=0
	PM=950.0
C
C   ***   Find environmental CAPE *** 
C
      TP=T(1)
      RP=R(1)
      PP=P(1)
      CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEA,ANA,TOA,INB,IFLAG)
      IF(IFLAG.NE.1)IFL=2
C
C   ***   Begin iteration to find mimimum pressure   ***
C
  100 CONTINUE
C
C   ***  Find CAPE at radius of maximum winds   ***
C
      TP=T(1)
      PP=PM
      RP=0.622*R(1)*PSL/(PM*(0.622+R(1))-R(1)*PSL)
      CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEM,AN,TOM,INBA,IFLAG) 
      IF(IFLAG.NE.1)IFL=2
C
C  ***  Find saturation CAPE at radius of maximum winds   ***
C
      TP=SSTK
      PP=PM
      RP=0.622*ES0/(PM-ES0)
      CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEMS,ANS,TOMS,INB,IFLAG)
      IF(IFLAG.NE.1)IFL=2
      RAT=SSTK/TOM
      IF(IDISS.EQ.0)RAT=1.0
C
C  ***  Initial estimate of minimum pressure   ***
C
      RS0=RP
      TV1=T(1)*(1.+R(1)/0.622)/(1.+R(1))
	TVAV=0.5*(TV1+SSTK*(1.+RS0/0.622)/(1.+RS0))
	CAT=CAPEM-CAPEA+0.5*CKCD*RAT*(CAPEMS-CAPEM)
	PNEW=PSL*EXP(-CAT/(287.04*TVAV))
C
C   ***  Test for convergence   ***
C
	IF(ABS(PNEW-PM).GT.0.2)THEN
	 PM=PNEW
	 NP=NP+1
	 IF(NP.GT.1000.OR.PM.LT.400.0)THEN
	  PMIN=400.0
	  IFL=0
	  GOTO 900
	 END IF
	 GOTO 100
	ELSE
	 CAT=CAPEM-CAPEA+CKCD*RAT*(CAPEMS-CAPEM)
	 PMIN=PSL*EXP(-CAT/(287.04*TVAV))
	END IF
  900	CONTINUE
	FAC=MAX(0.0,(CAPEMS-CAPEM))
	VMAX=0.8*SQRT(CKCD*RAT*FAC)
C
C   ***   Calculate genesis index     ***
C
	INBH=INBA/2
	IF(VMAX.LT.5.0.OR.INBH.LT.3)RETURN
C
C   ***   Calculate storm layer average wind shear   ***
C
	UBAR1=0.0
	VBAR1=0.0
	UBAR2=0.0
	VBAR2=0.0
	DO 1010 I=1,INBH-1
	 UBAR1=UBAR1+0.5*(U(I)+U(I+1))*(P(I)-P(I+1))
	 VBAR1=VBAR1+0.5*(V(I)+V(I+1))*(P(I)-P(I+1))
 1010	CONTINUE
	DO 1020 I=INBH,INBA-1
	 UBAR2=UBAR2+0.5*(U(I)+U(I+1))*(P(I)-P(I+1))
	 VBAR2=VBAR2+0.5*(V(I)+V(I+1))*(P(I)-P(I+1))
 1020	CONTINUE
	UBAR1=UBAR1/(P(1)-P(INBH))
	VBAR1=VBAR1/(P(1)-P(INBH))
	UBAR2=UBAR2/(P(INBH)-P(INBA))
	VBAR2=VBAR2/(P(INBH)-P(INBA))
C
C   ***  Average low level thetae deficit   ***
C
C
C   ***   Assign values of thermodynamic constants     ***
C
        CPD=1005.7
        CPV=1870.0
        CL=4190.0
        CPVMCL=2320.0
        RV=461.5
        RD=287.04
        EPS=RD/RV
        ALV0=2.501E6
C
        TPC=T(1)-273.15
        ESP=6.112*EXP(17.67*TPC/(243.5+TPC))
        EVP=R(1)*P(1)/(EPS+R(1))
        RH=EVP/ESP
	RH=MAX(RH,0.01)
        ALV=ALV0-CPVMCL*TPC
        SS=(CPD+R(1)*CL)*LOG(T(1))-RD*LOG(P(1)-EVP)+
     1   ALV*R(1)/T(1)-R(1)*RV*LOG(RH)	
	SBAR=0.0
	PBAR=0.0
	DO 1030 I=2,INBA
	 IF(P(I).GT.500.0.AND.P(I).LE.920.0)THEN
          TPC=T(I)-273.15
          ESP=6.112*EXP(17.67*TPC/(243.5+TPC))
          EVP=R(I)*P(I)/(EPS+R(I))
          RH=EVP/ESP
	  RH=MAX(RH,0.01)
          ALV=ALV0-CPVMCL*TPC
          SBAR=SBAR+((CPD+R(1)*CL)*LOG(T(I))-RD*LOG(P(I)-EVP)+
     I     ALV*R(I)/T(I)-R(I)*RV*LOG(RH))*(P(I-1)-P(I+1))	  
	  PBAR=PBAR+P(I-1)-P(I+1)
	 END IF
 1030	CONTINUE
	IF(SBAR.EQ.0.0)RETURN
	SBAR=SBAR/PBAR
C
C   ***  Calculate Index   ***
C
	F1=((SSTK-TOM)*(SS-SBAR)+ANA)/CPD
	F1=MAX(F1,0.0)
	F2=(1.0+SQRT((VBAR2-VBAR1)**2+(UBAR2-UBAR1)**2))/
     1    (1.0E-4+VMAX)
	F3=1.-C1*SQRT(F1)*(1.+C2*F2)
	F3=MAX(F3,0.0)
	F4=MIN(1.0,(ALAT/5.0)**2)
	GINDX=100.0*F3*F4
C
C   ***  Renormalize sounding arrays   ***
C
        DO 910 I=1,N
         R(I)=R(I)*1000.0
         T(I)=T(I)-273.15
  910   CONTINUE
C
	RETURN
	END
C        
      SUBROUTINE CAPE(TP,RP,PP,T,R,P,ND,N,SIG,CAPED,NA,TO,INB,IFLAG)
C
C     This subroutine calculates the CAPE of a parcel with pressure PP (mb), 
C       temperature TP (K) and mixing ratio RP (gm/gm) given a sounding
C       of temperature (T in K) and mixing ratio (R in gm/gm) as a function
C       of pressure (P in mb). ND is the dimension of the arrays T,R and P,
C       while N is the actual number of points in the sounding. CAPED is
C       the calculated value of CAPE, NA is the "negative area" in the
C       sounding, and TO is the temperature at the
C       level of neutral buoyancy.  IFLAG is a flag
C       integer. If IFLAG = 1, routine is successful; if it is 0, routine did
C       not run owing to improper sounding (e.g.no water vapor at parcel level).
C       IFLAG=2 indicates that routine did not converge.                 
C
      REAL T(ND),R(ND),P(ND),TVRDIF(100)   
      REAL NA
C
C   ***   Default values   ***
C      
      CAPED=0.0
      NA=0.0
      TO=T(1)
      INB=2
      IFLAG=1
C
C   ***   Check that sounding is suitable    ***
C
      IF(RP.LT.1.0E-6.OR.TP.LT.200.0)THEN
       IFLAG=0
       RETURN
      END IF            
C
C   ***   Assign values of thermodynamic constants     ***
C
      CPD=1005.7
      CPV=1870.0
      CL=4190.0
      CPVMCL=2320.0
      RV=461.5
      RD=287.04
      EPS=RD/RV
      ALV0=2.501E6
C
C   ***  Define various parcel quantities, including reversible   ***
C   ***                       entropy, S.                         ***
C                           
      TPC=TP-273.15
      ESP=6.112*EXP(17.67*TPC/(243.5+TPC))
      EVP=RP*PP/(EPS+RP)
      RH=EVP/ESP
      ALV=ALV0-CPVMCL*TPC
      S=(CPD+RP*CL)*LOG(TP)-RD*LOG(PP-EVP)+
     1   ALV*RP/TP-RP*RV*LOG(RH)            
C
C   ***  Find lifted condensation pressure, PLCL   ***
C     
	CHI=TP/(1669.0-122.0*RH-TP)
	PLCL=PP*(RH**CHI)
C
C   ***  Begin updraft loop   ***
C
	NCMAX=0
	DO 200 J=2,N
C
C    ***   Don't bother lifting parcel above 60 mb    ***
C
	 IF(P(J).LT.59.0)GOTO 200
C
C    ***  Parcel quantities below lifted condensation level   ***
C	 
	 IF(P(J).GE.PLCL)THEN
	  TG=TP*(P(J)/PP)**(RD/CPD)
	  RG=RP
C
C   ***   Calculate buoyancy   ***
C  
	  TLVR=TG*(1.+RG/EPS)/(1.+RG)
	  TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 ELSE
C
C   ***  Parcel quantities above lifted condensation level  ***
C	 
	  TG=T(J)          
	  TJC=T(J)-273.15 
	  ES=6.112*EXP(17.67*TJC/(243.5+TJC)) 
	  RG=EPS*ES/(P(J)-ES)
C
C   ***  Iteratively calculate lifted parcel temperature and mixing   ***
C   ***                ratio for reversible ascent                    ***
C
	  NC=0
  120	  CONTINUE
	  NC=NC+1
C
C   ***  Calculate estimates of the rates of change of the entropy    ***
C   ***           with temperature at constant pressure               ***
C  
	  ALV=ALV0-CPVMCL*(TG-273.15)
	  SL=(CPD+RP*CL+ALV*ALV*RG/(RV*TG*TG))/TG
	  EM=RG*P(J)/(EPS+RG)
	  SG=(CPD+RP*CL)*LOG(TG)-RD*LOG(P(J)-EM)+
     1      ALV*RG/TG
	  IF(NC.LT.3)THEN
	   AP=0.3
	  ELSE
	   AP=1.0
	  END IF
	  TGNEW=TG+AP*(S-SG)/SL  
C
C   ***   Test for convergence   ***
C
	  IF(ABS(TGNEW-TG).GT.0.01)THEN
	   TG=TGNEW
	   TC=TG-273.15
	   ENEW=6.112*EXP(17.67*TC/(243.5+TC))
C
C   ***   Bail out if things get out of hand   ***
C
	   IF(NC.GT.500.OR.ENEW.GT.(P(J)-1.0))THEN
            IFLAG=2
            RETURN
	   END IF
	   RG=EPS*ENEW/(P(J)-ENEW)           
	   GOTO 120
	  END IF
	  NCMAX=MAX(NC,NCMAX)
C
C   *** Calculate buoyancy   ***
C
          RMEAN=SIG*RG+(1.-SIG)*RP
	  TLVR=TG*(1.+RG/EPS)/(1.+RMEAN)
	  TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 END IF
  200	CONTINUE
C
C  ***  Begin loop to find NA, PA, and CAPE from reversible ascent ***
C
	NA=0.0
	PA=0.0
C
C   ***  Find maximum level of positive buoyancy, INB    ***
C
	INB=1
	DO 550 J=N,1,-1
	 IF(TVRDIF(J).GT.0.0)INB=MAX(INB,J)
  550	CONTINUE
C
C   ***  Find positive and negative areas and CAPE  ***
C
	IF(INB.GT.1)THEN
	 DO 600 J=2,INB
	  TVM=0.5*(TVRDIF(J)+TVRDIF(J-1))
	  PMA=0.5*(P(J)+P(J-1))
	  IF(TVM.LE.0.0)THEN
	   NA=NA-RD*TVM*(P(J-1)-P(J))/PMA
	  ELSE
	   PA=PA+RD*TVM*(P(J-1)-P(J))/PMA
	  END IF
  600	 CONTINUE
C
C   ***   Find residual positive area above INB, and TO  ***
C
       PAT=0.0
       TO=T(INB)
       IF(INB.LT.N)THEN
        PINB=(P(INB+1)*TVRDIF(INB)-P(INB)*TVRDIF(INB+1))/
     1   (TVRDIF(INB)-TVRDIF(INB+1))
        PAT=RD*TVRDIF(INB)*(P(INB)-PINB)/(P(INB)+PINB)
	TO=(T(INB)*(PINB-P(INB+1))+T(INB+1)*(P(INB)-PINB))/
     1    (P(INB)-P(INB+1))
       END IF
C
C   ***   Find CAPE  ***
C            
	 CAPED=PA+PAT-NA
	 CAPED=MAX(CAPED,0.0)
	END IF
C
	RETURN
	END
	 
	
	
	 
	
	
