   subroutine setupaod(lunin,mype,nchanl,nreal,nobs,&
     obstype,isis,is,aero_diagsave,init_pass)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    setupaod    compute rhs of oi equation for aod
!   prgmmr: hclin            org: ncar/mmm                date: 2010-10-20
!
! abstract: read in data, first guess, and obtain rhs of oi equation
!        for aod.
!
! program history log:
!   2010-10-20  hclin - modified from setuprad for aod
!
!  input argument list:
!     lunin   - unit from which to read radiance (brightness temperature, tb) obs
!     mype    - mpi task id
!     nchanl  - number of channels per obs
!     nreal   - number of pieces of non-tb information per obs
!     nobs    - number of tb observations to process
!     obstype - type of tb observation
!     isis    - sensor/instrument/satellite id  ex.amsua_n15
!     is      - integer counter for number of observation types to process
!     aero_diagsave - logical to switch on diagnostic output (.false.=no output)
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

  use radinfo, only: nsigradjac
  use aeroinfo, only: nsigaerojac
  use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex, &
      isatid,itime,ilon,ilat,iszen_ang,isazi_ang
!  use CRTM_Aerosolscatter, only : crtm_ke
  use mpeu_util, only: die,perr
  use kinds, only: r_kind,r_single,i_kind,r_double
  use crtm_spccoeff, only: sc
  use obsmod, only: ianldate,iadate,ndat,mype_diaghdr,nchan_total, &
           dplat,obsdiags,obsptr,lobsdiagsave,lobsdiag_allocated,&
           dirname,time_offset
  use obsmod, only: obs_diag
  use gsi_4dvar, only: nobs_bins,hr_obsbin
  use gridmod, only: nsig,regional,msig,get_ij
  USE constants, ONLY: tiny_r_kind,zero,one,three,r10
  use jfunc, only: jiter,miter
  use m_dtime, only: dtime_setup, dtime_check, dtime_show
  use chemmod, only: laeroana_gocart, l_aoderr_table
  use aeroinfo, only: jpch_aero, nusis_aero, nuchan_aero, iuse_aero, &
       error_aero, gross_aero
  use obsmod, only: i_aero_ob_type, aerohead, aerotail, aero_ob_type
  use obsmod, only: rmiss_single
  use qcmod, only: ifail_crtm_qc
!  use raqmsmod, only : aodobsgrd,wtadobsgrd
!  use crtmKEmod, only : ke

  implicit none

! Declare passed variables
  logical                           ,intent(in   ) :: aero_diagsave
  character(10)                     ,intent(in   ) :: obstype
  character(20)                     ,intent(in   ) :: isis
  integer(i_kind)                   ,intent(in   ) :: lunin,mype,nchanl,nreal,nobs,is
  logical                           ,intent(in   ) :: init_pass  ! state of "setup" processing

! Declare external calls for code analysis
  external:: stop2

! Declare local parameters
  integer(i_kind),parameter:: ipchan=4
  integer(i_kind),parameter:: ireal=5

  real(r_kind),parameter:: r1e10=1.0e10_r_kind

! Declare local variables
  character(128) diag_aero_file

  integer(i_kind) error_status,istat
  integer(i_kind) m,jc
  integer(i_kind) icc
  integer(i_kind) j,k,ncnt,i
  integer(i_kind) mm1
  integer(i_kind) kk,n,kk2,ibin,ioff,iii
  integer(i_kind) ii,jj,idiag

  real(r_single) freq4,pol4,wave4,varch4
  real(r_kind) obserr,errinv,useflag
  real(r_kind) trop5,pangs
  real(r_kind) cenlon,cenlat,slats,slons,dtime
  real(r_kind) val_obs

