module computeaodraqms
  use kinds, only : i_kind,r_kind
  implicit none
#define NORAQMSMOD
#ifndef RAQMSMOD
  logical masterproc
#endif
  integer(i_kind) nac_hlut,irh
  parameter (nac_hlut=14)
  real(r_kind) mw_aero(nac_hlut),mw_air,mw_w
  parameter(mw_w=18.)
  real(r_kind) den_aero(nac_hlut)
! Variables for Harvard Lookup table

  character(len=50)  infl_qext
  character(len=10)  spec_name(nac_hlut)
  integer(i_kind) ispec,ilamda,iterm
  integer(i_kind) nspec_hlut,nlamda_hlut,nterm_hlut
  parameter (nspec_hlut=42,nlamda_hlut=7,nterm_hlut=8)

  real(r_kind) rh_hlut(7)

  real(r_kind) q_hlut(nlamda_hlut,nspec_hlut)
  real(r_kind) reff_hlut(nlamda_hlut,nspec_hlut)
  real(r_kind) lamda_hlut(nlamda_hlut)

  real(r_kind) q_tmp_hlut(7),reff_tmp_hlut(7)
  character*100 infl,indir,date*10
  integer(i_kind) nrh
  parameter(nrh=99)
  integer(i_kind) rh_lut(nrh)
  real(r_kind) mole_frac_so4(nrh),den_mix_so4(nrh)
  real(r_kind) mole_frac_bc(nrh),den_mix_bc(nrh)
  real(r_kind) mole_frac_oc(nrh),den_mix_oc(nrh)
  real(r_kind) mole_frac_no3(nrh),den_mix_no3(nrh)
  integer(i_kind) im,im1,im2
  integer(i_kind) nmole_frac2,nmole_frac
  parameter(nmole_frac=100)
  real(r_kind) moledensend(2,4,nrh)
  real(r_kind) qrsend(2,4,nmole_frac)
  real(r_kind) qext_tmp,reff_tmp
  real(r_kind) mole_frac_lut(nmole_frac)
  real(r_kind) qext_so4(nmole_frac)
  real(r_kind) qext_bc(nmole_frac)
  real(r_kind) qext_oc(nmole_frac)
  real(r_kind) qext_no3(nmole_frac)
  real(r_kind) reff_so4(nmole_frac)
  real(r_kind) reff_bc(nmole_frac)
  real(r_kind) reff_oc(nmole_frac)
  real(r_kind) reff_no3(nmole_frac)
  integer(i_kind) nest,iest
  parameter(nest=150)
  real(r_kind) estbar(nest)

  data rh_hlut/0.,50.,70.,80.,90.,95.,99./
  character(len=30)  spec_line(nspec_hlut)
  integer(i_kind) :: nrhaod
!  parameter (nrhaod=3000)
  parameter (nrhaod=10000)
!  real(r_kind) :: KEaod(0:nrhaod,14),reff(0:nrhaod,14)
  real(r_kind),allocatable :: KEaod(:,:),reff(:,:)
  real(r_kind),allocatable :: KEcrtmrh(:,:)
  real(r_kind) :: divrhraq,rrmaxrhraq
  parameter (divrhraq=1./float(nrhaod),rrmaxrhraq=float(nrhaod))
end module computeaodraqms
  subroutine closecomputeaodraqms
  use computeaodraqms
  if(allocated(KEaod))then
    deallocate (KEaod,reff)
  endif
  if(allocated(KEcrtmrh))then
    deallocate (KEcrtmrh)
  endif
  return
  end subroutine closecomputeaodraqms
  subroutine initcomputeaodraqms
#ifdef RAQMSMODEL
  use kinds, only : i_kind,r_kind
  use pmgrid, only : iam,masterproc
  use mpishorthand
  use computeaodraqms
  integer(i_kind) nchem
  integer mype

