CCCCC   LA SMOG PROGRAM
CCCCC	MODIFICATIONS MADE NOV. 1994, NOV. 1995, AND OCT.1996
CCCCC   BY DAN JAFFE, GEOPHYSICAL INSTITUTE AND DEPARTMENT
CCCCC   OF CHEMISTRY, UNIVERSITY OF ALASKA FAIRBANKS
CCCCC	EMAIL:	 DJAFFE@GI.ALASKA.EDU
CCCCC
CCCCC   IMPORTANT NOTE:  THIS PROGRAM IS PURELY FOR TEACHING PURPOSES
CCCCC	IT SHOULD NOT BE USED FOR RESEARCH OR POLICY APPLICATIONS
CCCCC	AS-IS.	PROGRAM CALCULATES O3
CCCCC	CONCENTRATIONS IN A BOX USING SIMPLE HC CHEMISTRY.
CCCCC   INPUTS ARE NOX EMISSIONS, HC EMISSIONS, AND WIND SPEED
CCCCC	UNITS FOR NOX AND HC EMISSIONS ARE PPBV PER SECOND.
CCCCC   UNITS FOR WIND SPEED ARE M/S AND ARE UTILIZED ASSUMING A
CCCCC   10 x 10 KM BOX SIZE
CCCCC   ALL CALCS ARE IN UNITS OF MIXING RATIO SO FOR EX. 
CCCCC   1E-9 IS 1 PPBV
CCCCC     
	CHARACTER*20 DFILE
	DIMENSION R(50),X(50),A(50),B(50),XA(50),XS(50)
	WRITE(*,10)
 10     FORMAT('    ENTER Q_NO AND Q_NMHC (cr) ')
	WRITE(*,11)
 11	FORMAT('  UNITS ARE PPBV PER SECOND ')
	READ(*,*)QNO,QHC
	WRITE(*,*)QNO,QHC
	WRITE(*,12)
	QNO=QNO*1E-9
	QHC=QHC*1E-9
 12	FORMAT(' ENTER THE WIND SPEED (m/s)')
	READ(*,*)V
	WRITE(*,*)V
	IF(V.EQ.0.0)V=1.0E-10
	WRITE(*,15)
 15     FORMAT('   ENTER THE #DAYS ')
	READ(*,*)DAYS
	WRITE(*,*)DAYS
	DHR=0.25
	WRITE(*,18)
 18     FORMAT('   ENTER THE NAME OF THE OUTPUT FILE')
	READ(*,19)DFILE
 19     FORMAT(A)
	OPEN(88,FILE=DFILE,STATUS='NEW')
	QCO=10.0*QHC
	HHO=0.01
	NXS=19
	CONS=3.14159/6.0
	XN=(1.0*1000.0)/(0.0821*293.0)
	AHC=0.0
	AHO=0.0
	AHN=0.0
	AER=0.0
	TADV=10000.0/V
	RTADV=1.0/TADV
	THICKJ=0.0
	THINJ=0.0
	FLUX=0.1
	DO 20 NNN=1,NXS
	XA(NNN)=0.0
 20     CONTINUE
	XA(3)=3.0E-8
	XA(7)=1.0E-7
	IDAYS=DAYS*24/DHR
	NDAY=1
	HR=0.5
	DDT=DHR*720.
CCCCC  THIS NEXT LOOP WILL INITIALIZE THE X'S WITH "CLEAN" AIR
	DO 25 NNN=1,NXS
	X(NNN)=XA(NNN)
	A(NNN)=0.0
	B(NNN)=0.0
 25     CONTINUE
	CALL SMOG(R,THICKJ,THINJ,H20,TADV,RTADV,QHC,QNO,QCO,XA,X,FLUX,A,B)
	DO 30 II=1,NXS
CCCCC  DOUBLE THE A'S AND B'S SINCE THE SMOG SUBROUTINE GETS AN AVE A + B
	A(II)=A(II)*2.0
	B(II)=B(II)*2.0
	XS(II)=XA(II)
 30     CONTINUE
CCCCC  EVERYTHING UP TO THIS POINT IS INITIALIZATION
CCCCC  NOW BEGINS THE MAIN INTEGRATION LOOPS
	DO 60 I=1,IDAYS
	DO 46 IJJ=1,5
		DO 40 K=1,NXS
		X(K)=(XS(K)+A(K)*DDT)/(1.0+B(K)*DDT)
 40             CONTINUE
