 	PROGRAM MIXDRAFT
C
C   ***  This program accepts input from the file <<sounding>>    ***
C   ***   and calculates the properties of downdrafts formed by   ***
C   ***     mixing reversibly lifted near-surface air with        ***
C   ***        environmental air. Output is contained in          ***
C   ***               the file <<mixdraft.out>>                   ***
C
	PARAMETER (NA=60)
	REAL LV(NA),LV0
	REAL P(NA),TP(NA),TVP(NA),LW(NA)
	REAL T(NA),R(NA),THETA(NA),EV(NA),ES(NA)
	REAL RS(NA),RDCP(NA)
	OPEN(UNIT=11,FILE='sounding',STATUS='OLD')
	OPEN(UNIT=12,FILE='mixdraft.out',STATUS='NEW')
C
C   ***   Assign values of thermodynamic constants     ***
C
        CPD=1005.7
        RV=461.5
        RD=287.04
	EPS=RD/RV
        LV0=2.501E6
	CPV=1870.0
        CL=4190.0
        CPVMCL=CL-CPV
C
C   ***   Read in the sounding from the file <<sounding>>          ***
c
	READ(11,10)N
   10	FORMAT(4X,I3,//)
	DO 20 I=1,N
	 READ(11,15)P(I),T(I),R(I)
   15	 FORMAT(5X,F9.3,8X,F9.3,12X,F9.3)
	 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)))
	 RS(I)=EPS*ES(I)/(P(I)-ES(I))
	 RS(I)=MAX(RS(I),0.0)
	 R(I)=MIN(R(I),RS(I))
	 LV(I)=LV0-CPVMCL*T(I)
	 T(I)=T(I)+273.15
	 RDCP(I)=(RD+R(I)*RV)/(CPD+R(I)*CPV)
	 THETA(I)=T(I)*(1000.0/P(I))**RDCP(I)
   20	CONTINUE
	CLOSE(UNIT=11,DISPOSE='SAVE')
C
C   ***  Calculate lifted condensation pressure   ***
C
        TC=2840.0/(3.5*LOG(T(1))-LOG(EV(1))-4.805)+55.0
        PLCL=P(1)*((TC/T(1))**3.48)
C
C   ***  Calculate first level above the LCL (=IB)  ***
C
        IB=N-1
        DO 80 I=2,N-1
         IF(P(I).LT.PLCL)THEN
          IB=MIN(IB,I)
         END IF
   80   CONTINUE
	EM=MAX(EV(1),1.0E-6) 
C
C   ***  Calculate entropy of air at the surface   ***
C
	S=(CPD+R(1)*CL)*LOG(T(1))-RD*LOG(P(1)-EV(1))+
     1    LV(1)*R(1)/T(1)-R(1)*RV*LOG(EM/ES(1))
C
C   ***  Find properties of lifted surface air below cloud base  ***
C
	DO 100 I=1,IB-1
	 TP(I)=T(1)*(P(I)/P(1))**(RDCP(1))
	 TVP(I)=TP(I)*(1.+R(1)/EPS)/(1.+R(1))
	 LW(I)=0.0
  100	CONTINUE
C
C   ***  Find properties of lifted surface air above cloud base   ***
C
	DO 200 I=IB,N
C
C   ***   Calculate approximate rate of change of entropy with   ***
C   ***           temperature at constant pressure               ***
C
	 SL=(CPD+R(1)*CL+LV(I)*LV(I)*RS(I)/
     1     (RV*T(I)*T(I)))/T(I)
C
C   ***   Iterate to find lifted parcel temperature and mixing ratio  ***
C
	 TG=T(I)
	 RG=RS(I)
	 TOLER=0.1
	 DO 150 K=1,4 
	  EM=RG*P(I)/(EPS+RG)
	  ALV=LV0-CPVMCL*(TG-273.15)
	  SG=(CPD+R(1)*CL)*LOG(TG)-RD*LOG(P(I)-EM)+
     1      ALV*RG/TG
	  IF(ABS(S-SG).LT.TOLER)GOTO 160
	  TG=TG+(S-SG)/SL  
	  TC=TG-273.15
	  ENEW=6.112*EXP(17.67*TC/(243.5+TC))
	  RG=EPS*ENEW/(P(I)-ENEW)           
  150	 CONTINUE
  160	 CONTINUE
	 TP(I)=TG
	 TVP(I)=TG*(1.+RG/EPS)/(1.+R(1))
	 LW(I)=R(1)-RG
	 LW(I)=MAX(LW(I),0.0)
  200	CONTINUE
C
C   ***  Find the updraft velocity, the mixing fraction of   ***
C   ***   environmental air that creates a just-saturated    ***
C   ***     mixture, the maximum downdraft speed, and        ***
C   ***         the maximum displacement of the              ***
C   ***               downward moving air                    ***
C
	WRITE(12,210)
  210	FORMAT(1X,'Pressure',2X,' Mixing ',2X,' Max. downward   ',
     1   2X, 'Max. downdraft',2X,' Max. updraft ')
	WRITE(12,212)
  212	FORMAT(1X,'  (mb)  ',2X,'fraction',2X,'displacement (mb)',
     1   2X, 'velocity (m/s)',2X,'velocity (m/s)')
	WRITE(12,214)
  214	FORMAT(1X,'--------',2X,'--------',2X,'-----------------',
     1    2X, '-------------',2X,'--------------')
	THEVB=THETA(1)*(1.+1.608*R(1))/(1.+R(1))
	CAPE=0.0
	  DO 300 I=IB,N
	  THEVI=THETA(I)*(1.+1.608*R(I))/(1.+R(I))
	  TVI=THEVI*T(I)/THETA(I)
	  TVIM=T(I-1)*(1.+R(I-1)/EPS)/(1.+R(I-1))
	  TVM=0.5*(TVP(I)-TVI+TVP(I-1)-TVIM)
	  PM=0.5*(P(I)+P(I-1))
	  CAPE=CAPE+RD*(TVP(I)-TVI)*(P(I-1)-P(I))/PM  
	  CAPE=MAX(CAPE,0.0)
	  WUP=SQRT(CAPE)                      
	  SIGM=(RS(I)-R(I))/(R(1)-R(I)+LV(I)*RS(I)*(
     1     THEVI-THEVB)/(RV*THEVI*T(I)))
	  SIGM=MAX(SIGM,0.0)
	  SIGM=MIN(SIGM,1.0)
	  DELP=(PLCL-P(I))*SIGM
	  PNEW=P(I)+DELP
	  WD2=CPD*(THEVI-THEVB)*SIGM*((0.001*PNEW)**RDCP(IB)-
     1     T(I)/THETA(I))
	  WD2=MAX(WD2,0.0)
	  WD=-SQRT(WD2)
C
C   ***    Write to ouput file   ***
C
	  WRITE(12,250)P(I),SIGM,DELP,WD,WUP
  250	  FORMAT(1X,F6.1,5X,F5.2,8X,F5.1,
     1     13X,F6.1,10X,F6.1)
  300	  CONTINUE
	CLOSE(UNIT=12,DISPOSE='SAVE')
	STOP
	END	
