  subroutine callcompute_aod_raqms(pres,temp,prsi,rhi,aerosoli,aod,aodpart,aodlvl,KEcrtm,kgkg_kgm2)
  use mpimod, only : mype
  use kinds, only : r_kind,i_kind
  use gridmod, only : nsig
  use constants, only : rd,grav,r1000 
  implicit none
  real(r_kind) :: den(nsig),rh(nsig),dz(nsig),pres(nsig),temp(nsig),delp(nsig),presi(nsig+1)
  real(r_kind) :: aod,ext_3d(nsig),mass_ext_3d(nsig,14),mole_frac_3d(nsig,14),aodpart(14),aodlvl(nsig),KEcrtm(nsig,14)
  real(r_kind) :: aerosol(nsig,14),mass(nsig),prsi(nsig+1),rhi(nsig),aerosoli(nsig,14),kgkg_kgm2(nsig)
  logical ldiag
  common/dodiag/ldiag
! 1-5 sulf,bc12,oc12
! 6-9 du1234
! 11-14 ss1234
  integer(i_kind) :: k,m,nchem
  do k=1,nsig
!    if(rhi(k)<=0.0)then
!      write(6,*)'rh zero',k,rhi(k),'pres',pres(k)
!      call flush(6)
!    endif
    rh(k)=max(0.,min(99.,100.*rhi(k) ))
    mass(k)=(prsi(k)-prsi(k+1))*r1000/grav
    den(k)=pres(k)*r1000/(rd*temp(k))
!   dp=den*grav*dz
    dz(k)=mass(k)/den(k)
!   convert from ug/kg to kg/kg
    do m=1,14
      aerosol(k,m)=max(0.0,aerosoli(k,m))*1.e-9
    end do
  end do
  nchem=14
  aodlvl=0.0
  aodpart=0.0
  aod=0.0
  KEcrtm=0.0
  call compute_aod_raqms(nsig,nchem,den,rh,dz,aerosol,aod,nsig,ext_3d,mass_ext_3d,mole_frac_3d, &
  aodpart,aodlvl,KEcrtm,kgkg_kgm2)
  KEcrtm(:,10)=0.0 ! since raqms does not have dust 5
  return
  end subroutine callcompute_aod_raqms
      subroutine compute_aod_raqms(n1,nchem,den,rh,dz, &
                    chemp,aod,nlay,ext_3d, &
#define SAVEAOD
#ifdef SAVEAOD
                    mass_ext_3d,mole_frac_3d,aodpart, &
                    aodlvl,KEcrtm,kgkg_kgm2)
#else
                    mass_ext_3d,mole_frac_3d,KEcrtm,kgkg_kgm2)
#endif

! tks
! note subscripts on some arrays were changed.
!!!!!! arrays have different indices than in RAQMS
! rv is mixing ratio
! theta is theta here; in original code it was virtual theta

! This is based on 'compute_aod_RAQMS_ck.pro'

!  use funcphys
   use mpimod, only : mpi_comm_world,mpi_real8,mpi_real4,mpi_integer4,mpi_rtype
!   use module_assimaod, only : communicator,mpireal,mpiint,mpichar
!   USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
!   USE MODULE_GFS_PHYSCONS, only :  EPS => con_eps, EPSM1 => con_epsm1
  use mpimod, only : mype
  use mpeu_mpif, only : MPI_CHARACTER
  use kinds, only : r_kind,i_kind
  use computeaodraqms

  implicit none
  integer(i_kind) iam
  integer(i_kind) n1,nchem,nlay,mpicom,mpiint,mpireal,mpichar
! n1 is number of layers in atmosphere
! nlay is number of layers (from bottom) to do this calculation
  logical fexist,first
  save first
  data first/.true./
  real(r_kind) aod
#ifdef SAVEAOD
  real(r_kind) aodpart(14),aodlvl(n1),aodlvlsum,kgkg_kgm2(n1)
