	PROGRAM GATE
C
C Note: The parameter NA must be greater than or equal to the total number 
C       of model levels at which the state variables are defined.  This will
C       equal the number of tropospheric levels plus 9. NOTE THAT NA MUST BE
C       SET TO THE SAME VALUE IN SUBROUTINE CONVECT. NTIME is at least the
C       number of time steps that certain quanitities like precipitation are
C       output to the graphics arrays.
C
	PARAMETER (NA=19, NAP=NA+1, NTIME=5000, NCL=1)
	PARAMETER(NK=500,NI=19)
C 
C  Dimension various arrays
C
C  Temperatures:
C
	REAL T1(NA), T2(NA), T3(NA), TTEMP(NA)
	REAL TRAD(NA), TVP(NA)
C
C  Specific humidities:
C
	REAL R1(NA), R2(NA), R3(NA), RS(NA), RTEMP(NA) 
C
C  Pressures:
C
	REAL P(NA), PH(NAP)
C
C  Forcings:
C
	REAL FT(NA), FR(NA), FTNET(NA), FRNET(NA)
C
C  Other convective quantities used for diagnostics:
C
	REAL MP(NA), RP(NA), SIJ(NA,NA), WATER(NA)
	REAL FRSUB(NA)
	REAL MUP(NA), MDOWN(NA), ENT(NA), DET(NA), FRDET(NA)
C
C  Quantities from GATE Soundings
C
	REAL TIMER(NK),WT(NK,NI), RCOOL(NK,NI)
	REAL TR(NK,NI),RR(NK,NI),ADVT(NK,NI),ADVR(NK,NI)
	REAL U(NK,NI), V(NK,NI)
	REAL RSR(NK,NI), DUMMY(9,9,19)
C
C  Time-dependent quantities for graphics output:
C
	REAL RAIN_EVAP(NTIME,2), SST_TS(NTIME,2), GTIME(NTIME)
	REAL RHTIME(NTIME,NA),RHTIMR(NTIME,NA),RHDIF(NTIME,NA)
	REAL TPQTIME(NTIME,NA),MUTIME(NTIME,NA),MDTIME(NTIME,NA)
	REAL MPENTIME(NTIME,NA),TPMTIME(NTIME,NA)
	REAL OMEGA(NTIME,NA),PBUDGET(NTIME),EBUDGET(NTIME)
	REAL Q500(NTIME,2), Q850(NTIME,2),T200(NTIME,2)
	REAL T850(NTIME,2), ENTHALPY(NTIME), DOWNV(NTIME)
	REAL MSE(NTIME,2)
C
C  Mean quantities for output:
C
	REAL FTBAR(NA,5), FRBAR(NA,5), BUOYBAR(NA), RHBAR(NA,2)
	REAL HBAR(NA,2), MBAR(NA,3), MDET(NA,2), WATERBAR(NA,2)
	REAL TPBAR(NA,2)
C
C  Special graphics arrays:
C
	REAL PB1(NA), PB2(NA,2), PB3(NA,3), PB5(NA,5)
        CHARACTER*4 DSH2(2), DSH3(3), DSH5(5)
        CHARACTER*12 XLAB, YLAB, ZLAB
C
C  Latent heat at T=0 C
C
	REAL LV0
C
C   Pressure level data
C
	DATA P/ 991.25, 945.50, 893.79, 836.06, 772.82, 705.22, 635.05,
     *          564.48, 495.73, 430.71, 370.78, 316.72, 268.82, 226.98,
     *          190.82, 159.87, 133.55, 111.29,  92.56/
C
C  Common block with subroutine CONVECT:
C
	COMMON / CVT / TVP,MP,RP,SIJ,MUP,MDOWN,ENT,DET,FRDET,FRSUB,
     1     WATER 
C
C  Open input and output files:
C 
	OPEN(UNIT=8, FILE='params.in',STATUS='UNKNOWN')
	OPEN(UNIT=12,FILE='rhrms.out',ACCESS='APPEND',STATUS='UNKNOWN')
	OPEN(UNIT=9, FILE='GATE.T',STATUS='UNKNOWN')
	OPEN(UNIT=10, FILE='GATE.Q',STATUS='UNKNOWN')
	OPEN(UNIT=11, FILE='GATE.OMG',STATUS='UNKNOWN')
	OPEN(UNIT=13, FILE='gate.Htemp_sui',STATUS='UNKNOWN')
	OPEN(UNIT=14, FILE='gate.Hqv_sui', STATUS='UNKNOWN')
	OPEN(UNIT=15, FILE='gate.rad', STATUS='UNKNOWN')
	OPEN(UNIT=16, FILE='GATE.U', STATUS='UNKNOWN')
	OPEN(UNIT=17, FILE='GATE.V', STATUS='UNKNOWN')
        OPEN(UNIT=18,FILE='time.csv',STATUS='UNKNOWN')
        OPEN(UNIT=19,FILE='profile.csv',STATUS='UNKNOWN')
        OPEN(UNIT=20,FILE='contour.csv',STATUS='UNKNOWN')

C
C  Assign thermodynamic constants
C
	CPD=1005.7
	CPV=1870.0
	CL=4190.0
	RD=287.04
	RV=461.5
	LV0=2.501E6
	G=9.8
	ROWL=1000.0
	EPS=RD/RV
	CPVMCL=CL-CPV
	SIGD1=0.5
C
C  Asselin filter value
C
	AFILT=0.3
C
C  Headers for output files
C
      WRITE(18,3)
    3 FORMAT(1X,'Time (days),','Precipitation,','SST,','Surface T,',
     1 'Observed LHF,','Calculated LHF,','Observed Enthalpy,',
     2 'Predicted Enthalpy,','Calc. Enthalpy,','Radiative Cooling,',
     3 'Observed 500 mb q,','Predicted 500 mb q')
C
      WRITE(19,4)
    4 FORMAT(1X,'Pressure (mb),','Buoyancy,','Observed RH,',
     1 'Predicted RH,','Observed T,','Predicted T,',
     * 'Radiative Cooling Rate,','Upward Mass Flux,',
     2 'Penetrative Mass Flux,','Unsaturated downdraft M,',
     3 'Entrainment,','Detrainment')
C
      WRITE(20,5)
    5 FORMAT(1X,'Time (days),','Pressure (mb),',
     1 'Observed relative humidity,','Predicted relative humidity,',
     2 'Relative humidity difference,','Temperature Difference,',
     3 'Penetrative downdraft mass flux')
C

C
C  Read in the parameters of the model from the file params.in
C
	READ(8,10)
     1    ELCRIT,TLCRIT,SIGD,SIGS,SIGMD,
     2    WTRAIN,OMTSNOW,COEFFR,COEFFS,DTMAX,SLFRAC,ENTP,DEPTH,
     3    VS0,PFIX,AL1,BETA1,TAUR,TAUF,DT,AVTIME,PRINFREQ
