       SUBROUTINE PCMIN(SST,PSL,P,T,R,NA,N,PMIN,VMAX,IFL)
C
C   ***   This subroutine calculates the mimimum central pressure  ***
C   ***    achievable in tropical cyclones, given a sounding       ***
C   ***             and a sea surface temperature.                 ***
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.
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           IFL is a flag: A value of 1 means OK; A value of 0
C              indicates no convergence (hypercane)
C
C-----------------------------------------------------------------------------
	PARAMETER(NB=80)
	REAL T(NA), P(NA), R(NA)
	REAL EV(NB), ES(NB), TVRDIF(NB)
	REAL TLR(NB), TLVR(NB),  LW(NB)
	REAL NAR
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   ***   Adjustable constant        ***
C
	CKCD=1.0
C
C   ***   Normalize certain quantities   ***
C
	SSTK=SST+273.15
	IFL=1
	DO 40 I=1,N
	 R(I)=R(I)*0.001
	 EV(I)=R(I)*P(I)/(EPS+R(I))
	 ES(I)=6.112*EXP(17.67*T(I)/(243.5+T(I)))
	 T(I)=T(I)+273.15
C
C   ***   Default Values   ***
C
	 TLR(I)=T(I)
	 TLVR(I)=T(I)*(1.+R(I)/EPS)/(1.+R(I))
	 LW(I)=0.0
	 TVRDIF(I)=0.0
   40	CONTINUE
	ES0=6.112*EXP(17.67*SST/(243.5+SST))
	PMIN=PSL
	PC=950.0
	NP=0
C
C   ***   Begin iteration to find mimimum pressure   ***
C
  100	CONTINUE
C
C   ***  Define various conserved parcel quantities: reversible   ***
C   ***        entropy, S, and total water content, RS0           ***
C
	RS0=EPS*ES0/(PC-ES0)
	ALV=ALV0-CPVMCL*SST
	S=(CPD+RS0*CL)*LOG(SSTK)-RD*LOG(PC-ES0)+
     1    ALV*RS0/SSTK
	NCMAX=0
C
C   ***  Begin updraft loop   ***
C
	DO 200 J=2,N
