C***********************************************************************
C                                                                      *
C                                                                      *
C   PROGRAM TO GENERATE AEROSOL HEIGHT DISTRIBUTIO IN GAUSSIAN FORM    *
C                                                                      *
C   INPUT DATA: PEAK LOCATION AND WIDTH ARE INPUT    *
C                                                                      *
C	HISTORY: Modified 02/05/2001 to substitute SO2 with O2-O2  
C		Modified 04/20/2001 to make O2-O2 profile be consistent
C		with TOMRAD     
C	Modified 10/03/2001 to exclude O2-O2 and include SO2 again   
C       Modified 10/10/2001 to have the SO2 constant mixing ratio
C       withing the boundary layer: 
c     	Modified 2004	to add additional digits for surf. pressure 1013.2
c       MOdified 08/12/2015 [nk] to replace Gaussian with box-car aerosol 
c         vertical profie using peak[km] as the center height +/-width[km]
C***********************************************************************
C
C
c     subroutine aerprof(peak,width,so2max,so2width,totso2,nlvl,to2o2)
      subroutine aerprof(peak,width,so2max,so2width,totso2,nlvl)

C	totso2 - column So2 [atm cm]; so2max=0.0 and so2width=700.0 gives PBL pressure [mb]
      
      DIMENSION HS(100),PS(100),OZONE(100),AENUMD(100),tmp(100)
      real so2(100),so2max,so2width,totso2,loschmidt,coef1,buf
      real o2o2(100),to2o2,pn,tn,xo2,sum1,const
      parameter (pn=1013.25,tn=273.15,xo2=0.209476)
      DATA AENUMD/100*0.0/,so2/100*0.0/,loschmidt/2.6868E19/

      OPEN (10, FILE = 'profile.dat', STATUS = 'OLD')
      OPEN (11, FILE = 'pressure.prf', STATUS = 'unknown')
      OPEN (12, FILE = 'ozone.prf', STATUS = 'unknown')
      OPEN (13, FILE = 'temperature.prf', STATUS = 'unknown')
      OPEN (14, FILE = 'aerosol.prf', STATUS = 'unknown')
      OPEN (15, FILE = 'gas2.prf', STATUS = 'unknown')
      OPEN (16, FILE = 'water.prf', STATUS = 'unknown')
      
      pi=4.0*atan(1.0)

C
      WRITE(6,100) PEAK,WIDTH
  100 FORMAT(' GAUSSIAN PEAK AT ',F5.0,' KM WITH WIDTH ',F5.1,' KM ')

  110 FORMAT(/,I10,///)
      IF(NLVL.LT.0.OR.NLVL.GT.100) GO TO 999

      DO 10 I=1,NLVL
       READ(10,140) HS(I),PS(I),OZONE(I),tmp(i)
C140   FORMAT(F10.1,1P,4E10.3)
140   format(f10.1,1p,e12.5,2e10.3)	! add additional digits for surf. pressure 1013.25
   10 CONTINUE
C
C*** EVALUATE GAUSSIAN VALUE AT GIVEN HEIGHTS
C
      coef1=totso2*loschmidt/1.0E5
      const=34.163
      do i=1,nlvl
c     temp1=0.5*((hs(i)-so2max)/so2width)**2
      if(so2max.eq.0.0 .and. ps(i).ge.so2width) then
c---   p=nkT, k=1.381E-23 [J/K] Boltzmann constant
      so2(i)=1.0E-4*ps(i)/1.381E-23/tmp(i)
      else

c--   volcanic SO2 profile
c      temp1=0.5*((hs(i)-so2max)/so2width)**2
C----------------------o2o2 profile
C--------- n [#/cm^3] = 10^-6 * p[mb]*10^2 / k * T
C       k = 1.381*10^-23 [J/K] Boltzmann constant
c      o2o2(i)=(xo2*1.0E-4*ps(i)/1.381E-23/tmp(i))**2
c     if(temp1.le.180.0) so2(i)=exp(-temp1)
      so2(i)=0.0
      endif
      enddo

      sum=0.0
      sum1=0.0
      do i=1,nlvl-1
      sum=sum+0.5*(so2(i)+so2(i+1))*(hs(i+1)-hs(i))
c     sum1=sum1+0.5*(o2o2(i)+o2o2(i+1))*(hs(i+1)-hs(i))
      enddo

c     to2o2=sum1/loschmidt/loschmidt*1.0E5        !in cm*atm

      DO 20 I=1,NLVL
      AENUMD(I)=0.0
c     temp=0.5*((HS(I)-PEAK)/WIDTH)**2            
c     IF(TEMP.LE.180.0) AENUMD(I)=EXP(-TEMP)     ! Gaussian profile
      IF(ABS(HS(I)-PEAK).LE.WIDTH) AENUMD(I)=1.0           
      buf=coef1*so2(i)/sum
      WRITE(11,130) HS(I),PS(I)
      WRITE(12,130) HS(I),OZONE(I)
      WRITE(13,130) HS(I),tmp(i)
      WRITE(14,130) HS(I),AENUMD(I)
      write(15,130) hs(i),buf		! SO2 profile
      write(16,130) hs(i),buf		! SO2 profile
c     write(16,130) hs(i),o2o2(i)	! o2o2 profile
   20 CONTINUE
   
  130 FORMAT(F10.1,1p,E12.5)
      close(10)
      close(11)
      close(12)
      close(13)
      close(14)
      close(15)
      close(16)
      
c      call system('cp gas2.prf water.prf')
      return
      
  999 CONTINUE
      write(*,*) 'AEROSOL PROFILE ERROR: 
     &NUMBER OF LEVELS > 100, PROGRAM TERMINATED'
      close(10)
      close(11)
      close(12)
      close(13)
      close(14)
      close(15)

      stop
      END