	PROGRAM TRADE
C
C  ***    This program solves equations (13.4.12)-(13.4.15),    ***
C  ***     (13.4.21), and (13.4.28)-(13.4.33) for the           ***
C  ***       properties of equilibrium trade cumulus            ***
C  ***                   boundary layers                        ***
C
	PARAMETER (NP=60)
	REAL AR0(NP), ARB(NP), ARC(NP), ARI(NP), ARSI(NP), ARH(NP)
	REAL AAMU(NP), AAMD(NP)
	REAL HT,HB,HLVM,HTSM,HLVT,HTSI,HTST,HLV0P,HLVB,HLVI,HLV0
	INTEGER IFLAG(NP)
	OPEN(UNIT=8,FILE='trade.out',STATUS='NEW')
C
C  ***             Define thermodynamic constants               ***
C
	CPD=1005.7
	CPV=1870.0
	CL=4190.0
	RD=287.04
	RV=461.5
	ALV=2.5E6
	G=9.8
C
C  ***              Read in parameters from terminal            ***
C
	PRINT *, 'Enter altitude H (in meters)'
	READ *, ZM
	PRINT *, 'Enter virtual temperature at altitude H (degrees C)'
	READ *, TVM
	PRINT *, 'Enter mixing ratio at altitude H (g/kg)'
	READ *, RT
	PRINT *, 'Enter free atmosphere cooling rate (degrees C/day)'
	READ *, QC
	PRINT *, 'Enter boundary layer cooling rate (degrees C/day)'
	READ *, QCB
	PRINT *, 'Enter surface wind speed (m/s)'
	READ *, V0
	PRINT *, 'Enter lowest sea surface temperature to consider (C)'
	READ *, SST1
	PRINT *, 'Enter highest sea surface temperature to consider (C)'
	READ *, SST2
C
C   ***       Write out values of input parameters        ***
C
	WRITE(8,1)ZM
	WRITE(8,2)TVM
	WRITE(8,3)RT
	WRITE(8,4)QC
	WRITE(8,5)QCB
	WRITE(8,6)V0
    1   FORMAT(10X,'                           H =  ',F7.1,'  meters')
    2   FORMAT(10X,'    Virtual temperature at H =  ',F7.1,'  degrees'
     1   ' C')
    3   FORMAT(10X,'           Mixing ratio at H =  ',F7.1,'  g/kg')
    4   FORMAT(10X,'Free atmosphere cooling rate =  ',F7.1,'  degrees ',
     1    'C/day')
    5   FORMAT(10X,' Boundary layer cooling rate =  ',F7.1,'  degrees ',
     1    'C/day')
    6   FORMAT(10X,'          Surface wind speed =  ',F7.1,'  m/s')
C
C   ***   Normalize parameters, find T_m, find approximate    ***
C   ***              saturation mixing ratio at H             ***
C
	RT=0.001*RT
	TVM=TVM+273.15
	HLVM=CPD*TVM+G*ZM             
	TM=TVM*(1.+RT)/(1.+RT/0.622)
	PM=1000.0-0.0975*ZM
	CALL SAT(TM,PM,RSM)
	HTSM=TM+(ALV*RSM+G*ZM)/CPD
	RAD=CPD*1.157E-5*QC
	RADB=CPD*1.157E-5*QCB
C
C  ***       Find lapse rate of virtual static energy        ***
C
	DTDZ=-(G/CPD)*(1.+RSM*ALV/(RD*TM))/
     1   (1.+ALV*ALV*RSM/(RV*CPD*TM*TM))
	DHDZ=DTDZ*CPD*(1.+0.608*RSM*(1.+ALV/(RV*TM)))+
     1   G*(1.+0.608*RSM*CPD/RD)
	DELH=ZM*DHDZ