C
   10   FORMAT(//,14(T32,F12.6,/),7(T32,F12.6,/),T32,F12.6)
C
        PSURF=1008.0     
	ENDTIME=160.*3.*3600.
C
C  Renormalize certain parameters
C
	DT=DT*60.0
	AVTIME=AVTIME*3600.0*24.0
	PRINFREQ=PRINFREQ*3600.0
C
C  Read in and average the GATE data   
C
	DO 39 K=1,161
        TIMER(K)=3.0*3600.*FLOAT(K-1) 
	DO 38 NMI=1,19
	SUM=0.0
	DO 19 I=1,9
	 READ(9,21)(DUMMY(I,J,NMI),J=1,9)
	 DO 20 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   20	 CONTINUE
   19   CONTINUE
   21	FORMAT(9F10.3)
        TR(K,NMI)=SUM/81.
 	IF(TR(K,NMI).GE.0.0)THEN
  	 ER=6.112*EXP(17.67*TR(K,NMI)/(243.5+TR(K,NMI)))
 	ELSE
  	 TRA=TR(K,NMI)+273.15
  	 ER=EXP(23.33086-6111.72784/TRA+0.15215*LOG(TRA))
 	END IF
        RSR(K,NMI)=EPS*ER/(P(NMI)-ER*(1.-EPS))
        TR(K,NMI)=TR(K,NMI)+273.15
        SUM=0.0
        DO 22 I=1,9
         READ(10,21)(DUMMY(I,J,NMI),J=1,9)
         DO 23 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   23    CONTINUE
   22   CONTINUE
        RR(K,NMI)=0.001*SUM/81.
        SUM=0.0
        DO 24 I=1,9
         READ(11,21)(DUMMY(I,J,NMI),J=1,9)
         DO 25 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   25    CONTINUE
   24   CONTINUE
        WT(K,NMI)=1.0E-5*SUM/81.
        SUM=0.0
        DO 26 I=1,9
         READ(13,21)(DUMMY(I,J,NMI),J=1,9)
         DO 27 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   27    CONTINUE
   26   CONTINUE
        ADVT(K,NMI)=SUM/(3600.*24.*81.)
        SUM=0.0
        DO 28 I=1,9
         READ(14,21)(DUMMY(I,J,NMI),J=1,9)
         DO 29 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   29    CONTINUE
   28   CONTINUE
        ADVR(K,NMI)=0.001*SUM/(3600.*24.*81.)
        SUM=0.0
        DO 30 I=1,9
         READ(15,32)(DUMMY(I,J,NMI),J=1,9)
         DO 31 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   31    CONTINUE
   30   CONTINUE
   32	FORMAT(9F6.2)
        RCOOL(K,NMI)=-SUM/(3600.*24.*81.)
	IF(NMI.EQ.19)RCOOL(K,NMI)=0.0
        SUM=0.0
        DO 33 I=1,9
         READ(16,21)(DUMMY(I,J,NMI),J=1,9)
         DO 34 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   34    CONTINUE
   33   CONTINUE
        U(K,NMI)=SUM/81.
        SUM=0.0
        DO 35 I=1,9
         READ(17,21)(DUMMY(I,J,NMI),J=1,9)
         DO 36 J=1,9
          SUM=SUM+DUMMY(I,J,NMI)
   36    CONTINUE
   35   CONTINUE
        V(K,NMI)=SUM/81.
   38	CONTINUE
   39	CONTINUE
	K=161
	NP=19
	N=18
	DO 40 I=1,NP
	 T1(I)=TR(1,I)
	 R1(I)=RR(1,I)
	 R1(I)=MAX(R1(I),2.0E-6)
	 T2(I)=T1(I)
	 R2(I)=R1(I)
	 IF(I.GT.1)THEN
	  PH(I)=0.5*(P(I)+P(I-1))  
	 END IF
   40   CONTINUE 
c	PH(1)=2.*PH(2)-PH(3)
	PH(1)=PSURF
	PH(NP+1)=2.*PH(NP)-PH(NP-1)
	PH(NP+1)=MAX(0.1,PH(NP+1))
C
C  Initialize certain quantities and arrays
C     
	ICB=2                                     
	FRSUM=0.0
	FHSUM=0.0
	PRECBAR=0.0
	FTS=0.0 
	CBMF=0.0
	NCOUNT=0
C
	DO 60 I=1,NP
	 FTNET(I)=0.0
	 FRNET(I)=0.0
	 BUOYBAR(I)=0.0
	 RHBAR(I,1)=0.0
	 RHBAR(I,2)=0.0 
	 TPBAR(I,1)=0.0
	 TPBAR(I,2)=0.0  
	 PB1(I)=-P(I)
	 DO 52 J=1,2
	  HBAR(I,J)=0.0
	  MDET(I,J)=0.0
	  WATERBAR(I,J)=0.0
	  PB2(I,J)=-P(I)
   52    CONTINUE
	 DO 54 J=1,3
	  MBAR(I,J)=0.0
	  PB3(I,J)=-P(I)
   54    CONTINUE
	 DO 56 J=1,5
	  FTBAR(I,J)=0.0
	  FRBAR(I,J)=0.0
	  PB5(I,J)=-P(I)
   56    CONTINUE
   60   CONTINUE
C
C  Initialize counters
C
	TIME=0.0
	PRINTIME=0.0
	NT=0
	NPR=0
	NAV=0   
	NRH=0
	NTP=0
	RHERR=0.0
	TPERR=0.0
C
C  Begin time loop at next statement
C
  100   CONTINUE
C
	TIME=TIME+DT
	PRINTIME=PRINTIME+DT
	TS1=27.30+273.15
	TS2=TS1  
	NT=NT+1
	IF(TIME.GT.ENDTIME)GOTO 200
C
C  Calculate fraction of observational time step
C
        KT=1
  136   CONTINUE
        KT=KT+1
        IF(TIMER(KT).LT.TIME)GOTO 136
        FRAC=(TIME-TIMER(KT-1))/(TIMER(KT)-TIMER(KT-1))
        FRAC=MAX(FRAC,0.0)

C
C  Calculate saturation specific humidity and some temporary arrays
C   Do crude adjustment for supersaturation
C
        TC=TS2-273.15
        ES=6.112*EXP(17.67*TC/(243.5+TC))
        RSS=0.98*EPS*ES/(PH(1)-(1.-EPS)*ES)
        PRADJ=0.0
C
        DO 110 I=1,NP
         FTNET(I)=0.0
         FRNET(I)=0.0
         TC=T2(I)-273.15
         DENOM=243.5+TC
         DENOM=MAX(DENOM,1.0)
C
C  Calculate saturation specific humidity, assuming water saturation
C    for T> 0C and ice saturation otherwise
C
         IF(TC.GE.0.0)THEN
          ES=6.112*EXP(17.67*TC/DENOM)
         ELSE
          ES=EXP(23.33086-6111.72784/T2(I)+0.15215*LOG(T2(I)))
         END IF
         ES=MIN(ES,P(I))
         RS(I)=EPS*ES/(P(I)-(1.-EPS)*ES)
         TTEMP(I)=T2(I)
         RTEMP(I)=R2(I)
C
         IF(R2(I).GT.RS(I))THEN
          ALV=LV0-CPVMCL*TC
          CPN=CPD*(1.-R2(I))+CPV*R2(I)
          TNEW=(ALV*(R2(I)+RS(I)*(ALV/(RV*T2(I))-1.))+CPN*T2(I))/
     1     (CPN+ALV*ALV*RS(I)/(RV*T2(I)*T2(I)))
          RNEW=RS(I)*(1.+(TNEW-T2(I))*ALV/(RV*T2(I)*T2(I)))
          FRADJ=0.5*(RNEW-R2(I))/DT
          FRNET(I)=FRNET(I)+FRADJ
          FTADJ=-ALV*FRADJ/CPN
          FTNET(I)=FTNET(I)+FTADJ
          PRADJ=PRADJ-100.0*(PH(I)-PH(I+1))*FRADJ*1000.0*
     1     3600.0*24.0/(ROWL*G)
          IF(TIME.GT.(ENDTIME-AVTIME))THEN
           FTBAR(I,1)=FTBAR(I,1)+FTADJ
           FRBAR(I,2)=FRBAR(I,2)+FRADJ
          END IF
         END IF
  110   CONTINUE
      TG=TS2
      RG=RSS
      PG=PSURF
C      
C  Call convection subroutine
C
	CALL CONVECT(T2,R2,RS,P,PH,TG,RG,PG,NA,N,DT,IFLAG,FT,FR,
     1    PRECIP,WD,TPRIME,QPRIME,CBMF)
C
C  Add convective precip to precip from supersaturation adjustment
C
	PRECIP=PRECIP+PRADJ
C
C  Add convective tendencies to net tendencies
C
	DO 125 I=1,N
	 FTNET(I)=FTNET(I)+FT(I)
	 FRNET(I)=FRNET(I)+FR(I)
  125   CONTINUE
C
C  Add radiative cooling to net heating rate
C
        DO 135 I=1,NP
         FTNET(I)=FTNET(I)-RCOOL(KT-1,I)*(1.-FRAC)-RCOOL(KT,I)*FRAC
  135   CONTINUE
C
C  Calculate surface flux forcing
C
	UT=U(KT-1,1)*(1.-FRAC)+U(KT,1)*FRAC
	VT=V(KT-1,1)*(1.-FRAC)+V(KT,1)*FRAC
	WT1=10.*ABS(MP(ICB))*0.01*RD*T2(ICB)/(SIGD*P(ICB))
C       VSURF1=SQRT(VS0*VS0+UT*UT+VT*VT)
      VSU1=SQRT(VS0*VS0+UT*UT+VT*VT)
      VSU2=SQRT(VS0*VS0+UT*UT+VT*VT+WT1*WT1)
      VSURF=VSU1*(1.-SIGMD)+SIGMD*VSU2
      VSPRIME=SIGMD*(VSU2-VSU1)
c      VSURF=SQRT(VS0*VS0+UT*UT+VT*VT+SIGMD*WT1*WT1)
c      VSPRIME=SIGMD*WT1 
c      VSURF=(1.-SIGMD1)*VSURF1+SIGMD1*VSURF2
c       VSURF=VSURF1*(1.+MIN(WT1*WT1/(2.*VSURF1*VSURF1),2.0))
	CD=1.2E-3
	SFRAC=0.5
	TSA=T1(1)*(PSURF/P(1))**(RD/CPD)
	ROWS=PSURF/(RD*TSA*(1.+R1(1)*(1./EPS-1.)))
	FTSURF=G*ROWS*CD*(VSURF*(TS1-TSA)-SFRAC*VSPRIME*
     1 (LV0/CPD)*(RP(1)-R2(1)))/(PH(1)-PH(2))
	FTNET(1)=FTNET(1)+FTSURF
	FRSURF=G*ROWS*CD*(VSURF*(RSS-R1(1))-
     1  (1.-SFRAC)*VSPRIME*(RP(1)-R2(1)))/(PH(1)-PH(2))
C      FRSURF=G*ROWS*CD*((1.-SIGMD1)*VSURF1*(RSS-R1(1))+SIGMD1*
C     1  VSURF2*(RSS-RP(1)))/(PH(1)-PH(2))
	FRNET(1)=FRNET(1)+FRSURF
C
	PRECBAR=PRECBAR+PRECIP
C
C  Put values into time arrays
C
	IF(PRINTIME.GT.(PRINFREQ-10.0))THEN
	 PRINTIME=0.0
	 NPR=NPR+1
	 GTIME(NPR)=TIME/(3600.*24.)
C
	 RAIN_EVAP(NPR,1)=PRECBAR*DT/PRINFREQ
	 RAIN_EVAP(NPR,2)=100.*ROWS*CD*(VSURF*(RSS-R1(1))-(1.-SFRAC)* 
     1   VSPRIME*(RP(1)-R2(1)))*1000.0*3600.0*24./ROWL
C       RAIN_EVAP(NPR,2)=100.*ROWS*CD*((1.-SIGMD1)*VSURF1*(RSS-R1(1))+
C     1   SIGMD1*VSURF2*(RSS-RP(1)))*1000.*3600.*24./ROWL
C
	 SST_TS(NPR,1)=TSA-273.15
	 SST_TS(NPR,2)=TS1-273.15
	 Q500(NPR,1)=1000.*(RR(KT-1,9)*(1.-FRAC)+FRAC*RR(KT,9))
	 Q500(NPR,2)=1000.*R2(9)
	 Q850(NPR,1)=1000.*(RR(KT-1,4)*(1.-FRAC)+FRAC*RR(KT,4))
	 Q850(NPR,2)=1000.*R2(4)
	 T200(NPR,1)=TR(KT-1,15)*(1.-FRAC)+FRAC*TR(KT,15)-273.15
	 T200(NPR,2)=T2(15)-273.15
	 T850(NPR,1)=TR(KT-1,4)*(1.-FRAC)+FRAC*TR(KT,4)-273.15
	 T850(NPR,2)=T2(4)-273.15
	 PBUDGET(NPR)=-1.0E5*FRSUM*DT*3600.*24./(ROWL*G*PRINFREQ)
	 PBUDGET(NPR)=MAX(PBUDGET(NPR),0.0)
	 DOWNV(NPR)=WT1
	 FESUM=FHSUM*DT/PRINFREQ 
	 ALV=LV0-CPVMCL*SST_TS(NPR,1)
	 EBUDGET(NPR)=1.0E5*FESUM*3600.*24./(ROWL*G*ALV) 
	 FRSUM=0.0
	 PRECBAR=0.0
	 FHSUM=0.0   
	 ENTS=0.0
         MSE(NPR,1)=0.0
         MSE(NPR,2)=0.0	
	 DO 138 I=1,N
          TTEM=TR(KT-1,I)*(1.-FRAC)+TR(KT,I)*FRAC
          RTEM=RR(KT-1,I)*(1.-FRAC)+RR(KT,I)*FRAC
	  RHTIME(NPR,I)=100.*R2(I)/RS(I)
	  RHTIME(NPR,I)=MIN(RHTIME(NPR,I),100.0)
	  RHR=(RR(KT-1,I)/RSR(KT-1,I))*(1.-FRAC)+
     1     (RR(KT,I)/RSR(KT,I))*FRAC
	  RHR=MIN(RHR,1.0)
	  RHTIMR(NPR,I)=100.*RHR
	  RHDIF(NPR,I)=RHTIME(NPR,I)-RHTIMR(NPR,I)
	  OMEGA(NPR,I)=3600.*(WT(KT-1,I)*(1.-FRAC)+
     1     WT(KT,I)*FRAC)
          ALV=LV0-CPVMCL*(TTEM-273.15)
          ALV1=LV0-CPVMCL*(T2(I)-273.15)
          MSE(NPR,1)=MSE(NPR,1)+((CPD*(1.-RTEM)+RTEM*CL)*TTEM+
     1     ALV*RTEM)*(PH(I)-PH(I+1))/G
        MSE(NPR,2)=MSE(NPR,2)+((CPD*(1.-R2(I))+R2(I)*CL)*T2(I)
     1   +ALV1*R2(I))*(PH(I)-PH(I+1))/G
	TPQTIME(NPR,I)=(TR(KT-1,I)-TR(1,I))*(1.-FRAC)+
     1   (TR(KT,I)-TR(1,I))*FRAC
	TPMTIME(NPR,I)=T2(I)-TR(1,I)
	MUTIME(NPR,I)=MUP(I)*1000.
	MDTIME(NPR,I)=-MP(I)*1000.
	MPENTIME(NPR,I)=-MDOWN(I)*1000.   
          ALV=LV0-CPVMCL*(T2(I)-273.15)
          ENTS=ENTS+((CPD*(1.-R2(I))+CPV*R2(I))*FT(I)+(ALV+
     1     (CL-CPD)*(T2(I)-T2(1)))*FR(I))*(PH(I)-PH(I+1))
  138    CONTINUE
         ENTHALPY(NPR)=ENTS*3600.*24./(CPD*(PH(1)-PH(N+1)))
         MSE(NPR,1)=MSE(NPR,1)*G/(900.*CPD)
         MSE(NPR,2)=MSE(NPR,2)*G/(900.*CPD)
	END IF
C
C  Calculate vertical advection and add to net tendencies
C
	WOMEG=WT(KT-1,2)*(1.-FRAC)+WT(KT,2)*FRAC
	VVP=MAX(WOMEG,0.0)
	ALPHA=RD*T2(1)*(1.+R2(1)/EPS)/(P(1)*(1.+R2(1)))
	FTADV=VVP*ALPHA/CPD-VVP*(T2(1)-T2(2))/(P(1)-P(2))
	FTNET(1)=FTNET(1)+FTADV
	FRADV=-VVP*(R2(1)-R2(2))/(P(1)-P(2))
	RRB1=RR(KT-1,1)*(1.-FRAC)+RR(KT,1)*FRAC
	RRB2=RR(KT-1,2)*(1.-FRAC)+RR(KT,2)*FRAC
	TRB1=TR(KT-1,1)*(1.-FRAC)+TR(KT,1)*FRAC
	TRB2=TR(KT-1,2)*(1.-FRAC)+TR(KT,2)*FRAC
	FRRADV=-VVP*(RRB1-RRB2)/(P(1)-P(2))    
	FTRADV=-VVP*(TRB1-TRB2)/(P(1)-P(2))
	ALV=LV0-CPVMCL*(TRB1-273.15)
	FRNET(1)=FRNET(1)+FRADV
	FRADVNET=(-FRRADV+ADVR(KT-1,1)*(1.-FRAC)+ADVR(KT,1)*FRAC+
     1  (RR(KT,1)-RR(KT-1,1))/(TIMER(KT)-TIMER(KT-1)))*(PH(1)-PH(2))
	FRSUM=FRSUM+FRADVNET-G*ROWS*CD*(VSURF*(RSS-RRB1)-
     1 (1.-SFRAC)*VSPRIME*(RP(1)-R2(1)))    
C      FRSUM=FRSUM+FRADVNET-G*ROWS*CD*((1.-SIGMD1)*VSURF1*(RSS-RRB1)+
C     1  SIGMD1*VSURF2*(RSS-RP(1))) 
	RADC=RCOOL(KT-1,1)*(1.-FRAC)+RCOOL(KT,1)*FRAC
      FTADVNET=(-FTRADV+RADC+ADVT(KT-1,1)*(1.-FRAC)+ADVT(KT,1)*FRAC+
     1 (TR(KT,1)-TR(KT-1,1))/(TIMER(KT)-TIMER(KT-1)))*(PH(1)-PH(2)) 
      ESUM=FTADVNET*(CPD*(1.-RRB1)+CPV*RRB1)+ALV*FRADVNET-
     1  VVP*RD*0.5*(TRB1+TRB2)*(PH(1)-PH(2))/PH(2)
      FHSUM=FHSUM+ESUM 
	DO 140 I=2, NP-1
	 WOMEG=WT(KT-1,I)*(1.-FRAC)+WT(KT,I)*FRAC
	 WP=WT(KT-1,I+1)*(1.-FRAC)+WT(KT,I+1)*FRAC
	 VVP=MAX(WP,0.0)
	 VVM=MIN(WOMEG,0.0)
	 ALPHA=RD*T2(I)*(1.+R2(I)/EPS)/(P(I)*(1.+R2(I)))
	 FTADV=(VVP+VVM)*ALPHA/CPD-
     1    VVP*(T2(I)-T2(I+1))/(P(I)-P(I+1))-
     2    VVM*(T2(I-1)-T2(I))/(P(I-1)-P(I))
C
	 FTNET(I)=FTNET(I)+FTADV
C
	 FRADV=-VVP*(R2(I)-R2(I+1))/(P(I)-P(I+1))-    
     1     VVM*(R2(I-1)-R2(I))/(P(I-1)-P(I))
C                             
       RRB=RR(KT-1,I)*(1.-FRAC)+RR(KT,I)*FRAC
       RRBP=RR(KT-1,I+1)*(1.-FRAC)+RR(KT,I+1)*FRAC
       RRBM=RR(KT-1,I-1)*(1.-FRAC)+RR(KT,I-1)*FRAC
       TRB=TR(KT-1,I)*(1.-FRAC)+TR(KT,I)*FRAC
       TRBP=TR(KT-1,I+1)*(1.-FRAC)+TR(KT,I+1)*FRAC
       TRBM=TR(KT-1,I-1)*(1.-FRAC)+TR(KT,I-1)*FRAC
       FRRADV=-VVP*(RRB-RRBP)/(P(I)-P(I+1))-VVM*(RRBM-RRB)/(P(I-1)-P(I))
       FTRADV=-VVP*(TRB-TRBP)/(P(I)-P(I+1))-VVM*(TRBM-TRB)/(P(I-1)-P(I))
       FRADVNET=(-FRRADV+ADVR(KT-1,I)*(1.-FRAC)+ADVR(KT,I)*FRAC+
     1 (RR(KT,I)-RR(KT-1,I))/(TIMER(KT)-TIMER(KT-1)))*(PH(I)-PH(I+1))
	RADC=RCOOL(KT-1,I)*(1.-FRAC)+RCOOL(KT,I)*FRAC
       FTADVNET=(-FTRADV+RADC+ADVT(KT-1,I)*(1.-FRAC)+ADVT(KT,I)*
     1  FRAC+
     2 (TR(KT,I)-TR(KT-1,I))/(TIMER(KT)-TIMER(KT-1)))*(PH(I)-PH(I+1))
       FRSUM=FRSUM+FRADVNET
       ALV=LV0-CPVMCL*(TRB-273.15)
       TEMP=ALV*FRADVNET+(CPD*(1.-RRB)+CPV*RRB)*FTADVNET-
     1  (RD*VVP*0.5*(TRB+TRBP)/PH(I+1)+RD*VVM*0.5*(TRB+TRBM)/PH(I))*
     2  (PH(I)-PH(I+1))
       FHSUM=FHSUM+TEMP
       ESUM=ESUM+TEMP  
   73	 FORMAT(F8.2,1P,7E11.3)
   70	CONTINUE
C
       FRNET(I)=FRNET(I)+FRADV
C
C  Calculate certain average quantities for output graphs
C
	 IF(TIME.GT.(ENDTIME-AVTIME-10.0))THEN
	  FTBAR(I,4)=FTBAR(I,4)+FTADV
	  FRBAR(I,4)=FRBAR(I,4)+FRADV
	 END IF
  140   CONTINUE
	 RTR=RR(KT-1,1)*(1.-FRAC)+RR(KT,1)*FRAC                          
         FHSUM=FHSUM-ROWS*G*CD*(VSURF*(TS1-TSA)-SFRAC*VSPRIME*
     1   (LV0/CPD)*(RP(1)-R2(1)))*
     2   (CPD*(1.-RTR)+CPV*RTR) 
C
	TRB1=TR(KT-1,1)*(1.-FRAC)+TR(KT,1)*FRAC
	RRB1=RR(KT-1,1)*(1.-FRAC)+RR(KT,1)*FRAC 
	ALV=LV0-CPVMCL*(TRB1-273.15)
c       FRSURF=(ESUM-CD*G*VSURF*ROWS*(TS1-TRB1)*(CPD*(1.-RRB1)+
c     1 CPV*RRB1))/(ALV*(PH(1)-PH(2)))
c       FRNET(1)=FRNET(1)+FRSURF   
C
C  Add horizontal advections
C
	DO 145 I=1,NP
	 FADVT=ADVT(KT-1,I)*(1.-FRAC)+ADVT(KT,I)*FRAC
	 FADVR=ADVR(KT-1,I)*(1.-FRAC)+ADVR(KT,I)*FRAC
	 FTNET(I)=FTNET(I)-FADVT
	 FRNET(I)=FRNET(I)-FADVR
	 IF(TIME.GT.(ENDTIME-AVTIME-10.0))THEN
	  FTBAR(I,5)=FTBAR(I,5)-FADVT
	  FRBAR(I,5)=FRBAR(I,5)-FADVR
	 END IF
  145   CONTINUE
C
C  Calculate root mean square relative humidity error
C
	DO 148 I=1,NP
	 IF(P(I).GE.290.)THEN
	  RH=R2(I)/RS(I)
	  RHR=(RR(KT-1,I)/RSR(KT-1,I))*(1.-FRAC)+
     1     (RR(KT,I)/RSR(KT,I))*FRAC
	  RHERR=RHERR+(RH-RHR)**2
	  NRH=NRH+1
	 END IF
	 TPR=(TR(KT-1,I))*(1.-FRAC)+TR(KT,I)*FRAC
	 TPERR=TPERR+(T2(I)-TPR)**2
	 NTP=NTP+1   
  148   CONTINUE
C
C  Calculate certain quantities averaged at end of integration
C
	IF(TIME.GT.(ENDTIME-AVTIME))THEN
	 NAV=NAV+1
	 FTBAR(1,3)=FTBAR(1,3)+FTSURF
	 FRBAR(1,3)=FRBAR(1,3)+FRSURF
	 GZ=0.0 
	 COEF=0.5*(1.-2.*AFILT)/DT
	 DO 150 I=1,NP
	  FTBAR(I,1)=FTBAR(I,1)+FT(I)
	  FTBAR(I,3)=FTBAR(I,3)+COEF*(T2(I)-TTEMP(I))
	  FTBAR(I,2)=FTBAR(I,2)-RCOOL(KT-1,I)*(1.-FRAC)-FRAC*RCOOL(KT,I)
	  FRBAR(I,1)=FRBAR(I,1)+FRDET(I)
	  FRBAR(I,2)=FRBAR(I,2)+FRSUB(I)
	  FRBAR(I,3)=FRBAR(I,3)+COEF*(R2(I)-RTEMP(I))
	  BUOYBAR(I)=BUOYBAR(I)+TVP(I)-T2(I)*(1.+R2(I)/EPS)/
     1      (1.+R2(I))
	  RHBAR(I,1)=RHBAR(I,1)+R2(I)/MAX(RS(I),1.0E-6)
	  RHBAR(I,2)=RHBAR(I,2)+(RR(KT-1,I)/RSR(KT-1,I))*(1.-FRAC)+
     1     (RR(KT,I)/RSR(KT,I))*FRAC
	TPBAR(I,1)=TPBAR(I,1)+T2(I)-TR(1,I)
	TPBAR(I,2)=TPBAR(I,2)+(TR(KT-1,I))*(1.-FRAC)+FRAC*
     1   TR(KT,I)-TR(1,I)
	  IF(I.GT.1)GZ=GZ+0.5*RD*(T2(I-1)*(1.+R2(I-1)/EPS)/
     1      (1.+R2(I-1))+T2(I)*(1.+R2(I)/EPS)/(1.+R2(I)))*
     2      (P(I-1)-P(I))/PH(I)
	  ALV=LV0-CPVMCL*(T2(I)-273.15)
	  HBAR(I,1)=HBAR(I,1)+(CPD*(1.-R2(I))+CPV*R2(I))*T2(I)+
     1     ALV*R2(I)+GZ
	  HBAR(I,2)=HBAR(I,2)+(CPD*(1.-RS(I))+CPV*RS(I))*T2(I)+
     1     ALV*RS(I)+GZ
	  MBAR(I,1)=MBAR(I,1)+MUP(I)
	  MBAR(I,2)=MBAR(I,2)-MDOWN(I)
	  MBAR(I,3)=MBAR(I,3)-MP(I)
	  MDET(I,1)=MDET(I,1)+ENT(I)
	  MDET(I,2)=MDET(I,2)+DET(I)
	  WATERBAR(I,1)=WATERBAR(I,1)+WATER(I)
	  WATERBAR(I,2)=WATERBAR(I,2)+RP(I)-R2(I)
  150    CONTINUE
	END IF
C
C  Rectify temperature to force enthalpy conservation
C
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      SUM4=0.0
      SUM5=0.0
      TTA1=TR(KT-1,1)*(1.-FRAC)+FRAC*TR(KT,1)
      DO 155 I=1,N
       DELPH=PH(I)-PH(I+1)
       RTA=RR(KT-1,I)*(1.-FRAC)+FRAC*RR(KT,I)
       TTA=TR(KT-1,I)*(1.-FRAC)+FRAC*TR(KT,I)
       ALVA=LV0-CPVMCL*(TTA-273.15)
       ALV=LV0-CPVMCL*(T2(I)-273.15)
       SUM1=SUM1+(CPD*(1.-RTA)+CL*RTA)*TTA*DELPH
       SUM2=SUM2+ALVA*RTA*DELPH
       SUM3=SUM3+(CPD*(1.-R2(I))+CL*R2(I))*T2(I)*DELPH
       SUM4=SUM4+ALV*R2(I)*DELPH
       SUM5=SUM5+(CPD*(1.-R2(I))+CL*R2(I))*DELPH
  155 CONTINUE
C      ACOR=(SUM1+SUM2-SUM3-SUM4)/SUM5
      ACOR=(SUM1+SUM2-SUM3-SUM4)/((LV0-CPVMCL*(T2(1)-273.15))
     1  *(PH(1)-PH(2)))
      FRNET(1)=FRNET(1)+ACOR/(TAUF*24.*3600.)
      ALH1BAR=ALH1BAR+ACOR/(TAUF*24.*3600.)
      ACOR2=(SUM1+SUM2-SUM3-SUM4)/SUM5
      DO 158 I=1,N
       FTNET(I)=FTNET(I)+ACOR2/(TAUR*24.*3600.)
         IF(TIME.GT.(ENDTIME-AVTIME))THEN
        FTBAR(I,2)=FTBAR(I,2)+ACOR2/(TAUR*24.*3600.)
       END IF
  158 CONTINUE
C
C  Advance quantities one time step, applying Asselin filter
C
	DO 160 I=1,NP
	 IF(P(I).GE.PFIX)THEN
	  T3(I)=T1(I)+2.*DT*FTNET(I)
	  R3(I)=R1(I)+2.*DT*FRNET(I)
	  R3(I)=MAX(R3(I),2.0E-6)
	  T1(I)=T2(I)+AFILT*(T1(I)+T3(I)-2.*T2(I))
	  R1(I)=R2(I)+AFILT*(R1(I)+R3(I)-2.*R2(I))
	  T2(I)=T3(I)
	  R2(I)=R3(I)
	 END IF
  160   CONTINUE
C
C  Return to beginning of time loop
C
	GOTO 100
C
C  Output processing starts here
C
  200   CONTINUE
C
C  Write rms error of relative humidity to output file
C
	WRITE(12,209)ELCRIT,TLCRIT,SIGD,SIGS,VS0
  209   FORMAT(2X,'ELCRIT=',F11.5,' TLCRIT=',F9.3,' SIGD=',F6.4,
     1   ' SIGS=',F6.4,' VS0=',F5.1)
      WRITE(12,211)SIGMD,WTRAIN,OMTSNOW,COEFFR,COEFFS
  211 FORMAT(2X,'sig_md=',F6.4,'  wt_rain=',F6.3,' wt_snow=',
     1  F6.3,'  coeffr=',F6.2,'  coeffs=',F6.2)
      WRITE(12,214)DTMAX,SLFRAC,ENTP
  214 FORMAT(3X,'dTmax=',F6.4,'  slfrac=',F6.4,'  entp=',F6.4)        
	RHERR=100.*SQRT(RHERR/FLOAT(NRH))
	TPERR=SQRT(TPERR/FLOAT(NTP))
	WRITE(12,210)RHERR,TPERR
  210   FORMAT(5X,'Root mean square RH error =  ', F6.2,
     1    2X,'rms T error = ',F6.2)
C
C  NCAR graphics are used to make the plots. One may substitute one's own
C     graphics routines at this point.
C
C  Normalize averaged quantities
C
	FAC=1./FLOAT(NAV)
	DO 300 I=1,NP
	 BUOYBAR(I)=FAC*BUOYBAR(I)
	 BUOYBAR(I)=MAX(BUOYBAR(I),-4.0)
	 RHBAR(I,1)=FAC*RHBAR(I,1)
	 RHBAR(I,1)=MIN(RHBAR(I,1),1.0)
	 RHBAR(I,2)=FAC*RHBAR(I,2)
	 TPBAR(I,1)=FAC*TPBAR(I,1)
	 TPBAR(I,2)=FAC*TPBAR(I,2)  
	 DO 270 J=1,2
	  HBAR(I,J)=FAC*HBAR(I,J)/CPD
	  HBAR(I,J)=MIN(HBAR(I,J),HBAR(1,2))
	  MDET(I,J)=FAC*MDET(I,J)*1000.0
	  WATERBAR(I,J)=FAC*WATERBAR(I,J)*1000.0
  270 CONTINUE
	 DO 280 J=1,3
	  MBAR(I,J)=FAC*MBAR(I,J)*1000.0
  280 CONTINUE
	 DO 290 J=1,5
	  FTBAR(I,J)=FAC*FTBAR(I,J)*24.*3600.
	  FRBAR(I,J)=FAC*FRBAR(I,J)*1000.*24.*3600.
  290 CONTINUE
C
C  Clean up SIJ array for display purposes
C
	 DO 295 J=1,NP
	  IF(SIJ(I,J).GE.1.0)SIJ(I,J)=0.0
	  SIJ(I,J)=MAX(SIJ(I,J),0.0)
	  SIJ(I,J)=SIJ(I,J)-0.001
  295    CONTINUE
  300   CONTINUE
C
C  Prepare dash line patterns for NCAR graphics
C
        DSH5(1)='$$$$'
        DSH5(2)='''$$'''
        DSH5(3)='''*''*'
        DSH5(4)='''+'''''
        DSH5(5)='''H'''''
        DSH3(1)='$$$$'
        DSH3(2)=DSH5(2)
        DSH3(3)=DSH5(3)
        DSH2(1)=DSH5(1)
        DSH2(2)=DSH5(2)
        XLAB='Time (days)'
        YLAB='P (mb)'
C
C  Open and use NCAR graphics
C
        CALL OPNGKS
        CALL ANOTAT('TIME (DAYS)$', 
     1   'PRECIPITATION AND EVAPORATION (MM/DAY)$',1,1,2,DSH2)
        CALL DISPLA(1,1,0)
        CALL EZMXY(GTIME,RAIN_EVAP,NTIME,2,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$','SURFACE TEMP. AND SST(C)$',
     1   1,1,2,DSH2)
        CALL EZMXY(GTIME,SST_TS,NTIME,2,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$',
     1   'Integrated Enthalpy(K)$',1,1,2,DSH2)
        CALL DISPLA(1,1,0)
        CALL EZMXY(GTIME,MSE,NTIME,2,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$','Enthalpy Tendency$',
     1    1,1,2,DSH2)
        CALL DISPLA(1,1,0)
        CALL EZXY(GTIME,ENTHALPY,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$','Downdraft Velocity (m/s)$',
     1    1,1,2,DSH2)
        CALL DISPLA(1,1,0)
        CALL EZXY(GTIME,DOWNV,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$',
     1   '500 mb specific humidity (gm/kg)$',1,1,2,DSH2)
        CALL EZMXY(GTIME,Q500,NTIME,2,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$',
     1   '850 mb specific humidity (gm/kg)$',1,1,2,DSH2)
        CALL EZMXY(GTIME,Q850,NTIME,2,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$','200 mb temperature (C)$',
     1   1,1,2,DSH2)
        CALL EZMXY(GTIME,T200,NTIME,2,NPR,'$')
        CALL ANOTAT('TIME (DAYS)$','850 mb temperature (C)$',
     1   1,1,2,DSH2)
        CALL EZMXY(GTIME,T850,NTIME,2,NPR,'$')
        CALL ANOTAT('TEMP. TENDENCY (DEG C/DAY)$','PRESSURE (MB)$',
     1    1,1,5,DSH5)
        CALL DISPLA(1,2,0)
        CALL EZMXY(FTBAR,PB5,NA,5,NP,'$')
        CALL ANOTAT('SPECIFIC HUMIDITY TENDENCY (G/KG/DAY)$',
     1    'PRESSURE (MB)$',1,1,5,DSH5)
        CALL DISPLA(1,2,0)
        CALL EZMXY(FRBAR,PB5,NA,5,NP,'$')
        CALL ANOTAT('BUOYANCY (K)$','PRESSURE (MB)$',1,1,0,1)
        CALL DISPLA(1,1,0)
        CALL EZXY(BUOYBAR,PB1,N,'$')
        CALL ANOTAT('RELATIVE HUMIDITY$','PRESSURE (MB)$',1,1,0,1)
        CALL DISPLA(1,2,0)
        CALL EZMXY(RHBAR,PB2,NA,2,NP,'$')
        CALL ANOTAT('Temperature$','PRESSURE (MB)$',1,1,0,1)
        CALL DISPLA(1,2,0)
        CALL EZMXY(TPBAR,PB2,NA,2,NP,'$')
        CALL ANOTAT('h/c_p, h*/c_p (K)$','PRESSURE (MB)$',1,1,2,DSH2)
        CALL DISPLA(1,2,0)
        CALL EZMXY(HBAR,PB2,NA,2,N,'$')
        CALL ANOTAT('NET MASS FLUXES$','PRESSURE (MB)$',1,1,3,DSH3)
        CALL DISPLA(1,2,0)
        CALL EZMXY(MBAR,PB3,NA,3,N,'$')
        CALL ANOTAT('MASS ENTRAINMENT AND DETRAINMENT$',
     1   'PRESSURE (MB)$',1,1,2,DSH2)
        CALL DISPLA(1,2,0)
        CALL EZMXY(MDET,PB2,NA,2,N,'$')
c        CALL ANOTAT('RAIN WATER, RP(I)-R(I) (GM/KG)$',
c     1  'PRESSURE (MB)$',1,1,2,DSH2)
c        CALL DISPLA(1,2,0)
c        CALL EZMXY(WATERBAR,PB2,NA,2,N,'$')
c        CALL ANOTAT('SIGMA_U (*1000)$','PRESSURE (MB)$',1,1,0,1)
c        CALL DISPLA(1,1,0)
c        CALL EZXY(SIGBAR,PB1,N,'$')
        ZLAB='RH_calc'
        CALL PLREG(GTIME,PB1,RHTIME,NTIME,NA,NPR,N,
     1    1.0,1.0,XLAB,YLAB,ZLAB)
        ZLAB='RH_real'
        CALL PLREG(GTIME,PB1,RHTIMR,NTIME,NA,NPR,N,
     1    1.0,1.0,XLAB,YLAB,ZLAB)
        ZLAB='Delta_RH'
        CALL PLREG(GTIME,PB1,RHDIF,NTIME,NA,NPR,29,
     1    1.0,1.0,XLAB,YLAB,ZLAB)
        ZLAB='Omega (mb/hr)'
        CALL PLREG(GTIME,PB1,OMEGA,NTIME,NA,NPR,N,
     1    1.0,1.0,XLAB,YLAB,ZLAB)
        CALL CLSGKS
C
C  Write to output files
C
      DO 310 I=1,NPR
       WRITE(18,305)GTIME(I),RAIN_EVAP(I,1),SST_TS(I,1),SST_TS(I,2),
     2  Q500(I,1),Q500(I,2)
  305  FORMAT(1X,5(F9.3,','),F9.3)
       DO 308 J=1,NP
        WRITE(20,307)GTIME(I),P(J),RHTIMR(I,J),RHTIME(I,J),RHDIF(I,J),
     1   TPMTIME(I,J),MPENTIME(I,J)
  307  FORMAT(1X,F8.2,6(',',F8.2))
  308  CONTINUE
  310 CONTINUE
C
      DO 320 I=1,NP
       WRITE(19,315)P(I),BUOYBAR(I),RHBAR(I,2),RHBAR(I,1),TPBAR(I,2),
     1  TPBAR(I,1),FTBAR(I,2),MBAR(I,1),MBAR(I,2),MBAR(I,3),
     2  MDET(I,1),MDET(I,2)
  315  FORMAT(1X,11(F9.3,','),F9.3)
  320 CONTINUE

	STOP 
	END
C
***************************************************************************
*****                       SUBROUTINE CONVECT                        *****
*****                           VERSION 4.2                           *****
*****                          15 July, 1998                          *****
***************************************************************************
C
	SUBROUTINE CONVECT
     *    (T,       Q,   QS,     P,       PH,      TG,   QG,
     *     PG,      ND,  NL,     DELT,    IFLAG,   FT,   FQ,
     *     PRECIP,  WD,  TPRIME, QPRIME,  CBMF)
C
C-----------------------------------------------------------------------------
C    *** On input:      ***
C
C     T:   Array of absolute temperature (K) of dimension ND, with first
C           index corresponding to lowest model level. Note that this array
C           will be altered by the subroutine if dry convective adjustment
C           occurs
C
C     Q:   Array of specific humidity (gm/gm) of dimension ND, with first
C            index corresponding to lowest model level. Must be defined
C            at same grid levels as T. Note that this array will be altered
C            if dry convective adjustment occurs
C
C     QS:  Array of saturation specific humidity of dimension ND, with first
C            index corresponding to lowest model level. Must be defined
C            at same grid levels as T. Note that this array will be altered
C            if dry convective adjustment occurs
C
C     P:   Array of pressure (mb) of dimension ND, with first
C            index corresponding to lowest model level. Must be defined
C            at same grid levels as T
C
C     PH:  Array of pressure (mb) of dimension ND+1, with first index
C            corresponding to lowest level. These pressures are defined at
C            levels intermediate between those of P, T, Q and QS. The first
C             value of PH should be greater than (i.e. at a lower level than)
C             the first value of the array P
C
C     TG:  Scalar ground or sea surface temperature (K)
C
C     QG:  Scalar ground or sea surface value of saturation specific
C            humidity(gm/gm)
C
C     PG:  Scalar ground or sea surface value of pressure (mb). Can be the same
C           value as PH(1)
C
C     ND: The dimension of the arrays T,Q,QS,P,PH,FT and FQ
C
C     NL:  The maximum number of levels to which convection can
C              penetrate, plus 1
C              NL MUST be less than or equal to ND-1
C
C     DELT: The model time step (sec) between calls to CONVECT
C
C----------------------------------------------------------------------------
C    ***   On Output:         ***
C
C     IFLAG: An output integer whose value denotes the following:
C
C                VALUE                        INTERPRETATION
C                -----                        --------------
C                  0               No moist convection; atmosphere is not
C                                  unstable, or surface temperature is less
C                                  than 250 K or surface specific humidity
C                                  is non-positive.
C
C                  1               Moist convection occurs.
C
C                  2               No moist convection: lifted condensation
C                                  level is above the 200 mb level.
C
C                  3               No moist convection: cloud base is higher
C                                  then the level NL-1.
C
C                  4               Moist convection occurs, but a CFL condition
C                                  on the subsidence warming is violated. This
C                                  does not cause the scheme to terminate.
C
C     FT:   Array of temperature tendency (K/s) of dimension ND, defined at same
C             grid levels as T, Q, QS and P
C
C     FQ:   Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
C             defined at same grid levels as T, Q, QS and P
C
C     PRECIP: Scalar convective precipitation rate (mm/day)
C
C     WD:    A convective downdraft velocity scale. For use in surface
C       flux parameterizations. See convect.ps file for details
C
C     TPRIME: A convective downdraft temperature perturbation scale (K).
C             For use in surface flux parameterizations. See convect.ps
C             file for details
C
C     QPRIME: A convective downdraft specific humidity
C              perturbation scale (gm/gm).
C              For use in surface flux parameterizations. See convect.ps
C              file for details
C
C     CBMF:   The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
C              BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
C              ITS NEXT CALL. That is, the value of CBMF must be "remembered"
C              by the calling program between calls to CONVECT
C
C------------------------------------------------------------------------------
C
C    ***  THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN   ***
C    ***                OR EQUAL TO  ND + 1                    ***
C
      PARAMETER (NA=50)           
C
      INTEGER NENT(NA)
      REAL T(ND),Q(ND),QS(ND),P(ND),PH(ND+1)
      REAL FT(ND),FQ(ND)
      REAL M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA)
      REAL SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA)
      REAL QP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA)
      REAL SIGP(NA),TP(NA),TOLD(NA),CPN(NA)
      REAL LV(NA),LVCP(NA),LV0,H(NA),HP(NA),GZ(NA),HM(NA)
C
      REAL MUP(NA),MDOWN(NA),ENT(NA),DET(NA),FQDET(NA),FQSUB(NA)
      COMMON / CVT / TVP,MP,QP,SIJ,MUP,MDOWN,ENT,DET,FQDET,FQSUB,
     1     WATER
C
C------------------------------------------------------------------------------
C
C   ***                    SPECIFY PARAMETERS                        ***
C
C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
C   ***               BETWEEN 0 C AND TLCRIT)                        ***
C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
C   ***                       FORMULATION                            ***
C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
C   ***                        OF CLOUD                              ***
C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
C   ***                          OF RAIN                             ***
C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
c   ***                          OF SNOW                             ***
C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
C
        ELCRIT=.0011
        TLCRIT=-55.0
        ENTP=1.5
        SIGD=0.05
        SIGS=0.12
        OMTRAIN=50.0
        OMTSNOW=5.5 
        COEFFR=1.0
        COEFFS=0.8
	BETA=10.0
        DTMAX=0.85
        ALPHA=0.2  
        DAMP=0.1
C
C   ***        ASSIGN VALUES OF THERMODYNAMIC CONSTANTS,        ***
C   ***            GRAVITY, AND LIQUID WATER DENSITY.           ***
C   ***             THESE SHOULD BE CONSISTENT WITH             ***
C   ***              THOSE USED IN CALLING PROGRAM              ***
C   ***     NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT  ***
C
      CPD=1005.7
      CPV=1870.0
      CL=4190.0 
      RV=461.5
      RD=287.04
      LV0=2.501E6
      G=9.8  
      ROWL=1000.0
C
      CPVMCL=CL-CPV 
      EPS=RD/RV
      EPSI=1./EPS
      GINV=1.0/G
      DELTI=1.0/DELT
C
C           ***  INITIALIZE OUTPUT ARRAYS AND PARAMETERS  ***
C
	DO 5 I=1,ND
	 FT(I)=0.0
	 FQ(I)=0.0
	 MDOWN(I)=0.0
	 MUP(I)=0.0
	 ENT(I)=0.0
	 DET(I)=0.0
	 FQDET(I)=0.0
	 FQSUB(I)=0.0
	 WATER(I)=0.0
	 QP(I)=Q(I)
	 MP(I)=0.0
	 TVP(I)=T(I)*(1.+Q(I)*EPSI-Q(I))
    5   CONTINUE
	DO 7 I=1,NL+1
	 RDCP=(RD*(1.-Q(I))+Q(I)*RV)/
     1    (CPD*(1.-Q(I))+Q(I)*CPV)
	 TH(I)=T(I)*(1000.0/P(I))**RDCP
    7   CONTINUE
	PRECIP=0.0
        WD=0.0
        TPRIME=0.0
        QPRIME=0.0
	IFLAG=0
	INB=ND
C
C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
C
C     ***  DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A  ***
C     ***                BOUNDARY LAYER SCHEME !                 ***
C                        
      JC=0
	DO 30 I=NL-1,1,-1
	 JN=0
	   SUM=TH(I)
	 DO 10 J=I+1,NL
	    SUM=SUM+TH(J)
	    THBAR=SUM/FLOAT(J+1-I)
	  IF(TH(J).LT.THBAR)JN=J
   10    CONTINUE
         IF(I.EQ.1)JN=MAX(JN,2)
	 IF(JN.EQ.0)GOTO 30
   12    CONTINUE
	 AHM=0.0
	 RM=0.0
	 DO 15 J=I,JN
	  AHM=AHM+(CPD*(1.-Q(J))+Q(J)*CPV)*T(J)*(PH(J)-PH(J+1))
	  RM=RM+Q(J)*(PH(J)-PH(J+1))
   15    CONTINUE
	 RM=RM/(PH(I)-PH(JN+1))
	 A2=0.0
	 DO 20 J=I,JN
	  Q(J)=RM
	  RDCP=(RD*(1.-Q(J))+Q(J)*RV)/
     1     (CPD*(1.-Q(J))+Q(J)*CPV)  
	  X=(0.001*P(J))**RDCP
	  TOLD(J)=T(J)
	  T(J)=X
	  A2=A2+(CPD*(1.-Q(J))+Q(J)*CPV)*X*(PH(J)-PH(J+1))
   20    CONTINUE
	 DO 25 J=I,JN
	  TH(J)=AHM/A2
	  T(J)=T(J)*TH(J)
	  TC=TOLD(J)-273.15
	  ALV=LV0-CPVMCL*TC
	  QS(J)=QS(J)+QS(J)*(1.+QS(J)*0.608)*ALV*(T(J)-
     1     TOLD(J))/(RV*TOLD(J)*TOLD(J))
   25    CONTINUE
	 IF(TH(JN+1).LT.TH(JN))THEN
	  JN=JN+1
	  GOTO 12
	 END IF
	 IF(I.EQ.1)JC=JN 
   30   CONTINUE
C
C   ***   Remove any supersaturation that results from adjustment ***
C
      IF(JC.GT.1)THEN
       DO 38 J=1,JC
	  IF(QS(J).LT.Q(J))THEN 
	   ALV=LV0-CPVMCL*(T(J)-273.15)  
	   TNEW=T(J)+ALV*(Q(J)-QS(J))/(CPD*(1.-Q(J))+
     1      CL*Q(J)+QS(J)*(CPV-CL+ALV*ALV/(RV*T(J)*T(J))))
	   QNEW=Q(J)-(CPD*(1.-Q(J))+CPV*Q(J))*(TNEW-T(J))/ALV
	   PRECIP=PRECIP+24.*3600.*1.0E5*(PH(J)-PH(J+1))*
     1      (Q(J)-QNEW)/(G*DELT*ROWL)
	   T(J)=TNEW
	   Q(J)=QNEW
	  END IF     
   38  CONTINUE  
      END IF
C
C  *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY
C  
	GZ(1)=0.0
	CPN(1)=CPD*(1.-Q(1))+Q(1)*CPV
	H(1)=T(1)*CPN(1)
        LV(1)=LV0-CPVMCL*(T(1)-273.15)
        HM(1)=H(1)+LV(1)*Q(1)
        AHMIN=1.0E12
        IHMIN=NL
	DO 40 I=2,NL
	  TVX=T(I)*(1.+Q(I)*EPSI-Q(I))
	  TVY=T(I-1)*(1.+Q(I-1)*EPSI-Q(I-1))
	  GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(P(I-1)-P(I))/PH(I)
	  CPN(I)=CPD*(1.-Q(I))+CPV*Q(I)
	  H(I)=T(I)*CPN(I)+GZ(I)
          LV(I)=LV0-CPVMCL*(T(I)-273.15)
          HM(I)=H(I)+LV(I)*Q(I)
C
C  ***  Find level of minimum moist static energy    ***
C
          IF(HM(I).LT.AHMIN.AND.HM(I).LT.HM(I-1))THEN
           AHMIN=HM(I)
           IHMIN=I
          END IF
   40   CONTINUE
C
C  ***  CHECK WHETHER LOWEST LEVEL TEMPERATURE AND SPECIFIC HUMIDITY   ***
C  ***                          ARE REASONABLE                         ***
C  ***      Skip convection if HM increases monotonically upward       ***
C
        IF(T(1).LT.250.0.OR.Q(1).LE.0.0.OR.P(IHMIN).LT.400.)THEN
	 IFLAG=0
	 RETURN
	END IF
C
C  ***     Find that model level below the level of minimum moist       ***
C  ***  static energy that has the maximum value of moist static energy ***
C
        AHMAX=0.0
        DO 42 I=1,IHMIN
         IF(HM(I).GT.(AHMAX+100.0))THEN
          NK=I
          AHMAX=HM(I)
         END IF
   42   CONTINUE
C
C   ***  CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT PARCEL ORIGIN LEVEL ***
C   ***       (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980)      ***
C
	IF(T(NK).LT.250.0.OR.Q(NK).LE.0.0)GOTO 800
	RH=Q(NK)/QS(NK)
	CHI=T(NK)/(1669.0-122.0*RH-T(NK))
	PLCL=P(NK)*(RH**CHI)
	IF(PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN
	 IFLAG=2
	 GOTO 800
	END IF
C
C   ***  CALCULATE FIRST LEVEL ABOVE LCL (=ICB)  ***
C
	ICB=NL-1
	DO 50 I=NK+1,NL
	 IF(P(I).LT.PLCL)THEN
	  ICB=MIN(ICB,I)
	 END IF
   50   CONTINUE
	IF(ICB.GE.(NL-1))THEN
	 IFLAG=3
	 GOTO 800
	END IF
C
C   ***  CALCULATE MAXIMUM ALLOWED VALUE OF THE CLOUD BASE MASS FLUX  ***
C
	CBMFMAX=10.*(PH(ICB)-PH(ICB+1))*GINV*DELTI
C
C   *** FIND TEMPERATURE UP THROUGH ICB AND TEST FOR INSTABILITY           ***
C
C   *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL      ***
C   ***  TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC             ***
C   ***                   LIQUID WATER CONTENT                             ***
C
	CALL TLIFT(P,T,Q,QS,GZ,ICB,NK,TVP,TP,CLW,ND,NL,1)
	DO 54 I=NK,ICB
	 TV(I)=T(I)*(1.+Q(I)*EPSI-Q(I))
	 TVP(I)=TVP(I)-TP(I)*Q(NK)
   54   CONTINUE
C
C   ***  If there was no convection at last time step and parcel    ***
C   ***       is stable at ICB then skip rest of calculation        ***
C
        IF(CBMF.EQ.0.0.AND.TVP(ICB).LE.(TV(ICB)-DTMAX))THEN
	 IFLAG=0
	 GOTO 800
	END IF
C
C   ***  IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY ***
C
	IF(IFLAG.NE.4)IFLAG=1
C
C   ***  FIND THE REST OF THE LIFTED PARCEL TEMPERATURES          ***
C
	CALL TLIFT(P,T,Q,QS,GZ,ICB,NK,TVP,TP,CLW,ND,NL,2)
C
C   ***  SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF   ***
C   ***          PRECIPITATION FALLING OUTSIDE OF CLOUD           ***
C   ***      THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)     ***
C                 
        DO 57 I=1,NK
         EP(I)=0.0
         SIGP(I)=SIGS
         CLW(I)=0.0
   57   CONTINUE
        DO 60 I=NK+1,NL
	 TCA=TP(I)-273.15
	 IF(TCA.GE.0.0)THEN
	  ELACRIT=ELCRIT
	 ELSE
	  ELACRIT=ELCRIT*(1.0-TCA/TLCRIT)
	 END IF
         ELACRIT=MAX(ELACRIT,0.0)
	 EP(I)=1.0-ELACRIT/MAX(CLW(I),1.0E-8)
	 EP(I)=MAX(EP(I),0.0)
	 EP(I)=MIN(EP(I),1.0)
	 SIGP(I)=SIGS
   60   CONTINUE
C
C   ***       CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL     ***
C   ***                    VIRTUAL TEMPERATURE                    ***
C
	DO 64 I=ICB+1,NL
	 TV(I)=T(I)*(1.+Q(I)*EPSI-Q(I))
         TVP(I)=TVP(I)*(1.-Q(NK)+EP(I)*CLW(I))
   64   CONTINUE
C
C   ***        NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS       ***
C
	DO 70 I=1,NL
	 HP(I)=H(I)
	 QP(I)=Q(I)
	 NENT(I)=0
	 WATER(I)=0.0
	 EVAP(I)=0.0
	 MP(I)=0.0
	 M(I)=0.0
	 LVCP(I)=LV(I)/CPN(I)
	 DO 70 J=1,NL
	  QENT(I,J)=Q(J)
	  ELIJ(I,J)=0.0
	  MENT(I,J)=0.0
	  SIJ(I,J)=0.0
   70   CONTINUE
C
C  ***  FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S      ***
C  ***          HIGHEST LEVEL OF NEUTRAL BUOYANCY                 ***
C  ***     AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)           ***
C
        CAPE=0.0
	CAPEM=0.0
	INB=ICB+1
        INB1=INB
        DO 82 I=ICB+1,NL-1
         BY=(TVP(I)-TV(I))*(PH(I)-PH(I+1))/P(I)
         BYP=(TVP(I+1)-TV(I+1))*(PH(I+1)-PH(I+2))/P(I+1)
         CAPE=CAPE+BY
         IF(BY.GE.0.0)INB1=I+1
         IF(CAPE.GT.0.0)THEN
          INB=I+1
          CAPEM=CAPE
         END IF
   82 	CONTINUE
   81 	CONTINUE
        CAPE=CAPEM+BYP
        DEFRAC=CAPEM-CAPE
        DEFRAC=MAX(DEFRAC,0.001)
        FRAC=-CAPE/DEFRAC
        FRAC=MIN(FRAC,1.0)
        FRAC=MAX(FRAC,0.0)
C
C   ***   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL   ***
C
	DO 95 I=ICB,INB
	 HP(I)=H(NK)+(LV(I)+(CPD-CPV)*T(I))*EP(I)*CLW(I)
   95   CONTINUE                  
C
C   ***  CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),  ***
c   ***                   AT EACH MODEL LEVEL                       ***
C
	CBMFSUM=0.0
	DBOSUM=0.0
C   
C   ***     INTERPOLATE DIFFERENCE BETWEEN LIFTED PARCEL AND      ***
C   ***  ENVIRONMENTAL TEMPERATURES TO LIFTED CONDENSATION LEVEL  ***
C	
        TVPPLCL=TVP(ICB-1)-RD*TVP(ICB-1)*(P(ICB-1)-PLCL)/
     1    (CPN(ICB-1)*P(ICB-1))
        TVAPLCL=TV(ICB)+(TVP(ICB)-TVP(ICB+1))*(PLCL-P(ICB))/
     1    (P(ICB)-P(ICB+1))
        DTPBL=0.0
        DO 96 I=NK,ICB-1
         DTPBL=DTPBL+(TVP(I)-TV(I))*(PH(I)-PH(I+1))
   96   CONTINUE
        DTPBL=DTPBL/(PH(NK)-PH(ICB))
        DTMIN=TVPPLCL-TVAPLCL+DTMAX+DTPBL
        DTMA=DTMIN
C
C   ***  ADJUST CLOUD BASE MASS FLUX   ***
C
      CBMFOLD=CBMF
      CBMF=(1.-DAMP)*CBMF+0.1*ALPHA*DTMA
      CBMF=MAX(CBMF,0.0)
C
C   *** If cloud base mass flux is zero, skip rest of calculation  ***
C
      IF(CBMF.EQ.0.0.AND.CBMFOLD.EQ.0.0)THEN
       GOTO 800
      END IF
      CBMFSUM=0.5*(CBMF+CBMFOLD)
C
C   ***  LIMIT CLOUD BASE MASS FLUX TO AVOID VIOLATING CFL  ***
C
      IF(CBMFSUM.GT.CBMFMAX)THEN
       CBMF=CBMF*(CBMFMAX/CBMFSUM)
       CBMFSUM=CBMFMAX
      END IF        
C
C   ***   CALCULATE RATES OF MIXING,  M(I)   ***
C
      M(ICB)=0.0
      DO 103 I=ICB+1,INB
       K=MIN(I,INB1)
       DBO=ABS(TV(K+1)-TVP(K+1)-TV(K-1)+TVP(K-1))+
     1  ENTP*0.04*(PH(K)-PH(K+1))
       DBOSUM=DBOSUM+DBO
       M(I)=CBMF*DBO
  103 CONTINUE
      DO 110 I=ICB+1,INB
       M(I)=M(I)/DBOSUM  
  110 CONTINUE     
C
C   ***  CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING  ***
C   ***     RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING     ***
C   ***                        FRACTION (SIJ)                          ***
C
	DO 170 I=ICB+1,INB
	 QTI=Q(NK)-EP(I)*CLW(I)
	 DO 160 J=ICB,INB
	  BF2=1.+LV(J)*LV(J)*QS(J)/(RV*T(J)*T(J)*CPD)
	  ANUM=H(J)-HP(I)+(CPV-CPD)*T(J)*(QTI-Q(J))
	  DENOM=H(I)-HP(I)+(CPD-CPV)*(Q(I)-QTI)*T(J)
	  DEI=DENOM
	  IF(ABS(DEI).LT.0.01)DEI=0.01
	  SIJ(I,J)=ANUM/DEI
	  SIJ(I,I)=1.0
	  ALTEM=SIJ(I,J)*Q(I)+(1.-SIJ(I,J))*QTI-QS(J)
	  ALTEM=ALTEM/BF2
	  CWAT=CLW(J)*(1.-EP(J))
	  STEMP=SIJ(I,J)
	  IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR.
     1      ALTEM.GT.CWAT).AND.J.GT.I)THEN
	   ANUM=ANUM-LV(J)*(QTI-QS(J)-CWAT*BF2)
	   DENOM=DENOM+LV(J)*(Q(I)-QTI)
	   IF(ABS(DENOM).LT.0.01)DENOM=0.01
	   SIJ(I,J)=ANUM/DENOM
	   ALTEM=SIJ(I,J)*Q(I)+(1.-SIJ(I,J))*QTI-QS(J)
	   ALTEM=ALTEM-(BF2-1.)*CWAT
	  END IF
	  IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN
	   QENT(I,J)=SIJ(I,J)*Q(I)+(1.-SIJ(I,J))*QTI
	   ELIJ(I,J)=ALTEM
	   ELIJ(I,J)=MAX(0.0,ELIJ(I,J))
	   MENT(I,J)=M(I)/(1.-SIJ(I,J))
	   NENT(I)=NENT(I)+1
	  END IF
	  SIJ(I,J)=MAX(0.0,SIJ(I,J))
	  SIJ(I,J)=MIN(1.0,SIJ(I,J))
  160    CONTINUE
C
C   ***   IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS  ***
C   ***   AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES  ***
C
	 IF(NENT(I).EQ.0)THEN
	  MENT(I,I)=M(I)
	  QENT(I,I)=Q(NK)-EP(I)*CLW(I)
	  ELIJ(I,I)=CLW(I)
	  SIJ(I,I)=1.0
	 END IF 
  170   CONTINUE
	SIJ(INB,INB)=1.0
C
C   ***  NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL  ***
C   ***              PROBABILITIES OF MIXING                     ***
C
	DO 200 I=ICB+1,INB
	IF(NENT(I).NE.0)THEN
	 QP1=Q(NK)-EP(I)*CLW(I)
	 ANUM=H(I)-HP(I)-LV(I)*(QP1-QS(I))
	 DENOM=H(I)-HP(I)+LV(I)*(Q(I)-QP1)
	 IF(ABS(DENOM).LT.0.01)DENOM=0.01
	 SCRIT=ANUM/DENOM
	 ALT=QP1-QS(I)+SCRIT*(Q(I)-QP1)
	 IF(SCRIT.LT.0.0.OR.ALT.LT.0.0)SCRIT=1.0
	 ASIJ=0.0
	 SMIN=1.0
	 DO 175 J=ICB,INB
	  IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN
	   IF(J.GT.I)THEN
	    SMID=MIN(SIJ(I,J),SCRIT)
	    SJMAX=SMID
	    SJMIN=SMID
	    IF(SMID.LT.SMIN.AND.SIJ(I,J+1).LT.SMID)THEN
	     SMIN=SMID
	     SJMAX=MIN(SIJ(I,J+1),SIJ(I,J),SCRIT)
	     SJMIN=MAX(SIJ(I,J-1),SIJ(I,J))
	     SJMIN=MIN(SJMIN,SCRIT)
	    END IF
	   ELSE
	    SJMAX=MAX(SIJ(I,J+1),SCRIT)
	    SMID=MAX(SIJ(I,J),SCRIT)
	    SJMIN=0.0
	    IF(J.GT.1)SJMIN=SIJ(I,J-1)
	    SJMIN=MAX(SJMIN,SCRIT)
	   END IF
	   DELP=ABS(SJMAX-SMID)
	   DELM=ABS(SJMIN-SMID)
	   ASIJ=ASIJ+(DELP+DELM)*(PH(J)-PH(J+1))
	   MENT(I,J)=MENT(I,J)*(DELP+DELM)*(PH(J)-PH(J+1))
	  END IF
  175    CONTINUE
	 ASIJ=MAX(1.0E-21,ASIJ)
	 ASIJ=1.0/ASIJ
	 DO 180 J=ICB,INB
	  MENT(I,J)=MENT(I,J)*ASIJ
  180    CONTINUE
	 BSUM=0.0
	 DO 190 J=ICB,INB
	  BSUM=BSUM+MENT(I,J)
  190    CONTINUE
	 IF(BSUM.LT.1.0E-18)THEN
	  NENT(I)=0
	  MENT(I,I)=M(I)
	  QENT(I,I)=Q(NK)-EP(I)*CLW(I)
	  ELIJ(I,I)=CLW(I)
	  SIJ(I,I)=1.0
	 END IF
	END IF
  200   CONTINUE
C
C   ***  CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING    ***
C   ***             DOWNDRAFT CALCULATION                      ***
C
	IF(EP(INB).LT.0.0001)GOTO 405
C
C   ***  INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER   ***
C   ***                AND CONDENSED WATER FLUX                    ***
C
	QP(INB)=Q(INB-1)
	QP(INB+1)=Q(INB)
	WT(INB+1)=0.0
	JTT=2
C
C    ***                    BEGIN DOWNDRAFT LOOP                    ***
C
	DO 400 I=INB,1,-1
C
C    ***              CALCULATE DETRAINED PRECIPITATION             ***
C
	WDTRAIN=G*EP(I)*M(I)*CLW(I)
	IF(I.GT.1)THEN
	 DO 320 J=1,I-1
	 AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I)
	 AWAT=MAX(0.0,AWAT)
  320    WDTRAIN=WDTRAIN+G*AWAT*MENT(J,I)
	END IF
C
C    ***    FIND RAIN WATER AND EVAPORATION USING PROVISIONAL   ***
C    ***              ESTIMATES OF QP(I)AND QP(I-1)             ***
C     
c
c   ***  Value of terminal velocity and coeffecient of evaporation for snow   ***
c 
        COEFF=COEFFS
        WT(I)=OMTSNOW
c      
c  ***  Value of terminal velocity and coeffecient of evaporation for rain   ***
c
        IF(T(I).GT.273.0)THEN
         COEFF=COEFFR
         WT(I)=OMTRAIN
        END IF
        QSM=0.5*(Q(I)+QP(I+1))
        AFAC=COEFF*PH(I)*(QS(I)-QSM)/(1.0E4+2.0E3*PH(I)*QS(I))
        AFAC=MAX(AFAC,0.0)
        SIGT=SIGP(I)
        SIGT=MAX(0.0,SIGT)
        SIGT=MIN(1.0,SIGT)
        B6=100.*(PH(I)-PH(I+1))*SIGT*AFAC/WT(I)
        C6=(WATER(I+1)*WT(I+1)+WDTRAIN/SIGD)/WT(I)
        REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6))
        EVAP(I)=SIGT*AFAC*REVAP
        WATER(I)=REVAP*REVAP
C
C    ***  CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER     ***
C    ***              HYDROSTATIC APPROXIMATION                 ***
C   
	IF(I.EQ.1)GOTO 360
	DHDP=(H(I)-H(I-1))/(P(I-1)-P(I))
	DHDP=MAX(DHDP,10.0)
	MP(I)=100.*GINV*LV(I)*SIGD*EVAP(I)/DHDP
	MP(I)=MAX(MP(I),0.0)
C
C   ***   ADD SMALL AMOUNT OF INERTIA TO DOWNDRAFT              ***
C
	FAC=20.0/(PH(I-1)-PH(I))
	MP(I)=(FAC*MP(I+1)+MP(I))/(1.+FAC)
C
C   ***         LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION      ***
C
	AMPMAX=2.0*(PH(I)-PH(I+1))*DELTI
	AMP2=2.0*(PH(I-1)-PH(I))*DELTI
	AMPMAX=MIN(AMPMAX,AMP2)
	MP(I)=MIN(MP(I),AMPMAX)
C   
C    ***      FORCE MP TO DECREASE LINEARLY TO ZERO                 ***
C    ***      BETWEEN ABOUT 950 MB AND THE SURFACE                  ***
C
	  IF(P(I).GT.(0.949*P(1)))THEN
	   JTT=MAX(JTT,I)
	   MP(I)=MP(JTT)*(P(1)-P(I))/(P(1)-P(JTT))
	  END IF              
  360   CONTINUE
C
C    ***       FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT     ***
C
	IF(I.EQ.INB)GOTO 400
	IF(I.EQ.1)THEN
	 QTM=Q(1) 
	 QSTM=QS(1)
	ELSE
	 QTM=Q(I-1)
	 QSTM=QS(I-1)
	END IF
	QP(I)=QTM    
	IF(MP(I).GT.MP(I+1))THEN
	  RAT=MP(I+1)/MP(I)
	  QP(I)=QP(I+1)*RAT+Q(I)*(1.0-RAT)+100.*GINV*
     1       SIGD*(PH(I)-PH(I+1))*(EVAP(I)/MP(I))
	 ELSE
	  IF(MP(I+1).GT.0.0)THEN
	    QP(I)=(GZ(I+1)-GZ(I)+QP(I+1)*(LV(I+1)+T(I+1)*(
     1        CL-CPD))+CPD*(T(I+1)-T(I)))/(LV(I)+T(I)*(CL-CPD))
	  END IF
	END IF
	QP(I)=MIN(QP(I),QSTM)
	QP(I)=MAX(QP(I),0.0)
  400   CONTINUE
C
C   ***  CALCULATE SURFACE PRECIPITATION IN MM/DAY     ***
C
	PRECIP=PRECIP+WT(1)*SIGD*WATER(1)*3600.*24000./(ROWL*G)
C
  405   CONTINUE
C
C   ***  CALCULATE DOWNDRAFT VELOCITY SCALE AND SURFACE TEMPERATURE AND  ***
c   ***                    WATER VAPOR FLUCTUATIONS                      ***
C
      WD=BETA*ABS(MP(ICB))*0.01*RD*T(ICB)/(SIGD*P(ICB))
      QPRIME=0.5*(QP(1)-Q(1))
      TPRIME=LV0*QPRIME/CPD
C
C   ***  CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE  ***
C   ***                      AND MIXING RATIO                        ***
C
	DPINV=0.01/(PH(1)-PH(2))
	AM=0.0
	IF(NK.EQ.1)THEN
	 DO 410 K=2,INB
  410    AM=AM+M(K)
	END IF
	IF((2.*G*DPINV*AM).GE.DELTI)IFLAG=4
	FT(1)=FT(1)+G*DPINV*AM*(T(2)-T(1)+(GZ(2)-GZ(1))/CPN(1))
	FT(1)=FT(1)-LVCP(1)*SIGD*EVAP(1)
	FT(1)=FT(1)+SIGD*WT(2)*(CL-CPD)*WATER(2)*(T(2)-
     1   T(1))*DPINV/CPN(1)
	FQ(1)=FQ(1)+G*MP(2)*(QP(2)-Q(1))*
     1    DPINV+SIGD*EVAP(1)
	FQDET(1)=FQ(1)
	FQ(1)=FQ(1)+G*AM*(Q(2)-Q(1))*DPINV
	FQSUB(1)=G*AM*(Q(2)-Q(1))*DPINV
	AMDE=0.0
	DO 415 J=2,INB
	 FQ(1)=FQ(1)+G*DPINV*MENT(J,1)*(QENT(J,1)-Q(1))
	 FQSUB(1)=FQSUB(1)+G*DPINV*MENT(J,1)*(QENT(J,1)-Q(1))
  415   CONTINUE
C
C   ***  CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO  ***
C   ***               AT LEVELS ABOVE THE LOWEST LEVEL                   ***
C
C   ***  FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES  ***
C   ***                      THROUGH EACH LEVEL                          ***
C
	DO 500 I=2,INB
	DPINV=0.01/(PH(I)-PH(I+1))
	CPINV=1.0/CPN(I)
	AMP1=0.0
	AD=0.0
	IF(I.GE.NK)THEN
	 DO 440 K=I+1,INB+1
  440    AMP1=AMP1+M(K)
	END IF
	DO 450 K=1,I
	DO 450 J=I+1,INB+1
	 AMP1=AMP1+MENT(K,J)
	 MUP(I)=AMP1
  450   CONTINUE
	IF((2.*G*DPINV*AMP1).GE.DELTI)IFLAG=4
	DO 470 K=1,I-1
	DO 470 J=I,INB
	AD=AD+MENT(J,K)
	MDOWN(I)=AD
  470   CONTINUE
	FT(I)=FT(I)+G*DPINV*(AMP1*(T(I+1)-T(I)+(GZ(I+1)-GZ(I))*
     1   CPINV)-AD*(T(I)-T(I-1)+(GZ(I)-GZ(I-1))*CPINV))
     2   -SIGD*LVCP(I)*EVAP(I)
	FT(I)=FT(I)+G*DPINV*MENT(I,I)*(HP(I)-H(I)+
     1    T(I)*(CPV-CPD)*(Q(I)-QENT(I,I)))*CPINV
	FT(I)=FT(I)+SIGD*WT(I+1)*(CL-CPD)*WATER(I+1)*
     1    (T(I+1)-T(I))*DPINV*CPINV
	FQ(I)=FQ(I)+G*DPINV*(AMP1*(Q(I+1)-Q(I))-
     1    AD*(Q(I)-Q(I-1)))
	FQSUB(I)=FQ(I)
	DO 480 K=1,I-1
	 AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I)
	 AWAT=MAX(AWAT,0.0)
	 FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-Q(I))
	 FQDET(I)=G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-Q(I))
  480   CONTINUE
	DO 490 K=I,INB
	 FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-Q(I))
	 FQDET(I)=FQDET(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-Q(I))
  490   CONTINUE
	FQ(I)=FQ(I)+SIGD*EVAP(I)+G*(MP(I+1)*
     1    (QP(I+1)-Q(I))-MP(I)*(QP(I)-Q(I-1)))*DPINV
	FQDET(I)=FQDET(I)+SIGD*EVAP(I)+G*(MP(I+1)*
     1    (QP(I+1)-Q(I))-MP(I)*(QP(I)-Q(I-1)))*DPINV
  500   CONTINUE
C
C   *** Adjust tendencies at top of convection layer to reflect  ***
C   ***       actual position of the level zero CAPE             ***
C
        FQOLD=FQ(INB)
        FQ(INB)=FQ(INB)*(1.-FRAC)
        FQ(INB-1)=FQ(INB-1)+FRAC*FQOLD*((PH(INB)-PH(INB+1))/
     1   (PH(INB-1)-PH(INB)))*LV(INB)/LV(INB-1)
        FTOLD=FT(INB)
        FT(INB)=FT(INB)*(1.-FRAC)
        FT(INB-1)=FT(INB-1)+FRAC*FTOLD*((PH(INB)-PH(INB+1))/
     1   (PH(INB-1)-PH(INB)))*CPN(INB)/CPN(INB-1)
C
C   ***   Very slightly adjust temperature tendencies to force exact   ***
C   ***                     enthalpy conservation                      ***
C
	ENTS=0.0
	DO 680 I=1,INB
	 ENTS=ENTS+(CPN(I)*FT(I)+LV(I)*FQ(I))*(PH(I)-PH(I+1))	
  680	CONTINUE
	ENTS=ENTS/(PH(1)-PH(INB+1))
	DO 640 I=1,INB
	 FT(I)=FT(I)-ENTS/CPN(I)
  640	CONTINUE
	DO 750 I=1,INB
	 DO 710 K=1,INB
	  DET(I)=DET(I)+MENT(K,I)
  710    CONTINUE
	 DO 720 K=1,INB
	  ENT(I)=ENT(I)+MENT(I,K)
  720    CONTINUE
	 ENT(I)=ENT(I)-MENT(I,I)
  750   CONTINUE
C
C   ***           RETURN           ***
C
	RETURN
  800   CONTINUE
	RETURN
	END
C
C ---------------------------------------------------------------------------
C
	SUBROUTINE TLIFT(P,T,Q,QS,GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)
	REAL GZ(ND),TPK(ND),CLW(ND),P(ND)
	REAL T(ND),Q(ND),QS(ND),TVP(ND),LV0
C
C   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
C
	CPD=1005.7
	CPV=1870.0
	CL=4190.0
        RV=461.5
	RD=287.04
	LV0=2.501E6
C
	CPVMCL=CL-CPV
	EPS=RD/RV
	EPSI=1./EPS
C
C   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
C
	AH0=(CPD*(1.-Q(NK))+CL*Q(NK))*T(NK)+Q(NK)*(LV0-CPVMCL*(
     1   T(NK)-273.15))+GZ(NK)
	CPP=CPD*(1.-Q(NK))+Q(NK)*CPV
	CPINV=1./CPP
C
	IF(KK.EQ.1)THEN
C
C   ***   CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE   ***
C
	DO 50 I=1,ICB-1
	 CLW(I)=0.0
   50   CONTINUE
	DO 100 I=NK,ICB-1
	 TPK(I)=T(NK)-(GZ(I)-GZ(NK))*CPINV
	 TVP(I)=TPK(I)*(1.+Q(NK)*EPSI)
  100   CONTINUE
	END IF
C
C    ***  FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE    ***
C
	NST=ICB
	NSB=ICB
	IF(KK.EQ.2)THEN  
	 NST=NL
	 NSB=ICB+1
	END IF
	DO 300 I=NSB,NST
	 TG=T(I)
	 QG=QS(I)
	 ALV=LV0-CPVMCL*(T(I)-273.15)
	 DO 200 J=1,2
	  S=CPD+ALV*ALV*QG/(RV*T(I)*T(I))
	  S=1./S
	  AHG=CPD*TG+(CL-CPD)*Q(NK)*T(I)+ALV*QG+GZ(I)
	  TG=TG+S*(AH0-AHG)
	  TC=TG-273.15
	  DENOM=243.5+TC
	  DENOM=MAX(DENOM,1.0)
          IF(TC.GE.0.0)THEN  
	   ES=6.112*EXP(17.67*TC/DENOM)
          ELSE  
           ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG))
          END IF  
	  QG=EPS*ES/(P(I)-ES*(1.-EPS))
  200    CONTINUE
	 ALV=LV0-CPVMCL*(T(I)-273.15)
         TPK(I)=(AH0-(CL-CPD)*Q(NK)*T(I)-GZ(I)-ALV*QG)/CPD
	 CLW(I)=Q(NK)-QG
	 CLW(I)=MAX(0.0,CLW(I))
	 RG=QG/(1.-Q(NK))
	 TVP(I)=TPK(I)*(1.+RG*EPSI)
  300   CONTINUE
	RETURN
	END
        SUBROUTINE PLREG(X,Y,Z,NA,NB,NX,NY,XSCA,YSCA,XLAB,YLAB,ZLAB)
c
c   Subroutine PLREG contours the array Z(NA,NB) as a function of X(NA) and
c   Y(NB), using NX points in the X direction and NY points in the Y
c   direction. The aspect ratio of the plot can be set by changing the
c   scaling of the X axis, XSCA, and/or the Y axis, YSCA. (0<XSCA,YSCA<1).
c   The X and Y axes are labelled by XLAB and YLAB (ideal dimensions
c   CHARACTER*12), and the contoured quantity is called ZLAB.
c         NUMC is the number of contours.
c
        REAL X(NA), Y(NB), Z(NA,NB)
        CHARACTER*(*) XLAB,YLAB,ZLAB
        CHARACTER*12 CHXMIN,CHXMAX,CHYMIN,CHYMAX
        CHARACTER*50 CHZ
c
        NUMC=10
c
        XSCA1=MIN(XSCA,1.0)
        XSCA1=MAX(XSCA,0.2)
        YSCA1=MIN(YSCA,1.0)
        YSCA1=MAX(YSCA,0.2)
        XSCA1=0.8*XSCA
        YSCA1=0.8*YSCA
        XMAX=X(NX)
        YMAX=Y(NY)
        XMIN=X(1)
        YMIN=Y(1)
        ZMAX=Z(1,1)
        ZMIN=Z(1,1)
        DO 100 I=1,NX
         DO 100 J=1,NY
          ZMAX=MAX(ZMAX,Z(I,J))
          ZMIN=MIN(ZMIN,Z(I,J))
  100   CONTINUE
        DELZ=MAX((ZMAX-ZMIN),1.0E-16)
        AINC=DELZ/FLOAT(NUMC)
        WRITE(CHZ,210)ZLAB,ZMIN,ZMAX
  210   FORMAT(1X,A14,'  From  ',F10.4,'  to  ',F10.4)
        WRITE(CHXMIN,220)XMIN
  220   FORMAT(1X,F8.2)
        WRITE(CHXMAX,220)XMAX
        WRITE(CHYMIN,220)YMIN
        WRITE(CHYMAX,220)YMAX
        SXMAX=0.5*(1.+XSCA1)+0.05
        SXMIN=0.5*(1.-XSCA1)+0.05
        SXMIN2=SXMIN+0.02
        SYMAX=0.5*(1.+YSCA1)
        SYMIN=0.5*(1.-YSCA1)
        SYMAX2=SYMAX+0.05
        SYMIN2=SYMIN-0.05
        CALL GSELNT(0)
        CALL SET(SXMIN,SXMAX,SYMIN,SYMAX,0.0,1.,0.,1.,1)
        CALL PERIM(NX,1,NY,1)
        CALL CONREC(Z,NA,NX,NY,0.,0.,AINC,-1,0,-750)
        CALL SET(0.,1.,0.,1.,0.,1.,0.,1.,1)
        CALL GSELNT(0)
        CALL PLCHHQ(0.55,SYMAX2,CHZ,0.012,0.,0.)
        CALL PLCHHQ(SXMIN2,SYMIN,CHYMIN,0.012,0.,1.)
        CALL PLCHHQ(SXMIN2,0.5,YLAB,0.012,0.,1.)
        CALL PLCHHQ(SXMIN2,SYMAX,CHYMAX,0.012,0.,1.)
        CALL PLCHHQ(SXMIN,SYMIN2,CHXMIN,0.012,0.,0.)
        CALL PLCHHQ(SXMAX,SYMIN2,CHXMAX,0.012,0.,0.)
        CALL PLCHHQ(0.55,SYMIN2,XLAB,0.012,0.,0.)
        CALL FRAME
        RETURN
        END

