c *********************************************************************

c   *********************************************************************
c   * Program to calculate pollution in the neighborhood due to toxic *
c   * --------------------------------------------------------------- *
c   *            fumes vent out from stack heights: polc.f              *
c   *            -----------------------------------------              *
c   *********************************************************************
c
      parameter(h=100.0,q=1000.0,u=10.0,a=0.105,c=0.06)
      real wind(120,120),pol(120,120),max(120,120),th(120,120)
      real pi,const,lam,e,beta,theta,xb,yb,constw,ew1,ew2,rad
c
      pi=4.0*atan(1.0)
      rad=pi/180.0
c
      write(*,*)'What is the value of beta?'
      read(*,*)beta
c
      xb=sin(beta*rad)
      yb=cos(beta*rad)
c
      constm=q/(pi*u*a*c)
      const=constm*a/sqrt(2.0*pi)
c
      do 10 i=1,120
      do 20 j=1,120
c
      if (i .ge. 59 .and. i .le. 61) then
      if (j .ge. 59 .and. j .le. 61) then
      pol(i,j) = 0.0
      max(i,j) = 0.0
      goto 20
      else
      goto 1111
      endif
      else
      goto 1111
      endif
c
1111 lam = sqrt((i*100.0 - 6000.0)**2 + (j*100.0 - 6000.0)**2)
      fac = lam * c /sqrt((lam*c)**2 + (h*a)**2)
      e   = exp(-(h/(lam*c))**2/2.0)
c
      pol(i,j) = const * e * fac/(lam**2)
      max(i,j) = constm * e /(lam**2)
c
      ctheta = (xb*(j*100.0 - 6000.0) - yb*(i*100.0 - 6000.0))/lam
      theta = acos(ctheta)/rad
c
      if (theta .ge. 90.0) then
      wind(i,j) = 0.0
      elseif (theta .eq. 0.0) then
      wind(i,j) = max(i,j)
      else
      x = lam*sin((90-theta)*rad)
      y = lam*cos((90-theta)*rad)
      constw = q/(pi*u*a*c*(x**2))
      ew1 = exp(-(y**2)/(2.0*((a*x)**2)))
      ew2 = exp(-(h**2)/(2.0*((c*x)**2)))
      wind(i,j) = constw*ew1*ew2
      endif
c

20   continue
10   continue
9999 format(6(E9.3,x),/)
991 format(//)
992 format(' The value of of j and i are: I = ', I4,' to ', I4,
     c ' and J = ', I4,' to ',I4,//)

c
      write(50,*) ' The Flow of wind case '
      write(50,991)
      write(50, *) ' The value of Beta is ', beta
      write(50,991)
      write(50,992) 1,6,1,6
      write(50,9999) ((wind(i*10,j*10),j=1,6),i=1,6)
      write(50,991)
      write(50,992) 1,6,7,12
      write(50,9999) ((wind(i*10,j*10),j=1,6),i=7,12)
      write(50,991)
      write(50,992) 7,12,1,6
      write(50,9999) ((wind(i*10,j*10),j=7,12),i=1,6)
      write(50,991)
      write(50,992) 7,12,7,12
      write(50,9999) ((wind(i*10,j*10),j=7,12),i=7,12)

      write(51,*) ' Steady State Concentration Distribution '
      write(51,991)
      write(51,992) 1,6,1,6
      write(51,9999) ((pol(i*10,j*10),j=1,6),i=1,6)
      write(51,991)
      write(51,992) 1,6,7,12
      write(51,9999) ((pol(i*10,j*10),j=1,6),i=7,12)
      write(51,991)
      write(51,992) 7,12,1,6
      write(51,9999) ((pol(i*10,j*10),j=7,12),i=1,6)
      write(51,991)
      write(51,992) 7,12,7,12
      write(51,9999) ((pol(i*10,j*10),j=7,12),i=7,12)

      write(52, *) ' Maximum Concentration Case '
      write(52,991)
      write(52,992) 1,6,1,6
      write(52,9999) ((max(i*10,j*10),j=1,6),i=1,6)
      write(52,991)
      write(52,992) 1,6,7,12
      write(52,9999) ((max(i*10,j*10),j=1,6),i=7,12)
      write(52,991)
      write(52,992) 7,12,1,6
      write(52,9999) ((max(i*10,j*10),j=7,12),i=1,6)
      write(52,991)
      write(52,992) 7,12,7,12
      write(52,9999) ((max(i*10,j*10),j=7,12),i=7,12)
c
      stop
      end