#endif
  real(r_kind) aod_so4,aod_bc1,aod_bc2
  real(r_kind) aod_oc1,aod_oc2,aod_no3
  real(r_kind) aod_du1_hlut,aod_du2_hlut
  real(r_kind) aod_du3_hlut,aod_du4_hlut,aod_du5_hlut
  real(r_kind) aod_ss1_hlut,aod_ss2_hlut
  real(r_kind) aod_ss3_hlut,aod_ss4_hlut
  real(r_kind) ,dimension(n1) :: den,dz
  real(r_kind) ,dimension(n1) :: rh
  real(r_kind) ,dimension(n1,nchem) :: chemp 
  real(r_kind) Hterm,Gterm
  real(r_kind) KEcrtm(n1,14)
  real(r_kind) rmin,rms,rmax,sss,ttt,bbb
  integer(i_kind) irh1,irh2
  integer(i_kind) nac,iiac,iz
!  parameter(nac=6)
  parameter(nac=5) ! we don't have nh3

!  real(r_kind) cp,R,cpr, cappa
!  parameter(cp=1004.,R=287.,cpr=cp/R, cappa= R/cp)
!match values in UW model (CCM3 physics)
!  parameter(cp=1004.64,R=287.04,cpr=cp/R, cappa= R/cp)
  integer(i_kind) k
!  real(r_kind) es,es1,es3,es4,es5,qs

!tks
  real(r_kind) aero_min
!tks
  real(r_kind) aero_conc(n1,nac_hlut)
  integer(i_kind) iac,itau

  real(r_kind) molef_rh(n1,nac),vmxr_w(n1,nac)
  real(r_kind) mass_w(n1,nac),ext(n1,nac)
  real(r_kind) reff_molef(n1,nac),mole_frac(n1,nac)
  real(r_kind) den_mix_rh(n1,nac),mass_ext(n1,nac)
  real(r_kind) mole_fractemp
  real(r_kind) qext_molef(n1,nac)
!  real(r_kind) pi,dz_in
  real(r_kind) pi
  real(r_kind) q_column(nac)

!  parameter(pi=3.14159265,dz_in=400.)
  parameter(pi=3.14159265)

! Wavelength for which AOD is calculated
  integer(i_kind) n_atau
  parameter(n_atau=1)
  real(r_kind) a_wli
  parameter(a_wli=0.55)

! Output AOD array
  real(r_kind) a_taus(n_atau,nac)

  integer(i_kind) nc
  character*10 ct,ct_new
  character*200 fl_aod

  real(r_kind) vmax,vmin
  real(r_kind) deltax,deltay,centrallat,centrallon
  real(r_kind) ext_3d(n1,nac_hlut)
  real(r_kind) mass_ext_3d(n1,nac_hlut)
  
  real(r_kind) mole_frac_3d(n1,nac_hlut)

  real(r_kind) beta_aero_hlut(n1,nac_hlut),q_column_hlut(nac_hlut)
  real(r_kind) f_w,den_aerow,mass_aerow,mass_aerow_min

! Wavelength for which AOD is calculated
  integer(i_kind) n_atau_hlut,ierr
  parameter(n_atau_hlut=1)
  real(r_kind) a_wli_hlut
  parameter(a_wli_hlut=0.55)

! Output AOD array
  real(r_kind) a_taus_hlut(n_atau_hlut,nac_hlut)
  real(r_kind) q_rh_hlut(n1,nac_hlut),reff_rh_hlut(n1,nac_hlut)
  logical ldiag
  common/dodiag/ldiag


!this is for GFS physics
!#ifdef USEGFUNCPHYS
!      call gfuncphys
!#endif


!tks 5-24-07 can be Nan otherwise
mpicom=mpi_comm_world
mpiint=mpi_integer4
mpireal=mpi_rtype
mpichar=mpi_character
reff_rh_hlut=0.
q_rh_hlut=0

! aerosol physical/chemical properties
!20060223 den_aero(1)=1.7e+3   ! sulfate
iam=mype

if(mype.eq.0)then
  masterproc=.true.
else
  masterproc=.false.
endif

! Initialize mole_frac_3d as 1.
do iac=1,nac_hlut

do k=1,nlay
 mole_frac_3d(k,iac)=1.
enddo
enddo

! others