! Declare local arrays

  real(r_single),dimension(ireal):: diagbuf
  real(r_single),allocatable,dimension(:,:):: diagbufchan

  real(r_kind),dimension(nchanl):: varinv,varinv_use,error0
  real(r_kind),dimension(nchanl):: tnoise,errmax
  real(r_kind),dimension(nchanl):: var,ratio_aoderr,aodinv
  real(r_kind),dimension(nreal+nchanl,nobs)::data_s
  real(r_kind),dimension(nsig):: prsltmp
  real(r_kind),dimension(nsig):: qvp,tvp
  real(r_kind),dimension(nsig+1):: prsitmp
  real(r_kind) dtsavg
  real(r_single) :: psfc

  integer(i_kind),dimension(nchanl):: ich,id_qc
  integer(i_kind),dimension(msig):: klevel

  character(10) filex
  character(12) string

  logical toss,l_may_be_passive
  logical,dimension(nobs):: luse

  logical:: in_curbin, in_anybin
  integer(i_kind),dimension(nobs_bins) :: n_alloc
  integer(i_kind),dimension(nobs_bins) :: m_alloc
  type(aero_ob_type),pointer:: my_head
  type(obs_diag),pointer:: my_diag
  character(len=*),parameter:: myname="setupaod"

  real(r_kind), dimension(nchanl) :: total_aod, aod_obs, aod

  integer(i_kind) :: istyp, idbcf, ilone, ilate,kmax,km,nobsused
  real(r_kind)    :: styp, dbcf

  real(r_kind),dimension(nchanl):: emissivity,ts,emissivity_k
  real(r_kind),dimension(nchanl):: tsim
  real(r_kind),dimension(nsig,nchanl):: wmix,temp,ptau5
  real(r_kind),dimension(nsigradjac,nchanl):: jacobian
  real(r_kind),dimension(nsigaerojac,nchanl):: jacobian_aero
  real(r_kind),dimension(nsigaerojac,nchanl):: jacobian_aero_adjoint
!  real(r_kind),dimension(nsig,14) :: aodfrac
  real(r_kind),dimension(nsig,nchanl):: layer_od
  real(r_kind) :: clw_guess, tzbgr, sfc_speed
  real(r_kind) :: sum_jacobian_aero
  logical ldiag
  common /dodiag/ldiag
  kmax=-1
  nobsused=0
!  write(6,*)mype,'airsco ajl new code top setupaod crtm_ke',mype,laeroana_gocart 
  if(mype.eq.0)then
  write(6,*)'raqmsaod nreal',nreal,'nchanl',nchanl
   endif
!  write(6,*)'r_kind',r_kind,'r_single',r_single,'r_double',r_double
!  call flush(6)
!  write(300+mype,*)'top setupaod crtm_ke'
!  call flush(300+mype)

  if ( .not. laeroana_gocart ) then
     return
  endif
!  if(.not.allocated(crtm_ke))then
!    allocate (crtm_ke(nsig,14))
!!endif

  n_alloc(:)=0
  m_alloc(:)=0
!**************************************************************************************
! Initialize variables and constants.
  mm1        = mype+1
  ncnt       = 0
  icc   = 0

  isatid    = 1  ! index of satellite id
  itime     = 2  ! index of analysis relative obs time 
  ilon      = 3  ! index of grid relative obs location (x)
  ilat      = 4  ! index of grid relative obs location (y)
  ilone     = 5  ! index of earth relative longitude (degrees)
  ilate     = 6  ! index of earth relative latitude (degrees)
  iszen_ang = 8  ! index of solar zenith angle (degrees)
  isazi_ang = 9  ! index of solar azimuth angle (degrees)
  istyp     = 10 ! index of surface type
  idbcf     = 11 ! index of deep blue confidence flag

! Initialize channel related information
  tnoise = r1e10
  errmax = r1e10
  l_may_be_passive = .false.
  toss = .true.
  jc=0

  do j=1,jpch_aero
     if(trim(isis) == trim(nusis_aero(j))) then 
        jc=jc+1
        if(jc > nchanl)then
           write(6,*)'setupaod:  ***ERROR*** in channel numbers, jc,nchanl=',jc,nchanl,&
                '  ***STOP IN setupaod***'
           call stop2(71)
        end if

!       Load channel numbers into local array based on satellite type

        ich(jc)=j
!
!       Set error instrument channels
        tnoise(jc)=error_aero(j)
        errmax(jc)=gross_aero(j)
        if (iuse_aero(j)< -1 .or. (iuse_aero(j) == -1 .and.  &
              .not.aero_diagsave)) tnoise(jc)=r1e10
        if (iuse_aero(j)>-1) l_may_be_passive=.true.
        if (tnoise(jc) < 1.e4_r_kind) toss = .false.
     end if
  end do
  if ( mype == 0 .and. .not.l_may_be_passive) write(6,*)mype,'setupaod: passive obs',is,isis
  if(nchanl > jc) write(6,*)'setupaod:  channel number reduced for ', &
       obstype,nchanl,' --> ',jc
  if(jc == 0) then
     if(mype == 0) write(6,*)'setupaod: No channels found for ', &
          obstype,isis
     if(nobs > 0)read(lunin)
     return
  end if
  if (toss) then
     if(mype == 0)write(6,*)'setupaod: all obs var > 1e4.  do not use ',&
          'data from satellite is=',isis
     if(nobs >0)read(lunin)                    
     return
  endif

! Initialize radiative transfer
  call init_crtm(mype_diaghdr(is),mype,nchanl,isis,obstype)