#else
  use mpimod, only : mype
  use kinds, only : i_kind,r_kind
  use mpimod, only : mpi_comm_world,mpi_real8,mpi_real4,mpi_integer4,mpi_rtype
  use mpeu_mpif, only : MPI_CHARACTER
  use computeaodraqms
  integer(i_kind) nchem,mpicom,mpiint,mpireal, mpichar
#endif
  integer(i_kind) rh_tmp,ierr,k,m
  real(r_kind) ratio_tmp,mole_frac_tmp,den_mix_tmp
  real(r_kind) estbar_tmp
  character(len=30)  line
  real(r_kind) tmp_hlut(12),rh(0:nrhaod)
  if(.not.allocated(KEaod))then
    allocate(KEaod(0:nrhaod,14),reff(0:nrhaod,14))
  endif
#ifndef RAQMSMODEL
  mpicom=mpi_comm_world
  mpiint=mpi_integer4
  mpireal=mpi_rtype
  mpichar=mpi_character
  if(mype.eq.0)then
    masterproc=.true.
  else
    masterproc=.false.
  endif
#else
  mype=iam
#endif
  den_aero(1)=1.769e+3   ! sulfate
  den_aero(2)=1.0e+3   ! bc
  den_aero(3)=1.0e+3   ! bc
  den_aero(4)=1.4e+3   ! oc
  den_aero(5)=1.4e+3   ! oc
  den_aero(6)=1.725e+3   ! no3
  den_aero(7)=2500.
  den_aero(8)=2650.
  den_aero(9)=2650.
  den_aero(10)=2650.
  den_aero(11)=2200.
  den_aero(12)=2200.
  den_aero(13)=2200.
  den_aero(14)=2200.

!mw_aero(1)=96.   ! sulfate
  mw_aero(1)=132.   ! sulfate
  mw_aero(2)=12.   ! bc
  mw_aero(3)=12.   ! bc
  mw_aero(4)=16.8   ! oc
  mw_aero(5)=16.8   ! oc
  mw_aero(6)=79.    ! no3
  mw_aero(7)=28.97    !dust
  mw_aero(8)=28.97    !dust
  mw_aero(9)=28.97    !dust
  mw_aero(10)=28.97   !dust
!Chieko said to use 29 for dust, now consistent with calling program
!mw_aero(7)=60.1
!mw_aero(8)=60.1
!mw_aero(9)=60.1
!mw_aero(10)=60.1
!J said to use 28.97 (air) for seasalt
  mw_aero(11)=28.97   !seasalt
  mw_aero(12)=28.97   !seasalt
  mw_aero(13)=28.97   !seasalt
  mw_aero(14)=28.97   !seasalt

  spec_name(1)='sulfate'
  spec_name(2)='Black C'
  spec_name(3)='Black C'
  spec_name(4)='Organic C'
  spec_name(5)='Organic C'
  spec_name(6)='sulfate' ! use sulfate extinction efficiency for now
  spec_name(7)='Mdust 0.8'
  spec_name(8)='Mdust 1.5'
  spec_name(9)='Mdust 2.5'
  spec_name(10)='Mdust 4.0'
  spec_name(11)='SSa00 Sea Salt (accum)'
!spec_name(12)='SSa00 Sea Salt (accum)'
!changed 3-19-08
  spec_name(12)='SSc00 Sea Salt (coarse)'
  spec_name(13)='SSc00 Sea Salt (coarse)'
  spec_name(14)='SSc00 Sea Salt (coarse)'
!tks 20080328mw_air=28.9
  mw_air=28.97

  if(masterproc)then

!tks fix for now  8-25-11
  indir='/home/lenzen/raqms_data/INTEX/AOD_lut/'
!tks