aod=0.
aod_so4=0.
aod_bc1=0.
aod_bc2=0.
aod_oc1=0.
aod_oc2=0.
aod_no3=0.
!aod_so4_hlut=0.
!aod_bc1_hlut=0.
!aod_bc2_hlut=0.
!aod_oc1_hlut=0.
!aod_oc2_hlut=0.
!aod_no3_hlut=0.
aod_du1_hlut=0.
aod_du2_hlut=0.
aod_du3_hlut=0.
aod_du4_hlut=0.
aod_ss1_hlut=0.
aod_ss2_hlut=0.
aod_ss3_hlut=0.
aod_ss4_hlut=0.


!rh brought in now



!============================================================
do itau=1,n_atau
do iac=1,nac

!RAQMS Global
 iiac=iac
 q_column(iac)=0.

! Convert mixing ratio to mass concentration

 do k=1,nlay
 chemp(k,iiac)=max(chemp(k,iiac),0.)
! aero_conc(k,iac)=chemp(k,iiac)*mw_aero(iac)/ &
!                  mw_air*den(k) ! in kg/m3
! now in kG/KG
 aero_conc(k,iac)=chemp(k,iiac)*den(k)
 enddo


! Relative humidity (rh) dependency

do k=1,nlay
  ! Find two closest relative humidities from rh_lut
  do irh=1,nrh
   if(rh(k).le.rh_lut(irh)) goto 18
  enddo
! safer to say irh is nrh else could be nrh+1 which is out of bounds ajl
  irh=nrh
18 continue

  if(irh.eq.1) irh1=1
  if(irh.ge.2) irh1=irh-1
  if(irh.eq.1) irh2=1
  if(irh.ge.2) irh2=irh

  if(iac.eq.2) then
   molef_rh(k,iac)=1.
   den_mix_rh(k,iac)=den_mix_bc(1)
  endif
  if(iac.eq.4) then
   molef_rh(k,iac)=1.
   den_mix_rh(k,iac)=den_mix_oc(1)
  endif

!  if((rh_lut(irh2)-rh_lut(irh1))<1.) then
  if((rh_lut(irh2)-rh_lut(irh1))<1) then
    if(iac.eq.1) then
      molef_rh(k,iac)=mole_frac_so4(irh1)
      den_mix_rh(k,iac)=den_mix_so4(irh1)

    endif
    if(iac.eq.3) then
     molef_rh(k,iac)=mole_frac_bc(irh1)
     den_mix_rh(k,iac)=den_mix_bc(irh1)
    endif
    if(iac.eq.5) then
     molef_rh(k,iac)=mole_frac_oc(irh1)
     den_mix_rh(k,iac)=den_mix_oc(irh1)
    endif
    if(iac.eq.6) then
      molef_rh(k,iac)=mole_frac_no3(irh1)
      den_mix_rh(k,iac)=den_mix_no3(irh1)
    endif

!  endif

!  if((rh_lut(irh2)-rh_lut(irh1))>=1.) then
  else
    if(iac.eq.1) then 
      molef_rh(k,iac)=mole_frac_so4(irh1)+ &
      (mole_frac_so4(irh2)-mole_frac_so4(irh1))/ &
         (rh_lut(irh2)-rh_lut(irh1))* &
         (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.3) then
      molef_rh(k,iac)=mole_frac_bc(irh1)+ &
      (mole_frac_bc(irh2)-mole_frac_bc(irh1))/ &
         (rh_lut(irh2)-rh_lut(irh1))* &
         (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.5) then
      molef_rh(k,iac)=mole_frac_oc(irh1)+ &
     (mole_frac_oc(irh2)-mole_frac_oc(irh1))/ &
         (rh_lut(irh2)-rh_lut(irh1))* &
         (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.6) then
      molef_rh(k,iac)=mole_frac_no3(irh1)+ &
      (mole_frac_no3(irh2)-mole_frac_no3(irh1))/ &
         (rh_lut(irh2)-rh_lut(irh1))* &
         (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.1) then
      den_mix_rh(k,iac)=den_mix_so4(irh1)+ &
            (den_mix_so4(irh2)-den_mix_so4(irh1))/ &
            (rh_lut(irh2)-rh_lut(irh1))* &
            (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.3) then
      den_mix_rh(k,iac)=den_mix_bc(irh1)+ &
            (den_mix_bc(irh2)-den_mix_bc(irh1))/ &
            (rh_lut(irh2)-rh_lut(irh1))* &
            (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.5) then
      den_mix_rh(k,iac)=den_mix_oc(irh1)+ &
            (den_mix_oc(irh2)-den_mix_oc(irh1))/ &
            (rh_lut(irh2)-rh_lut(irh1))* &
            (rh(k)-rh_lut(irh1))
    endif
    if(iac.eq.6) then
      den_mix_rh(k,iac)=den_mix_no3(irh1)+ &
            (den_mix_no3(irh2)-den_mix_no3(irh1))/ &
            (rh_lut(irh2)-rh_lut(irh1))* &
            (rh(k)-rh_lut(irh1))
    endif

  endif ! end if for rh

  den_mix_rh(k,iac)=min(den_mix_rh(k,iac),1800.)
  den_mix_rh(k,iac)=max(den_mix_rh(k,iac),1000.)


  vmxr_w(k,iac)=(1.-molef_rh(k,iac))/molef_rh(k,iac)* &
                chemp(k,iiac)