! If diagnostic file requested, open unit to file and write header.
  if (aero_diagsave) then
     filex=obstype
     write(string,1976) jiter
1976 format('_',i2.2)
     diag_aero_file= trim(dirname) // trim(filex) // '_' // trim(dplat(is)) // trim(string)
     if(init_pass) then
        open(4,file=trim(diag_aero_file),form='unformatted',status='unknown',position='rewind')
     else
        open(4,file=trim(diag_aero_file),form='unformatted',status='old',position='append')
     endif

!    Initialize/write parameters for satellite diagnostic file on
!    first outer iteration.
     if (init_pass .and. mype==mype_diaghdr(is)) then
        write(4) isis,dplat(is),obstype,jiter,nchanl,ianldate,ireal,ipchan,nsig
        write(6,*)'setupaod:  write header record for ',&
             isis,ireal,' to file ',trim(diag_aero_file),' ',ianldate

        do i=1,nchanl
           n=ich(i)
           if( n < 1 )cycle
           varch4=error_aero(n)
           freq4=sc(sensorindex)%frequency(i)
           pol4=sc(sensorindex)%polarization(i)
           wave4=sc(sensorindex)%wavenumber(i)
           write(4)freq4,pol4,wave4,varch4,iuse_aero(n),&
                nuchan_aero(n),ich(i)
        end do
     endif
  endif

  idiag=ipchan
  if (lobsdiagsave) idiag=idiag+4*miter+1
  allocate(diagbufchan(idiag,nchanl))
  
! Load data array for current satellite
!  call wallt('data_s',0)
  read(lunin) data_s,luse
!  call wallt('data_s',1)

!  write(*,*)mype,'read in data',nobs
!  call flush(6)
! Loop over data in this block
  call dtime_setup()
! call wallt('nobs',0)
  do n = 1,nobs
!    Extract analysis relative observation time.
     dtime = data_s(itime,n)
     call dtime_check(dtime, in_curbin, in_anybin)
     if(.not.in_anybin) cycle
!     write(6,*)mype,' n ',n
!     call flush(6)

     if(in_curbin) then

        id_qc = 0

!       Extract lon and lat.
        slons  = data_s(ilon,n)    ! grid relative longitude
        slats  = data_s(ilat,n)    ! grid relative latitude                     
        cenlon = data_s(ilone,n)   ! earth relative longitude (degrees)
        cenlat = data_s(ilate,n)   ! earth relative latitude (degrees)                       
!              if(cenlat<-51.)then
!               write(300+mype,*)'top cenlat',cenlat,cenlon
!              endif

        pangs  = data_s(iszen_ang,n)
        styp   = data_s(istyp,n)
        dbcf   = data_s(idbcf,n)
 
!       Set relative weight value
        val_obs=one

!       Load channel data into work array.
        aod_obs = rmiss_single
        do i = 1, nchanl
           aod_obs(i) = data_s(i+nreal,n)
        end do
!        if(n.eq.1)write(6,*)'nchanl',nchanl
!        if(n<10)write(6,*)mype,'lon',cenlon,cenlat,'aod',aod_obs(4),'l_aoderr_table',l_aoderr_table, &
!        'styp',styp,'dbcf',dbcf

        if ( .not. l_aoderr_table ) then
!          set observation error
           select case ( nint(styp) )
              case ( 0 )        ! water
                 tnoise = 0.03_r_kind+0.05_r_kind*aod_obs
              case ( 1, 2, 3 )  ! coast, desert, land
                 tnoise = 0.05_r_kind+0.15_r_kind*aod_obs
              case ( 4 )        ! deep blue
                 if ( nint(dbcf) >= 0 .and. nint(dbcf) <= 3 ) then
                    tnoise = 0.05_r_kind+0.15_r_kind*aod_obs+0.01_r_kind*(three-dbcf)
                 end if
              case ( 5 )  ! nnr ocean
                 tnoise = 0.2_r_kind*(aod_obs+0.01)
              case ( 6 )  ! nnr land
                 tnoise = 0.2_r_kind*(aod_obs+0.01)
           end select
        end if
 
!       Interpolate model fields to observation location, call ctm and create jacobians
!        write(6,*)mype,'call crtm'
!        call flush(6)
!        write(6,*)mype,'zero crtm_ke'
!        call flush(6)
!        crtm_ke=0.0
         if(n>5)then
         ldiag=.false.
         else

        ldiag=.true.
        endif
        call call_crtm(obstype,dtime,data_s(1,n),nchanl,nreal,ich, &
             tvp,qvp,clw_guess,prsltmp,prsitmp, &
             trop5,tzbgr,dtsavg,sfc_speed, &
             tsim,emissivity,ptau5,ts,emissivity_k, &
             temp,wmix,jacobian,error_status,layer_od,jacobian_aero,jacobian_aero_adjoint)
