	PROGRAM STRATUS
C
C   ***  STRATUS solves equations (13.3.27), (13.3.28), (13.3.36),      ***
C   ***       (13.3.38) and (13.3.39) to find the properties            ***
C   ***      of stratocumulus-topped equilibrium mixed layers           ***
C   ***    Input is from the terminal and output is to the file         ***
C   ***                       stratus.out                               ***
C
	REAL HB,HVB,HT,HTS,HV0,HVT,HVM
	OPEN(UNIT=8,FILE='stratus.out',STATUS='NEW')
C
C   ***               Define thermodynamic constants                    ***
C
	CPD=1005.7
	CPV=1870.0
	CL=4190.0
	RD=287.04
	RV=461.5
	EP=0.622
	ALV=2.5E6
	G=9.8
C   
C   ***             Input paramters from terminal                      ***
C
	PRINT *, 'Enter altitude H (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 radiative cooling rate (degrees C per day)'
	READ *, QC
	PRINT *, 'Enter surface wind speed (m/s)'
	READ *, V0
	PRINT *, 'Enter radiative flux divergence at cloud top (W/m**2)'
	READ *, DELI
	PRINT *, 'Enter value of mixing constant k'
	READ *, AK
	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)V0
	WRITE(8,6)DELI
	WRITE(8,7)AK
    1   FORMAT(10X,'                         H =  ',F7.1,'  meters')
    2   FORMAT(10X,'  Virtual temperature at H =  ',F7.1,'  degrees C')
    3   FORMAT(10X,'         Mixing ratio at H =  ',F7.1,'  g/kg')
    4   FORMAT(10X,'    Radiative cooling rate =  ',F7.1,'  degrees ',
     1    'C/day')
    5   FORMAT(10X,'        Surface wind speed =  ',F7.1,'  m/s')
    6   FORMAT(10X,'Cloudtop infrared flux div.=  ',F7.1,'  W/m**2')
    7   FORMAT(10X,'           Mixing constant =  ',F7.1,/) 
C
C   ***   Normalize parameters, find T_m, find approximate    ***
C   ***              saturation mixing ratio at H             ***
C
	RT=0.001*RT
	TVM=TVM+273.15
	HVM=CPD*TVM+G*ZM             
	TM=TVM*(1.+RT)/(1.+RT/0.622)
	PM=1000.0-0.0975*ZM
	CALL SAT(TM,PM,RSM)
	RAD=CPD*1.157E-5*QC
	ROW=1.0      
	DELI=DELI/ROW
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',///)
	WRITE(8,10)
   10	FORMAT( 20X, '(ALL STATIC ENERGIES ARE DIVIDED BY c_p)',//)
	WRITE(8,11)
   11	FORMAT(2X, '     ','Surface  ',' Cloud   ','Cloud')
	WRITE(8,12)
   12	FORMAT(2X,' SST  ','  RH     ',' Base   ',' Top   ','h_lvb  ',
     1   'h_lvt  ',' h_b   ',' h_t   ','h_t^*  ','h_lv0-h_lvb')
	WRITE(8,13)
   13	FORMAT(2X,' (C)  ',' (%)     ',' (m)    ',' (m)   ',' (K)   ',
     1   ' (K)   ',' (K)   ',' (K)   ',' (K)   ','    (K)')
  	WRITE(8,14)
   14	FORMAT(2X,'----  ','------- ','------  ','------  ','-----  ',
     1   '-----  ','-----  ','-----  ','-----  ','-----------')
C  
C  ***         Specify initial paramters              ***
C
	TT=TM
	ZT=ZM
	RSI=RSM
	ISST1=SST1
	ISST2=SST2
C
C  ***      Loop through sea surface temperatures            ***
C
	DO 100 I=ISST1,ISST2
	SST=FLOAT(I)+273.15
	CALL SAT(SST,1000.0,R0)
	HV0=CPD*SST*(1.+R0/0.622)/(1.+R0)
C
C  ***     Outer loop of nested iteration to solve equation  ***
C  ***                 for boundary layer top                ***
C
	DO 20 J=1,3
	GMGD=(1.+ALV*RSI/(RD*TT))/(1.+ALV*ALV*RSI/(RV*TT*TT*CPD))
	GMGD=GMGD/(1.+RSI)
	GM=GMGD*G/CPD
	BETA=(CPV-CPD+G/GM)*TT/(1.+RSI)     
	A1=AK*GMGD/(1.+AK*GMGD*(1.+R0))
	A=A1*DELI/(ZM*RAD)
	B=A1*(ALV-BETA)*(R0-RT)/DELH-(HVM-HV0)/DELH
     1    -DELI/(W0*DELH*(1.+AK*GMGD*(1.+R0)))
	B=B+A1*R0*(CL-CPV)*(SST-TT)/DELH
	X=1
C
C  ***               Inner loop                     ***
C
	DO 15 K=1,10
	 X=0.5*(X+EXP(B+A/X))
	 X=MIN(10.0,X)
   15	CONTINUE
   25	CONTINUE
	ZT=X*ZM            
	WT=X*WM
	HVT=HVM+DELH*LOG(X)
	TVT=(HVT-G*ZT)/CPD
        TT=TVT*(1.+RT)/(1.+RT/0.622)
C
C  ***    Calculation of saturation mixing ratio is approximate   ***
C
	PT=1000.0-0.0975*ZT
	CALL SAT(TT,PT,RSI)
   20   CONTINUE
	RES=ZM*(X-EXP(B+A/X))
C
C  ***              Find all other variables           ***
C
	RB=(W0*R0+WT*RT)/(W0+WT)
	HVB=HVT+(W0*(HV0-HVT)-DELI)/(W0+WT)
	TVB=HVB/CPD
	TB=TVB*(1.+RB)/(1.+RB/0.622)
	CALL SAT(TB,1000.0,RBS)
	RH=100.0*RB/RBS
	TC=2840.0/(3.5*LOG(TB)-LOG(RB/(0.622+RB))-11.71)+55.0
	TVC=TC*(1.+RB/0.622)/(1.+RB)
	ZC=CPD*(TVB-TVC)/G
	RRM=MIN(RB,RBS)
	HB=TVB+ALV*RRM/CPD
	HTS=TVT+(ALV*RSI+G*ZT)/CPD
	HT=TVT+(ALV*RT+G*ZT)/CPD
	SSTC=SST-273.15
	HVB=HVB/CPD
	HVT=HVT/CPD
	DELVI=HV0/CPD-HVB
	ZT=MIN(ZT,9999.9)
C
C  ***   Write output variables to the file stratus.out      ***
C
	WRITE(8,30)SSTC,RH,ZC,ZT,HVB,HVT,HB,HT,HTS,DELVI
   30   FORMAT(X,F5.1,2X,F5.1,3X,F6.1,2X,F6.1,2X,F5.1,2X,
     1   F5.1,2X,F5.1,2X,F5.1,2X,F5.1,2X,X,F6.2)
C	WRITE(8,40)RES
C   40	FORMAT(20X,'Residual = ',F6.1,' meters')
  100	CONTINUE
	CLOSE(UNIT=8,STATUS='SAVE')
	STOP
	END
C
C  ***    Subroutine for calculating the saturation mixing ratio  ***
C
	SUBROUTINE SAT(T,P,R)
	TC=T-273.15
	ES=6.112*EXP(17.67*TC/(243.5+TC))
	R=0.622*ES/(P-ES)
	RETURN
	END