! Read in values from LUP

    infl='hygroscopic_growth_factors_nh4so4_0.1um_273K.txt'

    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')


    do irh=1,nrh
      read (1,"(i12,f13.5,f13.7,f13.2)",advance="yes") &
      rh_tmp,ratio_tmp,mole_frac_tmp,den_mix_tmp
      rh_lut(irh)=rh_tmp
      mole_frac_so4(irh)=mole_frac_tmp
     den_mix_so4(irh)=den_mix_tmp
    enddo
    close(1)

    infl='hygroscopic_growth_factors_bc_0.1um_273K.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do irh=1,nrh
      read (1,"(i12,f13.5,f13.7,f13.2)",advance="yes") &
      rh_tmp,ratio_tmp,mole_frac_tmp,den_mix_tmp
      mole_frac_bc(irh)=mole_frac_tmp
      den_mix_bc(irh)=den_mix_tmp
    enddo
    close(1)

    infl='hygroscopic_growth_factors_oc_0.1um_273K.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do irh=1,nrh
      read (1,"(i12,f13.5,f13.7,f13.2)",advance="yes") &
      rh_tmp,ratio_tmp,mole_frac_tmp,den_mix_tmp
      mole_frac_oc(irh)=mole_frac_tmp
      den_mix_oc(irh)=den_mix_tmp
    enddo
    close(1)

    infl='hygroscopic_growth_factors_nh4no3_0.1um_273K.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do irh=1,nrh
      read (1,"(i12,f13.5,f13.7,f13.2)",advance="yes") &
      rh_tmp,ratio_tmp,mole_frac_tmp,den_mix_tmp
      mole_frac_no3(irh)=mole_frac_tmp
      den_mix_no3(irh)=den_mix_tmp
    enddo
    close(1)

    open(1,file=trim(indir)//'estbar.dat',status='old',form='formatted')
    estbar_tmp=0.
    do iest=1,nest
      read(1,*) estbar_tmp
      estbar(iest)=estbar_tmp
    enddo
    close(1)

! Read in values from LUT
    infl='scat_effic_nh4so4_r0.075_s2.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do im=1,nmole_frac
      read (1,"(3f13.6)",advance="yes") &
      mole_frac_tmp,qext_tmp,reff_tmp
      mole_frac_lut(im)=mole_frac_tmp

      qext_so4(im)=qext_tmp
      reff_so4(im)=reff_tmp
!      write(500+mype,*)'mole_frac_lut',im,mole_frac_lut(im),'reffso4',reff_so4(im)
    enddo
    close(1)

!   Read in values from LUT
    infl='scat_effic_bc_r0.075_s2.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do im=1,nmole_frac
      read (1,"(3f13.6)",advance="yes") &
      mole_frac_tmp,qext_tmp,reff_tmp
      qext_bc(im)=qext_tmp
      reff_bc(im)=reff_tmp
!      write(500+mype,*)'mole_frac_lut',im,mole_frac_lut(im),'reffbc',reff_bc(im)
    enddo
    close(1)

!   Read in values from LUT
    infl='scat_effic_oc_r0.075_s2.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do im=1,nmole_frac
      read (1,"(3f13.6)",advance="yes") &
      mole_frac_tmp,qext_tmp,reff_tmp
      qext_oc(im)=qext_tmp
      reff_oc(im)=reff_tmp
!      write(500+mype,*)'mole_frac_lut',im,mole_frac_lut(im),'reffoc',reff_oc(im)
    enddo
    close(1)

!   Read in values from LUT
   
    infl='scat_effic_nh4no3_r0.075_s2.txt'
    open (1,file=trim(indir)//trim(infl),status='old', &
        form='formatted')
    do im=1,nmole_frac
      read (1,"(3f13.6)",advance="yes") &
     mole_frac_tmp,qext_tmp,reff_tmp
     qext_no3(im)=qext_tmp
     reff_no3(im)=reff_tmp
    enddo
    close(1)
  endif
! bcast parameters
  moledensend(1,1,:)=mole_frac_so4(:)
  moledensend(1,2,:)=mole_frac_bc(:)
  moledensend(1,3,:)=mole_frac_oc(:)
  moledensend(1,4,:)=mole_frac_no3(:)
  moledensend(2,1,:)=den_mix_so4(:)
  moledensend(2,2,:)=den_mix_bc(:)
  moledensend(2,3,:)=den_mix_oc(:)
  moledensend(2,4,:)=den_mix_no3(:)

  call mpi_bcast(moledensend,8*nrh,mpireal,0,mpicom,ierr)
  mpiint=mpi_integer4

  call mpi_bcast(rh_lut,nrh,mpiint,0,mpicom,ierr)
  mole_frac_so4(:)=moledensend(1,1,:)
  mole_frac_bc(:)=moledensend(1,2,:)
  mole_frac_oc(:)=moledensend(1,3,:)
  mole_frac_no3(:)=moledensend(1,4,:)
  den_mix_so4(:)=moledensend(2,1,:)
  den_mix_bc(:)=moledensend(2,2,:)
  den_mix_oc(:)=moledensend(2,3,:)
  den_mix_no3(:)=moledensend(2,4,:)
  call mpi_bcast(estbar,nest,mpireal,0,mpicom,ierr)
  call mpi_bcast(mole_frac_lut,nmole_frac,mpireal,0,mpicom,ierr)
  qrsend(1,1,:)=qext_so4(:)
  qrsend(1,2,:)=qext_bc(:)
  qrsend(1,3,:)=qext_oc(:)
  qrsend(1,4,:)=qext_no3(:)
  qrsend(2,1,:)=reff_so4(:)
  qrsend(2,2,:)=reff_bc(:)
  qrsend(2,3,:)=reff_oc(:)
  qrsend(2,4,:)=reff_no3(:)
  call mpi_bcast(qrsend,8*nmole_frac,mpireal,0,mpicom,ierr)
  qext_so4(:)=qrsend(1,1,:)
  qext_bc(:)=qrsend(1,2,:)
  qext_oc(:)=qrsend(1,3,:)
  qext_no3(:)=qrsend(1,4,:)
  reff_so4(:)=qrsend(2,1,:)
  reff_bc(:)=qrsend(2,2,:)
  reff_oc(:)=qrsend(2,3,:)
  reff_no3(:)=qrsend(2,4,:)

! um to m
  reff_so4=reff_so4*1.e-6
  reff_bc=reff_bc*1.e-6
  reff_oc=reff_oc*1.e-6
  reff_no3=reff_no3*1.e-6
! Compute AOD using Harvard Lookup table

! Read in values from LUP

!tks
  if(masterproc)then
    infl_qext='Qext_Reff_GEOS-CHEM.dat'

    open (1,file=trim(indir)//trim(infl_qext),status='old', &
        form='formatted')
    read (1,"(a20)",advance="yes") line
    do ispec=1,nspec_hlut
      read (1,"(a20)",advance="yes") line
      spec_line(ispec)=line
      do ilamda=1,nlamda_hlut
        read (1,*) tmp_hlut
        lamda_hlut(ilamda)=tmp_hlut(1)
        q_hlut(ilamda,ispec)=tmp_hlut(2)
        reff_hlut(ilamda,ispec)=tmp_hlut(3)*1.e-6
     enddo
    enddo
    close(1)
  endif
  call mpi_bcast(spec_line,30*nspec_hlut,mpichar,0,mpicom,ierr)
  call mpi_bcast(lamda_hlut,nlamda_hlut,mpireal,0,mpicom,ierr)
  call mpi_bcast(q_hlut,nlamda_hlut*nspec_hlut,mpireal,0,mpicom,ierr)
  call mpi_bcast(reff_hlut,nlamda_hlut*nspec_hlut,mpireal,0,mpicom,ierr)
  do k=0,nrhaod
    rh(k)=min(99.,float(k)/float(nrhaod)*100.)
  end do
  rh(0)=1.e-20
  call compute_aod_KE_raqms(nrhaod+1,14,rh)
  KEaod(:,10)=0.0 ! since raqms has no dust 5
!  go to 300
!  write(6,*)iam,mype,'did compute_aod_KE_raqms ',nrhaod
!  call flush(6)
!  go to 300
  if(mype.eq.0)then
    open(30,file='/home/lenzen/GSI/raqmske.txt.raqms',form='formatted')
    open(31,file='/home/lenzen/GSI/raqmsreff.txt.raqms',form='formatted')
    write(30,*)' mass extinction m2/g '
    write(30,'("   rh    sulf    bc1    bc2    oc1    oc2    du1    du2    du3    du4    ss1    ss2    ss3    ss4")')
    write(31,*)'effective radius microns '
    write(31,'("   rh    sulf    bc1    bc2    oc1    oc2    du1    du2    du3    du4    ss1    ss2    ss3    ss4")')
    do k=0,nrhaod,30
      write(30,'(f6.1,1x,2f7.3,f8.3,10f7.3)')rh(k),(KEaod(k,m)*.001,m=1,9),(KEaod(k,m)*.001,m=11,14)
!      write(31,'(f6.2,1x,13f8.1)')rh(k),(reff(k,m)*1.e8,m=1,9),(reff(k,m)*1.e8,m=11,14)
      write(31,'(f6.1,1x,13f7.3)')rh(k),(reff(k,m)*1.e6,m=1,9),(reff(k,m)*1.e6,m=11,14)
    end do
    close(30)
    close(31)
  endif
  go to 300
  if(mype.eq.0)then
    write(6,*)'KE  SO4 BC1   BC2   OC1  OC2 '
    do k=0,1000
      write(6,'(f8.2,2x,5f9.1)')rh(k),KEaod(k,1:5)
    end do
    write(6,*)'KE du 1 - 4'
    do k=0,1000
      write(6,'(f8.2,2x,4f9.1)')rh(k),KEaod(k,6:9)
    end do
    write(6,*)'KE SS 1 - 4'
    do k=0,1000
      write(6,'(f8.2,2x,4f9.1)')rh(k),KEaod(k,11:14)
    end do
    write(6,*)'reff'
    do k=0,1000
      write(6,'(f8.2,2x,5f9.1)')rh(k),reff(k,1:5)*1.e8
    end do
    do k=0,1000
      write(6,'(f8.2,2x,4f9.1)')rh(k),reff(k,6:9)*1.e8
    end do
    do k=0,1000
      write(6,'(f8.2,2x,4f9.1)')rh(k),reff(k,11:14)*1.e8
    end do
    
  endif
300 continue
    
  end subroutine initcomputeaodraqms
  subroutine compute_aod_KE_raqms(nrhi,nchem,rh)

! This is based on 'compute_aod_RAQMS_ck.pro'

!  use funcphys
#ifdef RAQMSMODEL
   use pmgrid, only : iam
   use mpishorthand
#else
   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
#endif
  use kinds, only : r_kind,i_kind
  use computeaodraqms

  implicit none
#ifndef RAQMSMODEL
  integer(i_kind) iam
  integer(i_kind) nchem,nrhi,mpicom,mpiint,mpireal,mpichar
#else
  integer(i_kind) nchem,nrhi,mype
#endif
! nrhi 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) ,dimension(nrhi) :: rh
  real(r_kind) Hterm,Gterm