CCCCC  CALCULATE NEW A'S AND B'S WITH THE SMOG SUBROUTINE
	CALL SMOG(R,THICKJ,THINJ,H2O,TADV,RTADV,QHC,QNO,QCO,XA,X,FLUX,A,B)
		DO 45 L=1,NXS
		X(L)=(XS(L)+A(L)*DDT)/(1.0+B(L)*DDT)
		XS(L)=X(L)
 45             CONTINUE
 46             CONTINUE
CCCCC  NOW IS THE TIME TO STORE SOME ANSWERS
	WRITE(88,901)NDAY,HR,FLUX,(X(INXS),INXS=1,NXS)
901     format(1x,I4,',',F5.1,',',1P,E10.3,50(',',E10.3))
	HR=HR+DHR
	IF(HR.GE.24.0)GO TO 50
	GO TO 55
 50     HR=0.0
	NDAY=NDAY+1
 55     CONTINUE
	FLUX=0.1
	IF(HR.GE.6.0)FLUX=1.0
	IF(HR.GE.9.0)FLUX=0.4
	IF(HR.GE.16.0)FLUX=1.0
	IF(HR.GE.19)FLUX=0.4
	IF(HR.GE.6.0)THINJ=(1+COS(CONS*(HR-12.0)))/2.0
	IF(HR.GE.18)THINJ=0.0
	THICKJ=THINJ*THINJ
CCC J(NO2) for example = thinj*R(10)= 0.01 s-1 at noon
 60     CONTINUE
	END
	SUBROUTINE SMOG(R,THICKJ,THINJ,HHO,TADV,RTADV,QHC,QNO,QCO,XA,
     1X,FLUX,A,B)
	DIMENSION R(50),X(50),XA(50),A(50),B(50)
	XM=2.5E19
CCC Units of rate constants are as follows:  
CCC Sec-1 for first order reactions
CCC Atm/sec for second order reactions
CCC To convert second order rate constants into cm3/molecule*sec
CCC Multipy by 4.06E-20
	R(1)=1.4E+9
	R(2)=3.25E+7
	R(3)=1.4E+7
	R(4)=6.5
	R(5)=8.0E-4
	R(6)=2.375E+6
	R(7)=2.85E+7
	R(8)=3.5E+6
	R(9)=2.025E+8