C
C  ***         Find and print vertical velocity at H and     ***
C  ***                  surface flux velocity                ***
C
	W0=1.2E-3*V0
	WM=RAD/DHDZ             
	WRITE(8,8) (100.0*W0)
	WRITE(8,9) (100.0*WM)
    8   FORMAT(10X,'        Surface flux velocity = ',F7.1,'  cm/s')
    9	FORMAT(10X,'       Vertical velocity at H = ',F7.1,'  cm/s',///)
C
C  ***              Write output headers                     ***
C
	WRITE(8,10)
   10	FORMAT(20X, '(ALL STATIC ENERGIES ARE DIVIDED BY c_p)',//)
	WRITE(8,11)
   11	FORMAT(1X,'   ','  Cloud', '  Inv. ', '  Cloud')
	WRITE(8,12)
   12	FORMAT(1X,'SST','  Base ', '  Alt. ', '   Top ',' h_lvb',' h_lvi',
     1   ' h_lvt','  h_b ',' h_ts ')
	WRITE(8,13)
   13	FORMAT(1X,'(C)','   (m) ', '  (m)  ', '   (m) ','  (K) ','  (K) ',
     1   '  (K) ','  (K) ','  (K) ')
	WRITE(8,14)
   14	FORMAT(1X,'---','  ---- ', '  ---- ', '  -----',' -----',' -----',
     1   ' -----','  --- ',' ---- ')
C
C  ***      Loop through sea surface temperatures            ***
C
	ISST1=SST1
	ISST2=SST2
	DO 100 I=ISST1,ISST2
	SST=FLOAT(I)+273.15
	IFLAG(I)=1
	N=0
	CALL SAT(SST,1000.0,R0)
	HLV0=CPD*SST*(1.+R0/0.622)/(1.+R0)
	AR0(I)=R0
C
C  ***    Set provisional value of ratio z_m/z_t (=X)    ***
C
	X=1.0
	DM=1000.0
   25	CONTINUE
	IF(X.LT.0.1)GOTO 28 
	ZT=X*ZM            
	WT=X*WM
	HLVT=HLVM+ZM*DHDZ*LOG(X)
C
C  ***    Find approximate value of saturation mixing ratio at cloud top  ***
C  ***              and set provisional values of ZC and DELZI            ***
C
	TVT=(HLVT-G*ZT)/CPD
	TT=TVT*(1.+RT)/(1.+RT/0.622)
	PT=1000.0-0.0975*ZT
	CALL SAT(TT,PT,RST)
	ZC=300.0
	DELZI=500.0
	DELZI=MIN(DELZI,(1.+ZT-ZC))
C
C  ***      This loop iterates to find values of DELZI and ZC given X    ***
C  ***              and solves equations (13.4.12) - (13.4.15),          ***      
C  ***                      and (13.4.30) - (13.4.33)                    ***
C
	DO 20 J=1,21 
	HLVB=HLV0-RADB*ZC/W0
	HLVI=2.*HLVT-HLVB-2.*RADB*(ZT-ZC)/WT
	IF(HLVT.EQ.HLVI)HLVI=HLVT-1.0
	AMD=WT-RADB*DELZI/(HLVT-HLVI)
	IF(HLVI.EQ.HLVB)HLVI=HLVB+1.0
	AMU=RADB*(ZT-DELZI-ZC)/(HLVI-HLVB)-0.5*WT
	AMP=MAX(AMU,0.0001)
	WE=AMU+0.5*WT
	A=AMU*W0*(WE+0.5*WT)
	B=WT*(WE*(W0+AMU)+0.5*WT*(WE+W0))
	RI=(A*R0+B*RT)/(A+B)
	RC=(W0*WE*R0+WT*(WE+W0)*(RT-0.5*RI))/(WE*W0+0.5*
     1    WT*(WE+W0))
	RB=(W0*R0+WE*RC)/(W0+WE)
	RBT=RB
	RB=MAX(RB,0.001)
	TVB=HLVB/CPD
	TB=TVB*(1.+RB)/(1.+RB/0.622)
	TB=MAX(TB,270.0)  
	CALL SAT(TB,1000.0,RBS)
	RB=MIN(RB,RBS)
	RH=RB/RBS
C
C  **      (Find lifted condensation level and temperature)   ***
C
	TC=2840.0/(3.5*LOG(TB)-LOG(RB/(0.622+RB))-11.71)+55.0
	TVC=TC*(1.+RB/0.622)/(1.+RB)
	ZCT=CPD*(TVB-TVC)/G
	ZC=0.5*(ZC+ZCT)
	ZC=MAX(ZC,0.0)
	ZC=MIN(ZC,700.0)
	SIG=(RST-RT)/(RB-RT+ALV*RST*(HLVT-HLVB)/(RV*TT*HLVT))
	DELZIT=0.5*SIG*(ZT-ZC)
	DELZI=0.5*(DELZI+DELZIT)
	DELZI=MAX(DELZI,10.0)
	DELZI=MIN(DELZI,(0.5*ZT))
   20	CONTINUE
C
C  ***    Now see how well (13.4.21) is satisfied, if it is,  ***
C  ***       remember values of quantities. In either         ***
C  ***      case, decrease value of X and keep going.         ***
C
	RRM=MIN(RB,RBS)
	HB=TVB+ALV*RRM/CPD
	TVI=(HLVI-G*(ZT-DELZI))/CPD
	TI=TVI*(1.+RI)/(1.+RI/0.622)
	PI=1000.0-0.0975*(ZT-DELZI)
	CALL SAT(TI,PI,RSI)
	HTSI=TVI+(ALV*RSI+G*(ZT-DELZI))/CPD
	HTST=TVT+(ALV*RST+G*ZT)/CPD
	AL=(ZT-DELZI-ZC)/DELZI
	HTS=(HTST+AL*HTSI)/(1.+AL)
	DIF=ABS(HB-HTS)
	X=X-0.001
	IF(DIF.LT.DM.AND.AMU.GT.0.0.AND.AMD.GT.0.0)THEN
	DM=DIF
	N=1
	AZT=ZT
	AZC=ZC
	ATT=TVT
	AHB=HB
	AHTS=HTST
	AHLVI=HLVI
	AHLVB=HLVB
	ARB(I)=RBT*1000.0
	ARC(I)=RC*1000.0
	ARI(I)=RI*1000.0
	ARH(I)=RH*100.0
	ARSI(I)=RSI*1000.0
	AHLVT=HLVT
	ADELZI=DELZI
        AAMU(I)=AMU*100.0
	AAMU(I)=MIN(AAMU(I),99.9)
	AAMD(I)=AMD*100.0
	AAMD(I)=MIN(AAMD(I),99.9)
	AAP=AMP
	END IF
	GOTO 25
   28	CONTINUE
C
C  ***    Write out first set of quantities   ***
C
	IF(DM.GT.1.0.OR.N.EQ.0)THEN    
	 IFLAG(I)=0
	 GOTO 100
	END IF
	HT=ATT+(ALV*RT+G*AZT)/CPD
	HTS=AHTS
	HB=AHB
	AZI=AZT-ADELZI
	JSSTC=SST-273.15
	HLVB=AHLVB/CPD
	HLVT=AHLVT/CPD
	HLVI=AHLVI/CPD
	HLV0P=HLV0/CPD
	WRITE(8,30)JSSTC,AZC,AZI,AZT,HLVB,HLVI,HLVT,HB,HTS
   30	FORMAT(1X,I3,1X,F6.1,1X,F6.1,1X,F6.1,1X,F5.1,1X,F5.1,         
     1   1X,F5.1,1X,F5.1,1X,F5.1)               
  100	CONTINUE
C
C  ***   Print other quantities   ***
C
	WRITE(8,110)
  110	FORMAT(2X,///)
	WRITE(8,111)
  111	FORMAT(1X,'SST','  RH  ','   r_b ','   r_c ','   r_i ','  rs_i ',
     1   '   M_u ','   M_d ')
	WRITE(8,112)
  112	FORMAT(1X,'(C)','  (%) ',' (g/kg)',' (g/kg)',' (g/kg)',' (g/kg)',
     1   ' (cm/s)',' (cm/s)')
	WRITE(8,113)
  113	FORMAT(1X,'---',' ---- ',' ------',' ------',' ------',' ------',
     1   ' ------',' ------')
	DO 200 I=ISST1,ISST2
	IF(IFLAG(I).EQ.0)GOTO 200
	WRITE(8,115)I,ARH(I),ARB(I),ARC(I),ARI(I),ARSI(I),AAMU(I),
     1   AAMD(I)
  115	FORMAT(1X,I3,1X,F5.1,2X,F4.1,3X,F4.1,3X,F4.1,3X,F4.1,2X,
     1   F5.2,3X,F5.2)
  200   CONTINUE
	CLOSE(UNIT=8,STATUS='SAVE')
	STOP
	END
C
C  ***    Subroutine SAT calculates an approximate value of the ***        
C  ***                 saturation mixing ratio                  ***
C
	SUBROUTINE SAT(T,P,R)
	TC=T-273.15
	ES=6.112*EXP(17.67*TC/(243.5+TC))
	DENOM=MAX(1.0,(P-ES))
	R=0.622*ES/DENOM 
	RETURN
	END