!  real(r_kind) KEcrtm(nrhi,14),reff(nrhi,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(nrhi,nac_hlut)
  integer(i_kind) iac,itau

  real(r_kind) molef_rh(nrhi,nac),vmxr_w(nrhi,nac)
  real(r_kind) mass_w(nrhi,nac),ext(nrhi,nac)
  real(r_kind) reff_molef(nrhi,nac),mole_frac(nrhi,nac)
  real(r_kind) den_mix_rh(nrhi,nac),mass_ext(nrhi,nac)
  real(r_kind) mole_fractemp
  real(r_kind) qext_molef(nrhi,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(nrhi,nac_hlut)
  real(r_kind) mass_ext_3d(nrhi,nac_hlut)
  
  real(r_kind) mole_frac_3d(nrhi,nac_hlut)

  real(r_kind) beta_aero_hlut(nrhi,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(nrhi,nac_hlut),reff_rh_hlut(nrhi,nac_hlut)
  logical ldiag
  common/dodiag/ldiag


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


#ifndef RAQMSMODEL
!tks 5-24-07 can be Nan otherwise
mpicom=mpi_comm_world
mpiint=mpi_integer4
mpireal=mpi_rtype
mpichar=mpi_character
iam=mype
if(mype.eq.0)then
  masterproc=.true.
else
  masterproc=.false.
endif
#else
 mype=iam
#endif
reff_rh_hlut=0.
q_rh_hlut=0

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


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

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



!rh brought in now



!============================================================
if(mype.eq.0)write(6,*)'nrhi',nrhi,'nac',nac,'n_atau',n_atau,'nchem',nchem
do itau=1,n_atau
do iac=1,nac

!RAQMS Global
 iiac=iac

! Convert mixing ratio to mass concentration



! Relative humidity (rh) dependency

do k=1,nrhi
  ! Find two closest relative humidities from rh_lut
  if(iac.eq.2) then
     molef_rh(k,iac)=1.
     den_mix_rh(k,iac)=den_mix_bc(1)
  elseif(iac.eq.4) then
     molef_rh(k,iac)=1.
     den_mix_rh(k,iac)=den_mix_oc(1)
  elseif(rh(k).le.rh_lut(1))then
! just use 1 for iac 1,3,5,6
    irh1=1
    if(iac.eq.1) then
      molef_rh(k,iac)=mole_frac_so4(irh1)
      den_mix_rh(k,iac)=den_mix_so4(irh1)
    elseif(iac.eq.3) then
      molef_rh(k,iac)=mole_frac_bc(irh1)
      den_mix_rh(k,iac)=den_mix_bc(irh1)
    elseif(iac.eq.5) then
      molef_rh(k,iac)=mole_frac_oc(irh1)
      den_mix_rh(k,iac)=den_mix_oc(irh1)
    elseif(iac.eq.6) then
      molef_rh(k,iac)=mole_frac_no3(irh1)
      den_mix_rh(k,iac)=den_mix_no3(irh1)
    endif
  else
!   need to search for rh interval

    do irh=1,nrhi
      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)then
      write(6,*)'irh should not be one'
       call stop2(444)
!      call killit('irh')
    endif
    irh1=irh-1
    irh2=irh

    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))
    elseif(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))
    elseif(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))
    elseif(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))
    elseif(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))
    elseif(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))
    elseif(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.)