CCCC  J(NO2) value = R(10)*thinj
CCCC  which reaches a maximum of 0.01 sec-1 at noon
	R(10)=1.0E-2
	R(11)=8.0E+4
	R(12)=5.5E-6
	R(13)=7.25E+8
	R(14)=5.75E+9
	R(15)=7.5E+8
	R(16)=3.0E+5
	R(17)=1.1E+6
	R(18)=2.025E+4
	R(19)=5.25E+2
	R(20)=5.0E-2
	R(21)=1.0E-3
	R(22)=1.8E-6
	R(23)=4.75E+8
	R(24)=3.0E+8
	R(25)=6.25E+7
	R(26)=1.6E-6
	R(27)=1.775E+7
	R(28)=1.4E+5
	R(29)=8.75E+8
	R(30)=3.0E+6
	R(31)=1.475E+8
	R(32)=1.625E+9
	R(33)=1.75E+9
	R(34)=3.75E+6
	R(35)=1.0E-3
	R(36)=1.0E-4
	HHO=1.0E-2
	A(1)=(A(1)+R(10)*THINJ*X(17)+R(13)*X(2)+R(21)*THINJ*X(3)+R(20)*T
     1HINJ*X(14)+XA(1)/TADV)/2.0
	A(2)=(A(2)+R(12)*THICKJ*X(3)+XA(2)/TADV)/2.
	A(3)=(A(3)+R(11)*X(1)+XA(3)/TADV)/2.
	A(4)=(A(4)+QNO*FLUX+R(10)*THINJ*X(17)+R(31)*X(1)*X(17)+XA(4)/TAD
     1V)/2.
	A(5)=(A(5)+R(9)*X(4)*X(6)+2.0*R(14)*X(2)*HHO+R(18)*X(6)*X(3)+2.0
     1*R(26)*THICKJ*X(15)+R(29)*X(1)*X(6)+R(30)*X(1)*X(9)+R(36)*THICKJ*X
     2(19)+XA(5)/TADV)/2.
	A(6)=(A(6)+R(4)*X(11)+R(6)*X(13)+R(7)*X(8)+R(17)*X(3)*X(5)+R(27)
     1*X(15)*X(5)+XA(6)/TADV)/2.
	A(7)=(A(7)+QCO*FLUX+R(7)*X(8)+R(22)*THICKJ*X(9)+XA(7)/TADV)/2.
	A(8)=(A(8)+R(5)*THINJ*X(9)+R(24)*X(5)*X(9)+R(30)*X(1)*X(9)+XA(8)
     1/TADV)/2.
	A(9)=(A(9)+R(4)*X(11)+XA(9)/TADV)/2.
	A(10)=(A(10)+R(1)*X(16)*X(5)+XA(10)/TADV)/2.
	A(11)=(A(11)+R(3)*X(12)*X(4)+XA(11)/TADV)/2.
	A(12)=(A(12)+R(2)*X(10)+XA(12)/TADV)/2.
	A(13)=(A(13)+R(5)*THINJ*X(9)+R(8)*X(7)*X(5)+XA(13)/TADV)/2.
	A(14)=(A(14)+R(19)*X(3)*X(17)+R(34)*X(5)*X(19)+R(35)*THINJ*X(18)
     1+XA(14)/TADV)/2.
	A(15)=(A(15)+R(25)*X(6)*X(6)+XA(15)/TADV)/2.
	A(16)=(A(16)+QHC*FLUX+XA(16)/TADV)/2.
	A(17)=(A(17)+R(3)*X(12)*X(4)+R(9)*X(4)*X(6)+R(16)*X(4)*X(3)+2.0*
     1R(23)*X(4)*X(14)+r(35)*THINJ*X(18)+QNO*FLUX*0.1+R(20)*X(14)*THINJ+
     2R(36)*THICKJ*X(19)+XA(17)/TADV)/2.
	A(18)=(A(18)+R(33)*X(17)*X(14)+XA(18)/TADV)/2.
	A(19)=(A(19)+R(32)*X(17)*X(5)+XA(19)/TADV)/2.
	B(1)=(B(1)+R(11)+R(28)*X(3)+R(29)*X(6)+R(30)*X(9)+R(31)*X(17)+RT
     1ADV)/2.
	B(2)=(B(2)+R(13)+R(14)*HHO+RTADV)/2.
	B(3)=(B(3)+R(12)*THICKJ+R(16)*X(4)+R(19)*X(17)+R(17)*X(5)+R(18)
     1*X(6)+R(21)*THINJ+R(28)*X(1)+RTADV)/2.
	B(4)=(B(4)+R(3)*X(12)+R(9)*X(6)+R(16)*X(3)+R(23)*X(14)+RTADV)/2.
	B(5)=(B(5)+R(8)*X(7)+R(1)*X(16)+R(15)*X(6)+R(17)*X(3)+R(24)*X(9
     1)+R(27)*X(15)+R(32)*X(17)+R(34)*X(19)+RTADV)/2.
	B(6)=(B(6)+R(9)*X(4)+R(15)*X(5)+R(18)*X(3)+2.0*R(25)*X(6)+R(29)*
     1X(1)+RTADV)/2.
	B(7)=(B(7)+R(8)*X(5)+RTADV)/2.
	B(8)=(B(8)+R(7)+RTADV)/2.
	B(9)=(B(9)+R(5)*THINJ+R(22)*THICKJ+R(24)*X(5)+R(30)*X(1)+RTADV)/
     12.
	B(10)=(B(10)+R(2)+RTADV)/2.
	B(11)=(B(11)+R(4)+RTADV)/2.
	B(12)=(B(12)+R(3)*X(4)+RTADV)/2.
	B(13)=(B(13)+R(6)+RTADV)/2.
	B(14)=(B(14)+R(23)*X(4)+R(33)*X(17)+R(20)*THINJ+RTADV)/2.
	B(15)=(B(15)+R(26)*THICKJ+R(27)*X(5)+RTADV)/2.
	B(16)=(B(16)+R(1)*X(5)+RTADV)/2.
	B(17)=(B(17)+R(10)*THINJ+R(19)*X(3)+R(31)*X(1)+R(33)*X(14)+R(32)
     1*X(5)+RTADV)/2.
	B(18)=(B(18)+R(35)*THINJ+RTADV)/2.
	B(19)=(B(19)+R(36)*THICKJ+R(34)*X(5)+RTADV)/2.
	RETURN
	END




			
	
