program gaussian
        implicit none
        real  stackht,stackdia,emisrate,gasexitv,gasexitt
	real temp,atmcond,windvel(20),f0,x0,effht
	real sigmay,sigmaz,x(20),pi,cmaxa,cmaxb,cmaxc,cmax
	real scratch1,scratch2,scratch3,fudge
	real m1,m2
        integer i,npoints,wnpts,big
	pi=4.*(atan(1.))
	read*,stackht
	read*,stackdia
	read*,emisrate
	read*,gasexitv
	read*,gasexitt
	read*,temp
	read*,atmcond
	read*,wnpts
	read(*,*) (windvel(i),i=1,wnpts)
	read*,npoints
	read(*,*) (x(i),i=1,npoints)
	gasexitt=gasexitt+273.15
	temp=temp+273.15
	if(atmcond.lt.6.0) then
	   scratch1=(9.806*.02)/temp
        else
	   scratch1=(9.806*.035)/temp
        endif
	scratch3=3.12*0.785*gasexitv*stackdia**2*(gasexitt-temp)/gasexitt
c
	do 40 big=1,wnpts
C  Calculate the effective height of the stack
	f0=3.12*0.785*gasexitv*(stackdia**2)*
     +        ((gasexitt-temp)/gasexitt)
	if(f0.gt.55.0) then
	   x0=34.0*exp(0.4*log(f0))
        else
	   x0=14.0*exp(0.625*log(f0))
        endif
	if(atmcond.lt.5.0) then
	    effht=stackht+((1.6*exp(log(f0)/3.0))*
     +      exp((2.0*log(3.5*x0))/3.0))/windvel(big)
	else
	   m1=2.4*exp(log(scratch3/(scratch1*windvel(big)))/3.0)
	   m2=5.0*exp(log(scratch3/(scratch1*windvel(big)))/3.0)
	   fudge=min(m1,m2)
	   effht=stackht+fudge
        endif
c
	scratch2=(1.6*exp(log(f0)/3.0))*exp((2.0*log(3.5*x0))/3.0)
C
C Calculate the lateral and vertical dispersions
c
c        open(unit=5,file="/var/www/htdocs/tmp/plume.dat")
	do 20 i=1,npoints
	if (atmcond.eq.1.0) then
	    sigmay=x(i)*0.22*1000.*(1.0/(sqrt(1.0+0.1*x(i))))
	    sigmaz=x(i)*0.20*1000.0
        elseif(atmcond.eq.2.0) then
	    sigmay=x(i)*0.16*1000.*(1.0/(sqrt(1.0+0.1*x(i))))
	    sigmaz=x(i)*0.12*1000.
        elseif(atmcond.eq.3.0) then
            sigmay=x(i)*0.11*1000.*(1.0/(sqrt(1.0+0.1*x(i))))
	    sigmaz=x(i)*0.08*1000.*(1.0/(sqrt(1.0+0.2*x(i))))
        elseif(atmcond.eq.4.0) then
	    sigmay=x(i)*0.08*1000.*(1.0/(sqrt(1.0+0.1*x(i))))
	    sigmaz=x(i)*0.06*1000.*(1.0/(sqrt(1.0+1.5*x(i))))
        elseif(atmcond.eq.5.0) then
	    sigmay=x(i)*0.06*1000.*(1.0/(sqrt(1.0+0.1*x(i))))
	    sigmaz=x(i)*0.03*1000.*(1.0/(1.0+0.3*x(i)))
        else 
	    sigmay=x(i)*0.04*1000.*(1.0/(sqrt(1.0+0.1*x(i))))
	    sigmaz=x(i)*0.016*1000.*(1.0/(1.0+0.3*x(i)))
        endif
c
c  Calculate the concentration Cmax
C
        if((sigmay.lt.0.01).or.(sigmaz.lt.0.01)) then
	     cmax=0.0
        else
	     cmaxa=1.e6*emisrate/(2.0*pi*sigmay*sigmaz*windvel(big))
	     cmaxb=exp(-0.5*(effht/sigmaz)**2.)
	     cmaxc=exp(-0.5*(effht/sigmaz)**2.)
	     cmax=cmaxa*(cmaxb+cmaxc)
	endif
	print*,x(i),windvel(big),cmax
c	write(5,30) x(i),cmax
  20    continue
  40    continue
c  30    format(f7.2,2x,f7.2)
c	close(5)
        stop
	end