!  Hterm=mw_w/mw_air*(1.-molef_rh(k,iac))/molef_rh(k,iac)
! fix ajl 4/26/2016
  Hterm=mw_w/mw_aero(iac)*(1.-molef_rh(k,iac))/molef_rh(k,iac)

  mole_frac(k,iac)=molef_rh(k,iac)
!    if(iac.eq.1)then
!    write(6,*)'rh',rh(k),'k',k,'molef_rh',molef_rh(k,iac)
!     endif


! Find two closest mole fractions from lut
!  if(k<=2.and.iac.eq.1.and.iam.eq.0)then
!     write(6,*)'mole_frac_lut',mole_frac_lut
!  endif
  if(iac.eq.2) then
   qext_molef(k,iac)=qext_bc(nmole_frac)
   reff_molef(k,iac)=reff_bc(nmole_frac)
  elseif(iac.eq.4) then
   qext_molef(k,iac)=qext_oc(nmole_frac)
   reff_molef(k,iac)=reff_oc(nmole_frac)
  elseif(mole_frac(k,iac).le.mole_frac_lut(1))then
    im1=1
    if(iac.eq.1) then
      qext_molef(k,iac)=qext_so4(im1)
      reff_molef(k,iac)=reff_so4(im1)
!     if(iac.eq.1.and.rh(k)>90.)then
!      write(500+mype,*)k,rh(k),' im1=im2',im1,im2,' qext ',qext_molef(k,iac),' reff ',reff_molef(k,iac)
    ! endif
