c=======================================================================

c file:                  tauaero.f
c
c
c external routines:     uaerden,phaers,pharw,aerbwi,aerswi,usraer,
c                        angstrom,tauaero,relhum,stdaer
c
c required routines:     reverse, locate
c
c internal routines:     aeroden,aervint,aerswi
c
c internal common:       extd,aerblk
c=======================================================================

      subroutine aeroden(zz,v,aerd) 3,2

c purpose:   find number density of boundary layer aerosols, aerd,
c            at a given altitude, zz, and for a specified visibility
c input:
c   zz       altitude   (km)
c   v        visibility for a horizontal surface path (km) 
c output:
c   aerd     aerosol density at altitude z

c the vertical distribution of the boundary layer aerosol density is 
c based on the 5s vertical profile models for 5 and 23 km visibility.
c above 5 km, the aden05 and aden23 models are the same 
c below 5 km, the models differ as follows;
c aden05     0.99 km scale height (94% of extinction occurs below 5 km)
c aden23     1.45 km scale heigth (80% of extinction occurs below 5 km)
c    

      parameter (mz=33)
      dimension alt(mz),aden05(mz),aden23(mz)
      dimension zbaer(*),dbaer(*)
      save alt,aden05,aden23,nz

      data nz/mz/

      data alt/
     &        0.0,        1.0,        2.0,        3.0,        4.0,  
     &        5.0,        6.0,        7.0,        8.0,        9.0,
     &       10.0,       11.0,       12.0,       13.0,       14.0, 
     &       15.0,       16.0,       17.0,       18.0,       19.0,
     &       20.0,       21.0,       22.0,       23.0,       24.0, 
     &       25.0,       30.0,       35.0,       40.0,       45.0,
     &       50.0,       70.0,      100.0/
      data aden05/
     &  1.378E+04,  5.030E+03,  1.844E+03,  6.731E+02,  2.453E+02,      
     &  8.987E+01,  6.337E+01,  5.890E+01,  6.069E+01,  5.818E+01,
     &  5.675E+01,  5.317E+01,  5.585E+01,  5.156E+01,  5.048E+01,
     &  4.744E+01,  4.511E+01,  4.458E+01,  4.314E+01,  3.634E+01,
     &  2.667E+01,  1.933E+01,  1.455E+01,  1.113E+01,  8.826E+00,
     &  7.429E+00,  2.238E+00,  5.890E-01,  1.550E-01,  4.082E-02,
     &  1.078E-02,  5.550E-05,  1.969E-08/
      data aden23/
     &  2.828E+03,  1.244E+03,  5.371E+02,  2.256E+02,  1.192E+02,
     &  8.987E+01,  6.337E+01,  5.890E+01,  6.069E+01,  5.818E+01,
     &  5.675E+01,  5.317E+01,  5.585E+01,  5.156E+01,  5.048E+01,
     &  4.744E+01,  4.511E+01,  4.458E+01,  4.314E+01,  3.634E+01,
     &  2.667E+01,  1.933E+01,  1.455E+01,  1.113E+01,  8.826E+00,
     &  7.429E+00,  2.238E+00,  5.890E-01,  1.550E-01,  4.082E-02,
     &  1.078E-02,  5.550E-05,  1.969E-08/
c
      z=max(0.,min(100.,zz))
      aerd=0.
      if(z.gt.alt(nz)) return
        
      call locate(alt,nz,z,k)
      kp=k+1
      f=(z-alt(k))/(alt(kp)-alt(k))
      
      if(min(aden05(k),aden05(kp),aden23(k),aden23(kp)).le.0.) then
        aer05=aden05(k)*(1.-f)+aden05(kp)*f
        aer23=aden23(k)*(1.-f)+aden23(kp)*f
      else
        aer05=aden05(k)*(aden05(kp)/aden05(k))**f
        aer23=aden23(k)*(aden23(kp)/aden23(k))**f
      endif

      wth=(1./v-1/5.)/(1./23.-1./5.)
      wth=max(0.,min(1.,wth))

      aerd=(1.-wth)*aer05+wth*aer23

c      write(*,*) 'aeroden k,kp,z,aer05(k),aer05(kp),f,aerd'
c      write(*,'(2i5,1p5e11.3)') k,kp,z,aden05(k),aden05(kp),f,aerd

      return

c.......................................................................

      entry uaerden(nzb,zbaer,dbaer)
c
c purpose 
c       set the boundary layer aerosol vertical distribution.
c       This routine replaces data in aer05 and aer23. 
c
c input:
c
c   zbaer   altitude grid for aerosol density interpolation
c
c   dbaer   aerosol density at zbaer altitudes. The aerosol density at the
c           computational grid altitudes are obtained by interpolation
c           on the (zbaer,dbaer) grid.  If no elements of zbaer are set,
c           then dbaer represents the aerosol density at each
c           computational layer (unset values of dbaer have a default
c           value of -1, these values are pushed up to zero).
c 
      do 10 i=1,nz
        if(zbaer(1).lt.0.) then 
          aden05(i)=max(dbaer(i),0.)
        elseif(alt(i).lt.zbaer(1).or.alt(i).gt.zbaer(nzb)) then
          aden05(i)=0.
        else
          call locate(zbaer,nzb,alt(i),j)
          if(zbaer(j+1) .le. zbaer(j)) then 
            write(*,*) 'error -- incorrect specification of zbaer'
          endif
          w=(alt(i)-zbaer(j))/(zbaer(j+1)-zbaer(j))
          aden05(i)=dbaer(j)*(dbaer(j+1)/dbaer(j))**w
        endif
        
        aden23(i)=aden05(i)

c        if(i.eq.1) write(*,'(6x,7a11)')
c     &   'alt(i)','zbaer(j)','zbaer(j+1)','dbaer(j)',
c     &       'dbaer(j+1)','w','aden05(i)'
c        write(*,'(2i3,1p7e11.3)')
c     &   i,j,alt(i),zbaer(j),zbaer(j+1),dbaer(j),dbaer(j+1),w,aden05(i)
 10   continue

      return

      end