C
C    ***   Don't bother lifting parcel above 40 mb    ***
C
C	 IF(P(J).LT.39.0)GOTO 200
	 IF(P(J).GT.PC)THEN
	  TG=(T(1)*(P(J)-PC)+SSTK*(P(1)-P(J)))/(P(1)-PC)
	  RG=(R(1)*(P(J)-PC)+RS0*(P(1)-P(J)))/(P(1)-PC)
	  TLR(J)=TG
	  TLVR(J)=TG*(1.+RG/EPS)/(1.+RG)
	  LW(J)=0.0
	  TVRDIF(J)=TLVR(J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 ELSE
	 RS=EPS*ES(J)/(P(J)-ES(J))
C
C   ***  Iteratively calculate lifted parcel temperature and mixing   ***
C   ***                ratio for reversible ascent                    ***
C
	 TG=T(J)
	 RG=RS
	 NC=0
C
C   ***  Begin iteration to find the updraft temperature   ***
C
  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+RS0*CL+ALV*ALV*RG/(RV*TG*TG))/TG
	  EM=RG*P(J)/(EPS+RG)
	  SG=(CPD+RS0*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
	    PMIN=400.0
	    IFL=0
	    GOTO 900
	   END IF
	   RG=EPS*ENEW/(P(J)-ENEW)           
	   GOTO 120
	  END IF
  150	 CONTINUE
	 NCMAX=MAX(NC,NCMAX)
	 TLR(J)=TG
	 TLVR(J)=TG*(1.+RG/EPS)/(1.+RS0)
	 LW(J)=RS0-RG
	 LW(J)=MAX(0.0,LW(J))
	 TVRDIF(J)=TLVR(J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 END IF
  200	CONTINUE
C
C  ***  Begin loop to find CAPE, PA, and NA from reversible ascent ***
C
	 CAPER=0.0
	 NAR=0.0
	 PAR=0.0
C
C   ***  Find maximum level of positive buoyancy    ***
C
	 INBR=1
	 DO 550 J=N,1,-1
	  IF(TVRDIF(J).GT.0.0)INBR=MAX(INBR,J)
  550	 CONTINUE
	 IF(INBR.LT.N)THEN
	  DO 555 J=INBR+1,N
	   TVRDIF(J)=0.0
  555	  CONTINUE
	 END IF
C
C   ***  Do updraft loop        ***
C
	 IF(INBR.GT.1)THEN
	  DO 600 J=2,INBR
	   TVM=0.5*(TVRDIF(J)+TVRDIF(J-1))
	   PM=0.5*(P(J)+P(J-1))
	   IF(TVM.LE.0.0)THEN
	    NAR=NAR-RD*TVM*(P(J-1)-P(J))/PM
	   ELSE
	    PAR=PAR+RD*TVM*(P(J-1)-P(J))/PM
	   END IF
  600	  CONTINUE
	  CAPER=PAR-NAR
	  CAPERM=MAX(CAPER,0.0)
	 END IF
  800	CONTINUE
C
C   ***  Define various conserved parcel quantities: reversible   ***
C   ***        entropy, S, and total water content, RS0           ***
C
        RS0=R(1)*(PSL-EV(1))/(PC-EV(1))
        ALV=ALV0-CPVMCL*(T(1)-273.15)
	EV(1)=MAX(EV(1),0.01)
        S=(CPD+RS0*CL)*LOG(T(1))-RD*LOG(PC-EV(1))+
     1   ALV*RS0/T(1)-RS0*RV*LOG(EV(1)/ES(1))
        NCMAX=0
C
C   ***  Begin updraft loop   ***
C
        DO 1200 J=2,N
C
C    ***   Don't bother lifting parcel above 40 mb    ***
C
C        IF(P(J).LT.39.0)GOTO 1200
         IF(P(J).GT.PC)THEN
          TG=T(1)
          RG=(R(1)*(P(J)-PC)+RS0*(P(1)-P(J)))/(P(1)-PC)
          TLR(J)=TG
          TLVR(J)=TG*(1.+RG/EPS)/(1.+RG)
          LW(J)=0.0
          TVRDIF(J)=TLVR(J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
         ELSE
         RS=EPS*ES(J)/(P(J)-ES(J))
C
C   ***  Iteratively calculate lifted parcel temperature and mixing   ***
C   ***                ratio for reversible ascent                    ***
C
         TG=T(J)
         RG=RS
         NC=0
C
C   ***  Begin iteration to find the updraft temperature   ***
C
 1120    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+RS0*CL+ALV*ALV*RG/(RV*TG*TG))/TG
          EM=RG*P(J)/(EPS+RG)
          SG=(CPD+RS0*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
            PMIN=400.0
            IFL=0
            GOTO 900
           END IF
           RG=EPS*ENEW/(P(J)-ENEW)
           GOTO 1120
          END IF
 1150    CONTINUE
         NCMAX=MAX(NC,NCMAX)
         TLR(J)=TG
         TLVR(J)=TG*(1.+RG/EPS)/(1.+RS0)
         LW(J)=RS0-RG
         LW(J)=MAX(0.0,LW(J))
         TVRDIF(J)=TLVR(J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
         END IF
 1200   CONTINUE
C
C  ***  Begin loop to find CAPE, PA, and NA from reversible ascent ***
C
         CAPER=0.0
         NAR=0.0
         PAR=0.0
C
C   ***  Find maximum level of positive buoyancy    ***
C
         INBR=1
         DO 1550 J=N,1,-1
          IF(TVRDIF(J).GT.0.0)INBR=MAX(INBR,J)
 1550    CONTINUE
         IF(INBR.LT.N)THEN
          DO 1555 J=INBR+1,N
           TVRDIF(J)=0.0
 1555     CONTINUE
         END IF
C
C   ***  Do updraft loop        ***
C
         IF(INBR.GT.1)THEN
          DO 1600 J=2,INBR
           TVM=0.5*(TVRDIF(J)+TVRDIF(J-1))
           PM=0.5*(P(J)+P(J-1))
           IF(TVM.LE.0.0)THEN
            NAR=NAR-RD*TVM*(P(J-1)-P(J))/PM
           ELSE
            PAR=PAR+RD*TVM*(P(J-1)-P(J))/PM
           END IF
 1600     CONTINUE
          CAPERA=PAR-NAR
          CAPERA=MAX(CAPERA,0.0)
         END IF
 1800   CONTINUE
C
C  ***  Initial estimate of minimum pressure   ***
C
	TVAV=0.5*(TLVR(1)+SSTK*(1.+RS0/EPS)/(1.+RS0))
	CAT=CAPERA+0.5*CKCD*(CAPERM-CAPERA)
	PNEW=PSL*EXP(-CAT/(RD*TVAV))
C
C   ***  Test for convergence   ***
C
	IF(ABS(PNEW-PC).GT.0.2)THEN
	 PC=PNEW
	 NP=NP+1
	 IF(NP.GT.1000.OR.PC.LT.400.0)THEN
	  PMIN=400.0
	  IFL=0
	  GOTO 900
	 END IF
	 GOTO 100
	ELSE
	 CAT=CAPERA+CKCD*(CAPERM-CAPERA)
	 CAT=MAX(CAT,0.0)
	 PMIN=PSL*EXP(-CAT/(RD*TVAV))
	END IF
  900	CONTINUE
        CDIF=CKCD*MAX((CAPERM-CAPERA),0.0)
        VMAX=0.9*SQRT(CDIF)
	DO 910 I=1,N
	 R(I)=R(I)*1000.0
	 T(I)=T(I)-273.15
  910	CONTINUE
	RETURN
	END
	 
	
	
	 
	
	