!     if(iac.eq.1)then
!       write(500+mype,*)'k',k,'im1',im1,reff_molef(k,iac)
!      endif
    elseif(iac.eq.3) then
      qext_molef(k,iac)=qext_bc(im1)
      reff_molef(k,iac)=reff_bc(im1)
    elseif(iac.eq.5) then
      qext_molef(k,iac)=qext_oc(im1)
      reff_molef(k,iac)=reff_oc(im1)
    elseif(iac.eq.6) then
      qext_molef(k,iac)=qext_no3(im1)
      reff_molef(k,iac)=reff_no3(im1)
     endif
  else
    do im=1,nmole_frac
!   if(k<=10)then
!     write(500+mype,*)'im',im,'k',k,'mole_frac',mole_frac(k,iac),'mole_frac im',mole_frac_lut(im)
!   endif
     if(mole_frac(k,iac).le.mole_frac_lut(im)) goto 19
    enddo
    im=nmole_frac
  
19  continue
!    if(iac.eq.1.and.iam.eq.0.and.k>=9200.and.k<=9300)then
!       write(6,*)k,'im',im
!    endif
    if(im.eq.1)then
       write(6,*)'error im one'
       call stop2(445)
       !call killit('im')
    endif
    im1=im-1
    im2=im

!   if(iac.eq.1)then
!     write(500+mype,*)'k',k,'im',im,'im1',im1,im2,'Molefrac',mole_frac_lut(im1:im2)
!   endif
!    if(iac.eq.1.and.iam.eq.0.and.k>=9200.and.k<=9300)then
!      write(6,*)k,'im1',im1,im2,'mole_frac_lut',mole_frac_lut(im1:im2)
!    endif
    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))
