subroutine alumaerc

c Subroutine alumaer ------- David B. Considine, 11/11/96

c This subroutine calculates the surface area density for 
c Al2O3 particles produced by solid rocket motor emissions.

      save
      include 'com2d.h'

	COMMON/ALOX/ALO25(9),ALO35(9),RKALO25(Z$),RKALO35(Z$),
     *  AERAL(L$,Z$)

        COMMON/ALOXA/AERAL1(L$,Z$),AERAL2(L$,Z$),AERAL3(L$,Z$)

      real vfallal(l$,z$),
     c     vfallalmax,dtsedal,
     c     alfluxcorr,delgrid,ral2o3,rhoal2o3,ndensal,tempal,
     c     prodal,lossal

      integer idtsedal

c initialization:

      data init /0/

      if(init.eq.0)then

         init=1
         alfluxcorr=20.0     ! allows to adjust flux depending on dist
         delgrid=2.e5     ! vertical grid resolution, centimeters 
         ral2o3=0.0125     ! mean radius of al2o3 aerosol parts, microns
         rhoal2o3=3.98  ! density of al2o3 parts, g/cm^3
           
      endif

c Sedimentation calculation.  The code assumes that the
c Al2O3 particles follow an exponential distribution, with particle
c size distribution n(r)=(N/ro)exp(-r/ro).  The sedimentation
c flux depends on the size distribution. According to my notes
c (hetchem research, 11/1/96), The sedimenting flux for an
c exponential size distribution is: J=20*n*v(ro), where n is
c the number of condensed phase molecules per unit volume of
c air, and v(ro) is the fall velocity of a particle of size
c ro, the mean radius of the particle size distribution.
c This results in the setting of the variable fluxcorr to 20,
c and is equivalent to increasing the fall velocity by that
c amount.
      
c first calculate the fall velocity field and maximum fall velocity

      vfallalmax=0.
      vfallalmaxlat=0.
      vfallalmaxalt=0.

      do 10 ij=1,l$
      do 10 ik=1,z$

         ndensal=m(ij,ik)
         tempal=temp(ij,ik)

         vfallal(ij,ik)=
     c       alfluxcorr*fallvel(ral2o3,rhoal2o3,ndensal,tempal,0)

         if(vfallal(ij,ik).gt.vfallalmax)vfallalmax=vfallal(ij,ik)

 10   continue

c calculate time step for sedimentation.  The time step is calculated
c so that no gridbox can fall through the gridbox below it in one
c sedimentation time step.

      dtsedal=delgrid/vfallalmax

      if(dt/dtsedal.le.1.)then

         idtsedal=1
         dtsedal=dt             ! no need to change timestep

      else

         idtsedal=int(dt/dtsedal)+1 ! number of iterations necessary
         if(idtsedal.gt.10)idtsedal=10 ! limit number of iterations
         dtsedal=dt/(idtsedal*1.)
         vfallalmax=delgrid/dtsedal 

      endif

         
c Do sedimentation. Loop idtsedal times, use time step of dtsedal:


      do 20 ii=1,idtsedal
      do 20 ij=1,l$
      do 20 ik=1,z$-1
         
         if(vfallal(ij,ik+1).lt.vfallalmax)then
            prodal=dtsedal/delgrid*vfallal(ij,ik+1)*cn(71,ij,ik+1)
         else
            prodal=dtsedal/delgrid*vfallalmax*cn(71,ij,ik+1)
         endif

         if(vfallal(ij,ik).lt.vfallalmax)then
            lossal=dtsedal/delgrid*vfallal(ij,ik)
         else
            lossal=dtsedal/delgrid*vfallalmax
         endif

         cn(71,ij,ik)=(cn(71,ij,ik)+prodal)/(1.+lossal)
         
 20   continue

c Surface area density calculation.  This code assumes that the
c Al2O3 particles follow an exponential distribution, with particle
c size distribution n(r)=(N/ro)exp(-r/ro).  From Beiting,
c "Characteristics of Alumina Particles From Solid Rocket Motor
c Exhaust in the Stratosphere," Aerospace Report No. TR-95(5231)-8,
c ro=0.00933 microns and alpha (the mass density of alumina particles)
c is 4. grams/cm3.  The conversion from condensed molecule volume
c number density to surface area density is accoring to notes of
c 11/12/96.

c number of Al2O3 molecules per cm^3 of condensed phase

      alphal2o3=rhoal2o3/(101.9612*1.66e-24)

      do 30 ij=1,l$
      do 30 ik=1,z$

         aeral3(ij,ik)=cn(71,ij,ik)/(alphal2o3*ral2o3*1.e-4)

c The 1e-4 converts from microns to cm

 30   continue

      return
      end