!             temp,wmix,jacobian,error_status,layer_od,jacobian_aero,jacobian_aero_adjoint,aodfrac)
!        write(6,*)'did call crtm'
!        call flush(6)


! If the CRTM returns an error flag, do not assimilate any channels for this ob
! and set the QC flag to ifail_crtm_qc.
! We currently go through the rest of the QC steps, ensuring that the diagnostic
! files are populated, but this could be changed if it causes problems.
        if (error_status /=0) then
           id_qc(1:nchanl) = ifail_crtm_qc
           varinv(1:nchanl) = zero
        endif

        total_aod = zero
        do i = 1, nchanl
           total_aod(i) =sum(layer_od(:,i))
        enddo 
!        write(600+mype,*)'total_aod',total_aod(4),'aod_obs',aod_obs(4)
!        do k=nsig,1,-1
!          write(600+mype,*)'layer_od',k,layer_od(k,4)
!        end do

        do i = 1, nchanl
           aod(i) = aod_obs(i) - total_aod(i)
           error0(i)     = tnoise(i)
           if(aod_obs(i)>zero .and. tnoise(i) < 1.e4_r_kind .or. (iuse_aero(ich(i))==-1  &
              .and. aero_diagsave))then
              varinv(i)     = val_obs/tnoise(i)**2
!              if(cenlat<-51.)then
!                   write(300+mype,*)'south ',cenlat,cenlon,'aodinv',aod(i),'obs',aod_obs(i),'total_aod',total_aod(i)
!                   write(300+mype,*)'tnoise',tnoise(i)
!              endif
           else
              if(id_qc(i) == 0)id_qc(i)=1
              varinv(i)     = zero
!              if(cenlat<-51..and.i.eq.4)then
!                write(300+mype,*)'thow out ',cenlat,cenlon,'tnoise',tnoise(i),'aod_obs',aod_obs(i)
!              endif
           endif
        end do

        !do i = 1, nchanl
        !   if ( aod_obs(i) > zero ) then
        !     write(6,'(A,3i6,4f8.3,2f8.2)') 'mype, iobs, ichan, aod_crtm, aod_obs, omb, err, lat, lon : ',  &
        !         mype, n, i, total_aod(i), aod_obs(i),aod(i), tnoise(i), cenlat, cenlon
        !   end if
        !end do

        icc = 0
        do i = 1, nchanl
           ! Only process observations to be assimilated
           if (varinv(i) > tiny_r_kind ) then
               m = ich(i)
               ! Only "good" obs are included in J calculation.
               if (iuse_aero(m) >= 1)then
!           write(300+mype,*)'aod ',i,aod(i),' obs ',aod_obs(i),' total',total_aod(i)
!           call flush(300+mype)
                  icc = icc + 1
                  aodinv(icc) = aod(i)            ! obs-ges innovation
                  var(icc) = one/error0(i)**2     ! 1/(obs error)**2  (original uninflated error)
                  ratio_aoderr(icc)=error0(i)**2*varinv(i) ! (original error)/(inflated error)
               endif
           endif
        end do
     endif ! (in_curbin)
!     write(6,*)'at zz'
     !call flush(6)

!    In principle, we want ALL obs in the diagnostics structure but for
!    passive obs (monitoring), it is difficult to do if aero_diagsave
!    is not on in the first outer loop. For now we use l_may_be_passive...
!     write(300+mype,*)'l_may_be_passive',l_may_be_passive,'icc',icc,'in_curbin',in_curbin,'n',n
!     call flush(300+mype)
     if (l_may_be_passive) then