!      if(iac.eq.1.and.iam.eq.0.and.k>=9200.and.k<=9300)then
!        write(6,*)k,'im1',im1,'im2',im2,qext_so4(im1:im2)
!        write(6,*)'mole_frac_lut',mole_frac_lut(im1:im2)
        !write(6,*)'mole_frac',mole_frac(k,iac)
!        write(6,*)'reff_so4',reff_so4(im1:im2)
        !write(6,*)'qext',qext_molef(k,iac),reff_molef(k,iac)
    
!      endif
!     if(iac.eq.1.and.rh(k)>90.)then
!       write(500+mype,*)'k',k,rh(k),'im1',im1,reff_molef(k,iac),'im2',im2,'reff_so4',reff_so4(im1:im2)
!       write(500+mype,*)' qext ',qext_molef(k,iac),' qext2o4',qext_so4(im1:im2)
!     endif
!     if(iac.eq.1)then
!       write(500+mype,*)'k',k,'im1',im1,reff_molef(k,iac),'im2',im2,'reff_so4',reff_so4(im1:im2)
!     endif
    elseif(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))
    elseif(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))
    elseif(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))


  KEaod(k-1,iac)=mass_ext(k,iac)*(1.+Hterm)
!  if(iac.eq.1.and.iam.eq.0)then
!     !write(6,*)'sulf KEAOD ',k-1,KEaod(k-1,iac),'mass ',mass_ext(k,iac),' Hterm ',Hterm
!     write(6,*)'qext ',qext_molef(k,iac), ' den ',den_mix_rh(k,iac),' reff ',reff_molef(k,iac)
!     write(6,*)'im1',im1,' im2 ',im2
!   endif
!    if(iac.eq.1.and.rh(k)>90.)then
!      write(500+mype,*)'mass -ext ', mass_ext(k,iac),' den ',den_mix_rh(k,iac)
!      write(500+mype,*)' H term ',Hterm ,'reff ',reff_molef(k,iac)
!      write(500+mype,*)' KEaod ',KEaod(k-1,iac)
!      !write(500+mype,*)'mw_w/mw_air',mw_w/mw_air,' molef_rh ',molef_rh(k,iac)
!      write(500+mype,*)' frac ',(1.-molef_rh(k,iac))/molef_rh(k,iac)
!      write(500+mype,*)'irh1 ',irh1,irh2
!    endif
  reff(k-1,iac)=reff_molef(k,iac)
