subroutine alumaera

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.5618     ! mean radius of al2o3 aerosol parts, microns
         rhoal2o3=2.7   ! 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(69,ij,ik+1)
         else
            prodal=dtsedal/delgrid*vfallalmax*cn(69,ij,ik+1)
         endif

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

         cn(69,ij,ik)=(cn(69,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.6 microns and alpha (the mass density of alumina particles)
c is 2.6 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$

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

c The 1e-4 converts from microns to cm

 30   continue

      return
      end

      function fallvel(rad,density,ndens,tem,iwrite)

c Function fallvel ---------- David B. Considine, 11/11/96

c Code to calculate the fall velocity for a particle of radius r,
c which is passed to the function.  The fall velocity is calculated
c according to Kasten [1968]. (Kasten, F., Falling Speed of Aerosol
c Particles, J. Appl. Met., Oct, 1968

c declare variables

      real radius,rho,ndens,tem,sigsq,mfp,bet,s,dynvis,a,b,cc,
     c     grav,fallvel,rad,density

      integer iwrite
c code assumes radius is in microns and rho is in grams/cm^3.
c convert radius to meters and rho to kg/m^3:

      
      radius=rad*1.e-6
      rho=density*1.e3 

c Calculate mean free path. According to the U.S. Standard Atmosphere,
c 1976, the mean free path is given by the relationship: mfp =
c sqrt(2)/(2*pi*sig**2*N), where sig is the effective collision diameter
c (3.65e-10 meters) and N is the number density (units of #/meter**3)
c mfp units are 1/meters

      sigsq=1.3323e-19
      mfp=.22508/(sigsq*ndens*1.e6)

c Calculate dynamic viscosity. The formula is:
c dynvis=bet*T**(3/2)/(T+S), where bet is ..., s is Sutherland's
c constant (110.4 Kelvin).  The units are Newtons/sec/meter^2.

      bet=1.458e-6
      s=110.4
      dynvis=bet*tem**1.5/(tem+s)

c Calculate fall velocity.  The factor of 100 results in cm/s
c output.

      a=1.249
      b=0.42
      cc=0.87
      grav=9.8

      fallvel=0.2222*rho*radius*radius*grav/dynvis*100.*
     c                    (1.+ mfp/radius*(a+b*exp(-cc*radius/mfp)))


      if(iwrite.eq.1)write(6,*)'in fallvel',rad,radius,rho,density,
     c ndens,tem,dynvis,mfp,a,b,cc,fallvel,s,bet,sigsq

      return
      end