!       Link observation to appropriate observation bin
        if (nobs_bins>1) then
           ibin = NINT( dtime/hr_obsbin ) + 1
        else
           ibin = 1
        endif
        if (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin

        if (in_curbin) then
!          Load data into output arrays
           if (icc > 0) then
              ncnt =ncnt+1
              nchan_total=nchan_total+icc
!              write(440+mype,*)' add head',ncnt,' icc ',icc
!              call flush(440+mype)

              if (.not. associated(aerohead(ibin)%head))then
                 allocate(aerohead(ibin)%head,stat=istat)
                 if (istat /= 0) write(6,*)' failure to write aerohead '
                 aerotail(ibin)%head => aerohead(ibin)%head
              else
                 allocate(aerotail(ibin)%head%llpoint,stat=istat)
                 if (istat /= 0) write(6,*)' failure to write aerotail%llpoint '
                 aerotail(ibin)%head => aerotail(ibin)%head%llpoint
              end if
 
              m_alloc(ibin) = m_alloc(ibin) +1
              my_head => aerotail(ibin)%head
              my_head%idv = is
              my_head%iob = n
 
              allocate(aerotail(ibin)%head%res(icc),aerotail(ibin)%head%err2(icc), &
                       aerotail(ibin)%head%diags(icc),&
                       aerotail(ibin)%head%raterr2(icc), &
                       aerotail(ibin)%head%daod_dvar(nsigaerojac,icc), &
                       aerotail(ibin)%head%daod_dvar_adjoint(nsigaerojac,icc), &
                       aerotail(ibin)%head%ich(icc),&
                       aerotail(ibin)%head%icx(icc))

              aerotail(ibin)%head%nlaero  = icc         ! profile observation count
              call get_ij(mm1,slats,slons,aerotail(ibin)%head%ij(1),aerotail(ibin)%head%wij(1))
                     aerotail(ibin)%head%time=dtime
                     aerotail(ibin)%head%luse=luse(n)
                     aerotail(ibin)%head%ich(:)=-1
              iii=0
              do ii=1,nchanl
                 m=ich(ii)
                 if (varinv(ii)>tiny_r_kind .and. iuse_aero(m)>=1) then
                    iii=iii+1
                    aerotail(ibin)%head%res(iii)=aodinv(iii)
                    aerotail(ibin)%head%err2(iii)=var(iii)
                    aerotail(ibin)%head%raterr2(iii)=ratio_aoderr(iii)
                    aerotail(ibin)%head%icx(iii)=m
                    do k = 1, nsigaerojac
                       aerotail(ibin)%head%daod_dvar(k,iii)=jacobian_aero(k,ii)
                       km=(k-1)/nsig+1
                       kk=k-(km-1)*nsig
!                       aerotail(ibin)%head%daod_dvar_aodfrac(k,iii)=jacobian_aero(k,ii)*aodfrac(kk,km)
                       aerotail(ibin)%head%daod_dvar_adjoint(k,iii)=jacobian_aero_adjoint(k,ii)
!                       if(abs(aerotail(ibin)%head%daod_dvar_aodfrac(k,iii)-jacobian_aero_adjoint(k,ii))>1.e-20)then
                         !write(1000+mype,*)mype,'jacobian ',k,'km ',km,'kk',kk,'iii',iii,aerotail(ibin)%head%daod_dvar_aodfrac(k,iii),jacobian_aero_adjoint(k,ii),' aodfrac ',aodfrac(kk,km),' jacobain_aero ',jacobian_aero(k,ii)
                         !call flush(1000+mype)
!                       endif
                       if(aerotail(ibin)%head%daod_dvar_adjoint(k,iii)>0.0)then
                         if(kk>kmax)kmax=kk
!                         write(300+mype,*)'dvar_aodfrac',kk,km,aerotail(ibin)%head%daod_dvar_aodfrac(k,iii),' inv ',aodinv(iii)
!                         call flush(300+mype)
                       endif
!                       if(jacobian_aero(k,ii)>0.0)then
!                         if(kk>kmax)kmax=kk
                       !else
!                         if(mype.eq.0.and.ii.eq.4.and.km.ne.10)then
!                           write(6,*)n,'aodinv',aodinv(iii),'jacob',k,km,jacobian_aero(k,ii),'var',var(iii)
!                         endif
!                       endif
                    end do
                    nobsused=nobsused+1
                    my_head%ich(iii)=ii
                 end if
              end do
!if(mype==0) then
!   DO ii=1,iii
!   DO k=1,nsigaerojac
!      write(*,*) ii,k,aerotail(ibin)%head%daod_dvar(k,ii)
!   ENDDO
!   ENDDO
!              write(440+mype,*)' did head and fill'
!              call flush(400+mype)
!endif
           end if ! icc
        endif ! (in_curbin)

!       Link obs to diagnostics structure
        iii=0
!        write(440+mype,*)' before allocates '
!        call flush(440+mype)
        do ii=1,nchanl
           if (.not.lobsdiag_allocated) then
              if (.not.associated(obsdiags(i_aero_ob_type,ibin)%head)) then
                 allocate(obsdiags(i_aero_ob_type,ibin)%head,stat=istat)
                 if (istat/=0) then
                    write(6,*)'setupaod: failure to allocate obsdiags',istat
                    call stop2(276)
                 end if
                 obsdiags(i_aero_ob_type,ibin)%tail => obsdiags(i_aero_ob_type,ibin)%head
              else
                 allocate(obsdiags(i_aero_ob_type,ibin)%tail%next,stat=istat)
                 if (istat/=0) then
                    write(6,*)'setupaod: failure to allocate obsdiags',istat
                    call stop2(277)
                 end if
                 obsdiags(i_aero_ob_type,ibin)%tail => obsdiags(i_aero_ob_type,ibin)%tail%next
              end if
              allocate(obsdiags(i_aero_ob_type,ibin)%tail%muse(miter+1))
              allocate(obsdiags(i_aero_ob_type,ibin)%tail%nldepart(miter+1))
              allocate(obsdiags(i_aero_ob_type,ibin)%tail%tldepart(miter))
              allocate(obsdiags(i_aero_ob_type,ibin)%tail%obssen(miter))
              obsdiags(i_aero_ob_type,ibin)%tail%indxglb=(n-1)*nchanl+ii
              obsdiags(i_aero_ob_type,ibin)%tail%nchnperobs=-99999
              obsdiags(i_aero_ob_type,ibin)%tail%luse=.false.
              obsdiags(i_aero_ob_type,ibin)%tail%muse(:)=.false.
              obsdiags(i_aero_ob_type,ibin)%tail%nldepart(:)=-huge(zero)
              obsdiags(i_aero_ob_type,ibin)%tail%tldepart(:)=zero
              obsdiags(i_aero_ob_type,ibin)%tail%wgtjo=-huge(zero)
              obsdiags(i_aero_ob_type,ibin)%tail%obssen(:)=zero
 
              n_alloc(ibin) = n_alloc(ibin) +1
              my_diag => obsdiags(i_aero_ob_type,ibin)%tail
              my_diag%idv = is
              my_diag%iob = n
              my_diag%ich = ii
           else
              if (.not.associated(obsdiags(i_aero_ob_type,ibin)%tail)) then
                 obsdiags(i_aero_ob_type,ibin)%tail => obsdiags(i_aero_ob_type,ibin)%head
              else
                 obsdiags(i_aero_ob_type,ibin)%tail => obsdiags(i_aero_ob_type,ibin)%tail%next
              end if
              if (obsdiags(i_aero_ob_type,ibin)%tail%indxglb/=(n-1)*nchanl+ii) then
                 write(6,*)'setupaod: index error'
                 call stop2(278)
              endif
           endif

           if (in_curbin) then
              if (ii==1) obsptr => obsdiags(i_aero_ob_type,ibin)%tail
              if (ii==1) obsdiags(i_aero_ob_type,ibin)%tail%nchnperobs = nchanl
              obsdiags(i_aero_ob_type,ibin)%tail%luse = luse(n)
              obsdiags(i_aero_ob_type,ibin)%tail%nldepart(jiter) = aod(ii)
              obsdiags(i_aero_ob_type,ibin)%tail%wgtjo=varinv(ii)
!              if(ii.eq.4)then
!                call aodobsgrid(aod_obs(ii),slats,slons)
!              endif
 
!             Load data into output arrays
              m=ich(ii)
              if (varinv(ii)>tiny_r_kind .and. iuse_aero(m)>=1) then
                 iii=iii+1
                 aerotail(ibin)%head%diags(iii)%ptr => obsdiags(i_aero_ob_type,ibin)%tail
                 obsdiags(i_aero_ob_type,ibin)%tail%muse(jiter) = .true.
  
                 ! verify the pointer to obsdiags
 
                 my_head => aerotail(ibin)%head
                 my_diag => aerotail(ibin)%head%diags(iii)%ptr
 
                 if (my_head%idv      /= my_diag%idv .or. &
                     my_head%iob      /= my_diag%iob .or. &
                     my_head%ich(iii) /= my_diag%ich ) then
                    call perr(myname,'mismatching %[head,diags]%(idv,iob,ich,ibin) =', &
                          (/is,i,ii,ibin/))
                    call perr(myname,'my_head%(idv,iob,ich) =',(/my_head%idv,my_head%iob,my_head%ich(iii)/))
                    call perr(myname,'my_diag%(idv,iob,ich) =',(/my_diag%idv,my_diag%iob,my_diag%ich/))
                    call die(myname)
                 endif
              endif
           endif ! (in_curbin)
        enddo
        if (in_curbin) then
           if( iii/=icc ) then
              write(6,*)'setupaod: error iii icc',iii,icc
              call stop2(279)
           endif
        endif ! (in_curbin)
 
!    End of l_may_be_passive block
     endif

     if(in_curbin) then
!       Write diagnostics to output file.
        if (aero_diagsave .and. luse(n)) then
           diagbuf(1)  = cenlat                         ! observation latitude (degrees)
           diagbuf(2)  = cenlon                         ! observation longitude (degrees)
           diagbuf(3)  = dtime-time_offset              ! observation time (hours relative to analysis time)
           diagbuf(4)  = pangs                          ! solar zenith angle (degrees)
           diagbuf(5)  = data_s(isazi_ang,n)            ! solar azimuth angle (degrees)
 
           do i=1,nchanl
              diagbufchan(1,i)=aod_obs(i)      ! observed brightness temperature (K)
              diagbufchan(2,i)=aod(i)    ! observed - simulated Tb with no bias corrrection (K)
              errinv = sqrt(varinv(i))
              diagbufchan(3,i)=errinv          ! inverse observation error
              useflag=one
              if (iuse_aero(ich(i)) < 1) useflag=-one
              diagbufchan(4,i)= id_qc(i)*useflag! quality control mark or event indicator
           end do

           if (lobsdiagsave) then
              if (l_may_be_passive) then
                 do ii=1,nchanl
                    if (.not.associated(obsptr)) then
                       write(6,*)'setupaod: error obsptr'
                       call stop2(280)
                    end if
                    if (obsptr%indxglb/=(n-1)*nchanl+ii) then
                       write(6,*)'setupaod: error writing diagnostics'
                       call stop2(281)
                    end if

                    ioff=4
                    do jj=1,miter
                       ioff=ioff+1
                       if (obsptr%muse(jj)) then
                          diagbufchan(ioff,ii) = one
                       else
                          diagbufchan(ioff,ii) = -one
                       endif
                    enddo
                    do jj=1,miter+1
                       ioff=ioff+1
                       diagbufchan(ioff,ii) = obsptr%nldepart(jj)
                    enddo
                    do jj=1,miter
                       ioff=ioff+1
                       diagbufchan(ioff,ii) = obsptr%tldepart(jj)
                    enddo
                    do jj=1,miter
                       ioff=ioff+1
                       diagbufchan(ioff,ii) = obsptr%obssen(jj)
                    enddo
 
                    obsptr => obsptr%next
                 enddo
              else
                 ioff=5
                 diagbufchan(ioff+1:ioff+4*miter+1,1:nchanl) = zero
              endif
           endif
           psfc=prsitmp(1)*r10 ! convert tp hPa
           WRITE(4) psfc,diagbuf,diagbufchan

        end if
     endif ! (in_curbin)

100  continue

! End of n-loop over obs
  end do
! call wallt('nobs',1)

! Jump here when there is no data to process for current satellite
! Deallocate arrays
  deallocate(diagbufchan)

  if (aero_diagsave) then
     call dtime_show(myname,'diagsave:aero',i_aero_ob_type)
     close(4)
  endif

  call destroy_crtm

! End of routine
!  write(6,*)'bottom setup aod'
!  call flush(6)
  if(kmax.ne.35)then
  write(6,*)mype,'kmax',kmax,' nobs ',nobs,' nobsused ',nobsused
  call flush(6)
  endif
!  write(440+mype,*)'bottom of setupaod'
!  call flush(440+mype)

  return

end subroutine setupaod
subroutine calcaod(nchanl,nreali,obstypei,isisi,isi,icall)
  use raqmsmod, only : aodb,aodpartb
  use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex, &
      isatid,itime,ilon,ilat,iszen_ang,isazi_ang
  use mpeu_util, only: die,perr
  use kinds, only: r_kind,r_single,i_kind
  use crtm_spccoeff, only: sc
  use obsmod, only: ianldate,iadate,ndat,mype_diaghdr,nchan_total, &
           dplat,obsdiags,obsptr,lobsdiagsave,lobsdiag_allocated,&
           dirname,time_offset
  use mpimod,only : mype
  use aeroinfo, only: jpch_aero, nusis_aero, nuchan_aero, iuse_aero
  use gridmod, only : istart,jstart,ilat1,jlon1,nsig
  use radinfo, only: nsigradjac
  use aeroinfo, only: nsigaerojac
  implicit none
  character(10)                     ,intent(in   ) :: obstypei
  character(10)                      :: obstype
  character(20)                     ,intent(in   ) :: isisi
  character(20)                      :: isis
!  logical init_pass,final_pass,last
  integer(i_kind) :: nchanl,nreal,jc,j,i,n,is,ii,icall,ifill
  integer(i_kind) :: isi
  integer(i_kind) :: nreali
  integer(i_kind) error_status,istat,mm1
  real(r_kind), dimension(nchanl) :: total_aod
  real(r_kind) :: dtime
  real(r_kind),dimension(4) :: datagrid
  integer(i_kind),dimension(nchanl):: ich
  real(r_kind),dimension(nsig):: prsltmp
  real(r_kind),dimension(nsig):: qvp,tvp
  real(r_kind),dimension(nsig+1):: prsitmp
  real(r_kind),dimension(nchanl):: emissivity,ts,emissivity_k
  real(r_kind),dimension(nchanl):: tsim
  real(r_kind),dimension(nsig,nchanl):: wmix,temp,ptau5
  real(r_kind),dimension(nsigradjac,nchanl):: jacobian
  real(r_kind),dimension(nsigaerojac,nchanl):: jacobian_aero
  real(r_kind),dimension(nsig,nchanl):: layer_od
  real(r_kind) :: clw_guess, tzbgr, sfc_speed
  real(r_kind) trop5,pangs
  real(r_kind) dtsavg
! ajl add aodparts
!  real sumaod
  real :: aodaerosol(14),sumaerosol(14)
  common /doaodpart/aodaerosol,sumaerosol
  logical first
  logical ldiag
  common /dodiag/ldiag
  data first/.true./
  save first,is,isis,obstype
  if(.not.first.and.icall.eq.0)return
  if(first)then
    first=.false.
    is=isi
    isis=isisi
    obstype=obstypei
  endif
  if(icall.eq.0)then
    ifill=1
  elseif(icall.eq.1)then
    ifill=2
  endif
!  if(mype.eq.16)write(6,*)'calcaod init_pass',init_pass,final_pass,'last',last
  mm1=mype+1 
  if(.not.allocated(aodb))then
    allocate (aodb(ilat1(mm1),jlon1(mm1),2),aodpartb(ilat1(mm1),jlon1(mm1),14,2))
!    write(6,*)'mype',mype,'aodb',shape(aodb),'istart',istart(mm1),jstart(mm1)
!    call flush(6)
  endif
!  if(mype.eq.16)then
!    write(6,*)'calcaod nchanl',nchanl,nreal,'obstype',obstype,'isis',isis,'is',is
!    write(6,*)'jpch_aero',jpch_aero,'isis',isis
!  endif
! Initialize radiative transfer
  jc=0

  do j=1,jpch_aero
!     if(mype.eq.16)then
!       write(6,*)'nusis_aero',j,nusis_aero(j)
!     endif
     if(trim(isis) == trim(nusis_aero(j))) then 
        jc=jc+1
        if(jc > nchanl)then
           write(6,*)'setupaod:  ***ERROR*** in channel numbers, jc,nchanl=',jc,nchanl,&
                '  ***STOP IN setupaod***'
           call stop2(71)
        end if

!       Load channel numbers into local array based on satellite type

        ich(jc)=j
!
     end if
  end do
!  write(6,*)mype,'nchanl',nchanl,isis,obstype
  call init_crtm(mype_diaghdr(is),mype,nchanl,isis,obstype)
  dtime=0.
  nreal=4
  do j=1,jlon1(mm1)
    do i=1,ilat1(mm1)
      datagrid(ilat)=istart(mm1)+i-1 ! lat
      datagrid(ilon)=jstart(mm1)+j-1 ! lon

!      !if(mype.eq.16)then
!        write(6,*)'datagrid',datagrid,'ich',ich,'nchanl',nchanl,nreal
!        write(6,*)'obstype',obstype
!        call flush(6)
!      endif
      ldiag=.false.
      call call_crtm(obstype,dtime,datagrid,nchanl,nreal,ich, &
             tvp,qvp,clw_guess,prsltmp,prsitmp, &
             trop5,tzbgr,dtsavg,sfc_speed, &
             tsim,emissivity,ptau5,ts,emissivity_k, &
             temp,wmix,jacobian,error_status,layer_od,jacobian_aero)
             total_aod(4) =sum(layer_od(:,4))
             aodb(i,j,ifill)=total_aod(4)
!      if(datagrid(ilon).eq.1)then
!        write(6,*)mype,'jlon',datagrid(ilon),'ilat',datagrid(ilat),'aod',aodb(i,j),'i',i,j,'istart',istart(mm1),jstart(mm1)
!        call flush(6)
!      endif
             aodpartb(i,j,:,ifill)=aodaerosol(:)
!     if(mype.eq.16)then
!       write(6,*)mype,'aod',aodaerosol(:)
!       write(6,*)mype,'sum',sumaerosol(:)
!           call flush(6)
!          sumaod=0.

!          do ii=1,14
!            sumaod=sumaod+aodaerosol(ii)
!          end do
!          write(6,*)'sumaod',sumaod
!          write(6,*)'total aod',total_aod(4)
!      endif
     end do
   end do
!  if(mype.eq.16)then
!    write(6,*)'aodb',maxval(aodb),minval(aodb),'shape',shape(aodb)
!  endif
  call destroy_crtm
  return
end subroutine calcaod