!  write(500+mype,*)'reff',reff(k-1,iac),'k-1',k-1,'iac',iac,'KE',Keaod(k-1,iac)
!  call flush(500+mype)
!  if(mype.eq.0)then
!     write(6,*)'reff ',k-1,'iac',iac,'reff',reff(k-1,iac),'KE',keaod(k-1,iac)
!     call flush(6)
!  endif


 enddo   ! end loop over layers

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

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



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


 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

!-------------------------------------------!
!  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,nrhi
  ! Find two closest relative humidities from rh_lut
!  !if(rh(k)<2..and.iac.eq.7)then
!   write(500+mype,*)'rh_hlut',rh_hlut

!  endif
! handle dust first since does not depend on rh
  if(iac.ge.7.and.iac.le.10)then
    q_rh_hlut(k,iac)=q_tmp_hlut(1)
    reff_rh_hlut(k,iac)=reff_tmp_hlut(1)
    Gterm=1.0 ! since f_w=1.0
    den_aerow=den_aero(iac)
  elseif(iac>=11)then
!   handle seasalt
    do irh=1,7
!   if(k<10.and.iam.eq.0)then
!     write(6,*)'rh k ',rh(k),'rh_hlut',irh,rh_hlut(irh)
!   endif
     if(rh(k).le.rh_hlut(irh)) goto 8
    enddo
    irh=7
8   continue
    if(irh.eq.1)then
      q_rh_hlut(k,iac)=q_tmp_hlut(1)
      reff_rh_hlut(k,iac)=reff_tmp_hlut(1)
    else
      irh1=irh-1
      irh2=irh
      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)




      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)
    endif
! clips reff for sea salt at 1.e-10
    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

!  Gterm=1./f_w*den_aerow/den_aero(iac)
    Gterm=1./f_w*den_aerow/den_aero(iac)
  endif

  ! beta: specific extinction (m2/kg)

  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.
  if(iac.ge.7.and.iac.le.14) then 
    if(iac.le.10)then
      KEaod(k-1,iac-1)=beta_aero_hlut(k,iac)*(1.+Gterm)
      reff(k-1,iac-1)=reff_rh_hlut(k,iac)
    else
      KEaod(k-1,iac)=beta_aero_hlut(k,iac)*(1.+Gterm)
      reff(k-1,iac)=reff_rh_hlut(k,iac)
    endif
  endif

 enddo   ! end loop over layers



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


first=.false.

return
end subroutine compute_aod_KE_raqms 
#ifdef ADDCODEGSI
      subroutine compute_aod_raqmsgsi(nchem,rh,dzaod,chemp,kgkg_kgm2,aod,ext_3d,aodpart)
      use pmgrid, only : nc,beglat,endlat,plev,iam,masterproc
      use mpishorthand
      use computeaodraqms, only : KEaod,rrmaxrhraq,nrhaod,KEcrtmrh
      implicit none
      integer nac_hlut,i,j,k,m
      parameter (nac_hlut=14)
      real aod(nc,beglat:endlat),aodpart(nc,beglat:endlat,nac_hlut),ext_3d(nc,beglat:endlat,plev,nac_hlut)
      real aero_conc(nc,beglat:endlat,plev,nac_hlut),kgkg_kgm2(nc,beglat:endlat,plev)
      real dzaod(nc,beglat:endlat,plev)
      real mw(14),mwair
      data mw/132.,12.,12.,16.8,16.8,9*28.97/
      data mwair/28.97/

      aodpart=0.0
      do m=1,14
        do k=1,plev
          do lat=beglat,endlat
            do i=1,nc
              aero_conc(i,lat,k,m)=chemp(i,lat,k,m)*mw(m)/mwair*kgkg_kgm2(i,j,k)
              aodpart(i,lat,m)=aodpart(i,lat,m)+KEcrtm*aero_conc(i,lat,k,m)
              ext_3d(i,lat,k,m)=KEcrtm*aero_conc(i,lat,k,m)/dzaod(i,j,k)
            end do
          end do
        end do
        aod=aod+aodpart(:,:,m)
      end do
      return
      end subroutine compute_aod_raqmsgsi
#endif