!tks   mass_w(k,iac)=vmxr_w(k,iac)*mw_w/mw_air*den(k,i,j)

  mass_w(k,iac)=vmxr_w(k,iac)*mw_w/mw_air*den(k)
  Hterm=mw_w/mw_air*(1.-molef_rh(k,iac))/molef_rh(k,iac)
  mole_frac(k,iac)=1.

    mole_frac(k,iac)=molef_rh(k,iac)


! Find two closest mole fractions from lut
  do im=1,nmole_frac
   if(mole_frac(k,iac).le.mole_frac_lut(im)) goto 19
  enddo
  im=nmole_frac
  
19 continue
  if(im.eq.1) im1=1
  if(im.ge.2) im1=im-1
  if(im.eq.1) im2=1
  if(im.ge.2) im2=im

  if(iac.eq.2) then
   qext_molef(k,iac)=qext_bc(nmole_frac)
   reff_molef(k,iac)=reff_bc(nmole_frac)
  endif
  if(iac.eq.4) then
   qext_molef(k,iac)=qext_oc(nmole_frac)
   reff_molef(k,iac)=reff_oc(nmole_frac)
  endif
  
  if((mole_frac_lut(im2)-mole_frac_lut(im1))<0.01) then
  if(iac.eq.1) then
   qext_molef(k,iac)=qext_so4(im1)
   reff_molef(k,iac)=reff_so4(im1)
  endif
  if(iac.eq.3) then
   qext_molef(k,iac)=qext_bc(im1)
   reff_molef(k,iac)=reff_bc(im1)
  endif
  if(iac.eq.5) then
   qext_molef(k,iac)=qext_oc(im1)
   reff_molef(k,iac)=reff_oc(im1)
  endif
  if(iac.eq.6) then
   qext_molef(k,iac)=qext_no3(im1)
   reff_molef(k,iac)=reff_no3(im1)
  endif
  endif


  if((mole_frac_lut(im2)-mole_frac_lut(im1))>=0.01) then

  if(iac.eq.1) then
   qext_molef(k,iac)=qext_so4(im1)+ &
            (qext_so4(im2)-qext_so4(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))

   reff_molef(k,iac)=reff_so4(im1)+ &
            (reff_so4(im2)-reff_so4(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
  endif

  if(iac.eq.3) then
   qext_molef(k,iac)=qext_bc(im1)+ &
            (qext_bc(im2)-qext_bc(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
   reff_molef(k,iac)=reff_bc(im1)+ &
            (reff_bc(im2)-reff_bc(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
  endif

  if(iac.eq.5) then
   qext_molef(k,iac)=qext_oc(im1)+ &
            (qext_oc(im2)-qext_oc(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
   reff_molef(k,iac)=reff_oc(im1)+ &
            (reff_oc(im2)-reff_oc(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
  endif

  if(iac.eq.6) then
   qext_molef(k,iac)=qext_no3(im1)+ &
            (qext_no3(im2)-qext_no3(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
   reff_molef(k,iac)=reff_no3(im1)+ &
            (reff_no3(im2)-reff_no3(im1))/ &
            (mole_frac_lut(im2)-mole_frac_lut(im1))* &
            (mole_frac(k,iac)-mole_frac_lut(im1))
  endif

  endif ! end of mole_frac


  mass_ext(k,iac)=3./4.*qext_molef(k,iac)/ &
                  (den_mix_rh(k,iac)*reff_molef(k,iac))
  KEcrtm(k,iac)=mass_ext(k,iac)*(1.+Hterm)

  ext(k,iac)=mass_ext(k,iac)*(aero_conc(k,iac)+mass_w(k,iac))
  q_column(iac)=q_column(iac)+ext(k,iac)*dz(k)
#ifdef SAVEAOD
  aodlvl(k)=aodlvl(k)+ext(k,iac)*dz(k)
#endif
  
 
  ext_3d(k,iac)=ext(k,iac)
  mass_ext_3d(k,iac)=mass_ext(k,iac)
  mole_frac_3d(k,iac)=mole_frac(k,iac)

 enddo   ! end loop over layers

 a_taus(itau,iac)=q_column(iac)


enddo  ! end loop over species
enddo  ! end loop over wavelength

do itau=1,n_atau
 do iac=1,nac
  aod=aod+a_taus(itau,iac)
  if(iac.eq.1) aod_so4=aod_so4+a_taus(itau,iac)
  if(iac.eq.2) aod_bc1=aod_bc1+a_taus(itau,iac)
  if(iac.eq.3) aod_bc2=aod_bc2+a_taus(itau,iac)
  if(iac.eq.4) aod_oc1=aod_oc1+a_taus(itau,iac)
  if(iac.eq.5) aod_oc2=aod_oc2+a_taus(itau,iac)
  if(iac.eq.6) aod_no3=aod_no3+a_taus(itau,iac)
 enddo
enddo









!-------------------------------------------------------------------------



do itau=1,n_atau_hlut  ! loop ever wavelengths

   do ilamda=1,nlamda_hlut
     if(abs(lamda_hlut(ilamda)-a_wli_hlut*1000.).lt.1.) goto 9
   enddo
   write(6,*)'lamda is not found'
   call flush(6)
!   call killit('no lamda')
! stop
9 continue
!print *,'ilamda=',ilamda


!do iac=1,nac_hlut
do iac=7,nac_hlut
!Chieko's

!RAQMS Global (index of chemp)
 if(iac<=10)then
   iiac=iac-1 ! since one missing in chemp
 else
   iiac=iac ! sind no dust 5
 endif

 q_column_hlut(iac)=0.

 do ispec=1,nspec_hlut
   if(index(spec_line(ispec),trim(spec_name(iac))).gt.0) goto 99
 enddo
 write(100+iam,*)'spec is not found',spec_line(ispec)
 call flush(100+iam)
 do ispec=1,nspec_hlut
   write(100+iam,*)'spec_name',spec_name(iac)
   call flush(100+iam)
 end do
! call killit('spec not found')
99 continue

!============================================================



! Convert mixing ratio to mass concentration
! do k=1,n1
 do k=1,nlay

! aero_conc(k,iac)=chemp(k,iiac)*mw_aero(iac)/ &
!                  mw_air*den(k) ! in kg/m3
! now in kg/kg and not ppv
 aero_conc(k,iac)=chemp(k,iiac)*den(k)
 enddo

!-------------------------------------------!
!  q(nlamda,nspec),  reff(nlamda,nspec)
!  ssa(nlamda,nspec),pf(nterm,nlamda,nspec)
!-------------------------------------------!


! Exitinction coefficient at lamda

! ispec:ispec+6 code was put in place of original
! ispec:ispec+7 from Chieko.  Chieko and Brad 
! reviewed this code and recommended the change 4-19-07

if(iac.le.6.or.iac.ge.11) & 
 q_tmp_hlut(1:7)=q_hlut(ilamda,ispec:ispec+6)
if(iac.ge.7.and.iac.le.10) &
 q_tmp_hlut(1)=q_hlut(ilamda,ispec)

! Effective radius
if(iac.le.6.or.iac.ge.11) &
 reff_tmp_hlut(1:7)=reff_hlut(ilamda,ispec:ispec+6)
if(iac.ge.7.and.iac.le.10) &
 reff_tmp_hlut(1)=reff_hlut(ilamda,ispec)


! Relative humidity (rh) dependency
!do k=1,n1
do k=1,nlay
  ! Find two closest relative humidities from rh_lut
  do irh=1,7
   if(rh(k).le.rh_hlut(irh)) goto 8
  enddo
  irh=7
8 continue
  if(irh.eq.1) irh1=1
  if(irh.ge.2) irh1=irh-1
  if(irh.eq.1) irh2=1
  if(irh.ge.2) irh2=irh
  if(irh.eq.1)then
      write(700+mype,*)'rh small ',rh(k),'irh1',irh1,irh2
  endif
  if((rh_hlut(irh2)-rh_hlut(irh1))>=1.) &
  q_rh_hlut(k,iac)=q_tmp_hlut(irh1)+ &
         (q_tmp_hlut(irh2)-q_tmp_hlut(irh1))/ &
         (rh_hlut(irh2)-rh_hlut(irh1))* &
         (rh(k)-rh_hlut(irh1))
  if(iac.eq.2.or.iac.eq.4) q_rh_hlut(k,iac)=q_tmp_hlut(1)

  if(iac.ge.7.and.iac.le.10) q_rh_hlut(k,iac)=q_tmp_hlut(1)



  if((rh_hlut(irh2)-rh_hlut(irh1))>=1.) then
  reff_rh_hlut(k,iac)=reff_tmp_hlut(irh1)+ &
            (reff_tmp_hlut(irh2)-reff_tmp_hlut(irh1))/ &
            (rh_hlut(irh2)-rh_hlut(irh1))* &
            (rh(k)-rh_hlut(irh1))
!
! rbp fix hydrophillic bc, oc and dust reff bug (no re-initialize)
!
!   may not want inside of this loop
    if(iac.eq.2.or.iac.eq.4) reff_rh_hlut(k,iac)=reff_tmp_hlut(1)

    if(iac.ge.7.and.iac.le.10) reff_rh_hlut(k,iac)=reff_tmp_hlut(1)
   endif

  reff_rh_hlut(k,iac)=max(reff_rh_hlut(k,iac),1.e-10)
  reff_rh_hlut(k,iac)=min(reff_rh_hlut(k,iac),10.e-6)




!
! rbp fix hydrophillic bc, oc and dust ssa bug (no re-initialize)
!



!tks    reff_rh_hlut can't be less that 1e-10 because of 
!of the max(), min() 20 lines above

  if(reff_rh_hlut(k,iac)>=1.e-15) &
   f_w=reff_tmp_hlut(1)**3/reff_rh_hlut(k,iac)**3
  if(reff_rh_hlut(k,iac)<1.e-15) &
   f_w=1.
  den_aerow=f_w*den_aero(iac)+(1.-f_w)*1.e+3

!  mass_aerow=aero_conc(k,iac)
!   +4./3.*pi*(reff_rh(k,iac)**3)*den_aerow

  mass_aerow=1./f_w*den_aerow/den_aero(iac)*aero_conc(k,iac)
  Gterm=1./f_w*den_aerow/den_aero(iac)

  ! beta: specific extinction (m2/kg)
!tks 4-24-07

!   aero_min=1.e-20
!   aero_min=1.e-30
   aero_min=0.0

!    aero_min=1.e-12

!tks 3-27-08 added after discussion with Brad, avoids possible division
! by zero

  mass_aerow_min=1./f_w*den_aerow/den_aero(iac)*aero_min
  if(mass_aerow.lt.mass_aerow_min)mass_aerow=mass_aerow_min
! tks 3-27-08


!tks  if(aero_conc(k,iac)>=1.e-12) &
!tks 3-27-08  if(aero_conc(k,iac)>=aero_min) &

  beta_aero_hlut(k,iac)=3.*q_rh_hlut(k,iac) &
   /(4.*den_aerow*reff_rh_hlut(k,iac))


!tks  if(aero_conc(k,iac) < 1.e-12) beta_aero_hlut(k,iac)=0.
!tks 3-27-08  if(aero_conc(k,iac) < aero_min) beta_aero_hlut(k,iac)=0.
  q_column_hlut(iac)=q_column_hlut(iac)+ &
    beta_aero_hlut(k,iac)*(aero_conc(k,iac)+mass_aerow)*dz(k)
  if(iac.ge.7.and.iac.le.14) then 
    ext_3d(k,iac)=beta_aero_hlut(k,iac)*(aero_conc(k,iac)+mass_aerow)
    mass_ext_3d(k,iac)=beta_aero_hlut(k,iac)
#ifdef SAVEAOD
    aodlvl(k)=aodlvl(k)+ &
    beta_aero_hlut(k,iac)*(aero_conc(k,iac)+mass_aerow)*dz(k)
    if(iac.le.10)then
      KEcrtm(k,iac-1)=beta_aero_hlut(k,iac)*(1.+Gterm)
      if(rh(k)<=.1.or.irh1.eq.irh2)then
         write(800+mype,*)'ddirh1',irh1,irh2,'rh',rh(k),'KEcrtm',KEcrtm(k,iac),'iac',iac,'q',q_rh_hlut(k,iac), &
        'den',den_aerow,'reff',reff_rh_hlut(k,iac)
      endif
    else
      KEcrtm(k,iac)=beta_aero_hlut(k,iac)*(1.+Gterm)
      if(rh(k)<=.1.or.irh1.eq.irh2)then
         write(800+mype,*)'ssirh1',irh1,irh2,'rh',rh(k),'KEcrtm',KEcrtm(k,iac),'iac',iac,'q',q_rh_hlut(k,iac), &
        'den',den_aerow,'reff',reff_rh_hlut(k,iac)
      endif
    endif
#endif
  endif

 enddo   ! end loop over layers

 a_taus_hlut(itau,iac)=q_column_hlut(iac)


enddo  ! end loop over species
enddo  ! end loop over wavelengths

do itau=1,n_atau_hlut
!  do iac=1,nac_hlut
  do iac=7,nac_hlut
!20880327  if(iac.ge.7.and.iac.le.10) &
  if(iac.ge.7.and.iac.le.14)then
    aod=aod+a_taus_hlut(itau,iac)
  endif
  if(iac.eq.7) aod_du1_hlut=aod_du1_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.8) aod_du2_hlut=aod_du2_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.9) aod_du3_hlut=aod_du3_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.10) aod_du4_hlut=aod_du4_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.11) aod_ss1_hlut=aod_ss1_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.12) aod_ss2_hlut=aod_ss2_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.13) aod_ss3_hlut=aod_ss3_hlut+a_taus_hlut(itau,iac)
  if(iac.eq.14) aod_ss4_hlut=aod_ss4_hlut+a_taus_hlut(itau,iac)
 enddo
enddo
#ifdef SAVEAOD
aodlvlsum=0.0
do k=1,nlay
  aodlvlsum=aodlvlsum+aodlvl(k)
end do
if(abs(aodlvlsum-aod)>.0001)then
   write(6,*)'aodlvlsum',aodlvlsum,'aod',aod
aodlvlsum=0.0
   do k=1,nlay
     aodlvlsum=aodlvlsum+aodlvl(k)
     write(6,*)'aodlvl',k,aodlvl(k),' int ',aodlvlsum
   end do
endif
#endif



!!!!!!!!!!!!!!!!!!!!!!!!1147  continue

!============================================================
!!!tks change
!!!fl_aod='/idea_data1/Aero_fx/Output/AOD/uwnms_aod_'// &
!!!trim(ct_new)//'.dat'

#ifdef SAVEAOD
 aodpart(1)=aod_so4
 aodpart(2)=aod_bc1
 aodpart(3)=aod_bc2
 aodpart(4)=aod_oc1
 aodpart(5)=aod_oc2
 aodpart(6)=aod_du1_hlut
 aodpart(7)=aod_du2_hlut
 aodpart(8)=aod_du3_hlut
 aodpart(9)=aod_du4_hlut
! for GSI
 aodpart(10)=0.0
 aodpart(11)=aod_ss1_hlut
 aodpart(12)=aod_ss2_hlut
 aodpart(13)=aod_ss3_hlut
 aodpart(14)=aod_ss4_hlut
! aodpart(14)=aod_no3
! call mxmn2(aodpart(10),'ss1',0,1,rmax,rmin,rms)
#endif

first=.false.

return
end subroutine compute_aod_raqms 
