SUBROUTINE aerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
     pplay, paprs, t_seri, rhcl, &
     maerosol, maerosol_pi, &
     tau_aero, piz_aero, cg_aero )
  
  USE dimphy
  IMPLICIT NONE

! Input arguments
  LOGICAL, INTENT(IN)                      :: debut
  LOGICAL, INTENT(IN)                      :: new_aod
  INTEGER, INTENT(IN)                      :: flag_aerosol
  REAL, INTENT(IN)                         :: rjourvrai
  REAL, INTENT(IN)                         :: pdtphys
  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
  REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_seri
  REAL, DIMENSION(klon,klev), INTENT(IN)   :: rhcl   ! humidite relative ciel clair

! Output arguments
  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: maerosol
  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: maerosol_pi
  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero
  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero
  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero

! Local variables
  REAL, DIMENSION(klon)        :: aerindex ! POLDER aerosol index 
  REAL, DIMENSION(klon,10)     :: fractnat_allaer
  REAL, DIMENSION(klon,klev)   :: sulfate  ! SO4 aerosol concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: bcsol    ! BC soluble concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: bcins    ! BC insoluble concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: pomsol   ! POM soluble concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: pomins   ! POM insoluble concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: sulfate_pi
  REAL, DIMENSION(klon,klev)   :: bcsol_pi
  REAL, DIMENSION(klon,klev)   :: bcins_pi
  REAL, DIMENSION(klon,klev)   :: pomsol_pi
  REAL, DIMENSION(klon,klev)   :: pomins_pi
  REAL, DIMENSION(klon,klev)   :: pdel
  REAL, DIMENSION(klon,klev,8) :: m_allaer

  INTEGER :: k
  
!
! Read aersosols from file
!
  maerosol(:,:) = 0. 
  maerosol_pi(:,:) = 0. 

  bcsol(:,:) = 0. 
  bcins(:,:) = 0. 
  pomsol(:,:) = 0. 
  pomins(:,:) = 0. 
  sulfate(:,:) = 0. 

! Read sulfate
  IF ( flag_aerosol .EQ. 1 .OR. &
       flag_aerosol .EQ. 4 .OR. &
       flag_aerosol .EQ. 6 ) THEN 

     CALL readaerosol(5,rjourvrai, debut, sulfate)
     CALL readaerosol_preind(5,rjourvrai, &
          debut, sulfate_pi)

     maerosol(:,:) = maerosol(:,:) + sulfate(:,:) 
     maerosol_pi(:,:) = maerosol_pi(:,:) + sulfate_pi(:,:) 
  ENDIF


! Read bcscol and bcins
  IF ( flag_aerosol .EQ. 2 .OR. &
       flag_aerosol .EQ. 4 .OR. &
       flag_aerosol .EQ. 5 ) THEN 

     ! Get bc aerosol distribution 
     CALL readaerosol(3,rjourvrai, debut,bcsol )
     CALL readaerosol_preind(3,rjourvrai, debut, bcsol_pi)

     CALL readaerosol(7,rjourvrai, debut,bcins )
     CALL readaerosol_preind(7,rjourvrai, debut, bcins_pi)

     maerosol(:,:) = maerosol(:,:) + bcsol(:,:) 
     maerosol_pi(:,:) = maerosol_pi(:,:) + bcsol_pi(:,:) 
  ENDIF


! Read pomsol and pomins
  IF ( flag_aerosol .EQ. 3 .OR. &
       flag_aerosol .EQ. 4 .OR. &
       flag_aerosol .EQ. 5 .OR. &
       flag_aerosol .EQ. 6 ) THEN

     CALL readaerosol(4,rjourvrai, debut,pomsol)
     CALL readaerosol_preind(4,rjourvrai,debut,pomsol_pi )

     CALL readaerosol(8,rjourvrai, debut,pomins)
     CALL readaerosol_preind(8,rjourvrai,debut,pomins_pi )

     maerosol(:,:) = maerosol(:,:) + pomsol
     maerosol_pi(:,:) = maerosol_pi(:,:) + pomsol_pi
  ENDIF


! Save all aerosols in one variable
!
! ACo pour couplage aerosol offline 07/04/2009
! Tableau contenant les masses pour tous les aerosols 
! les valeurs a zero seront a remplacer par les bons
! tableaux lorsque les routines de lectures seront
! ajoutees. 
  m_allaer(:,:,1) = 0.                ! SSSSM || CSSSM
  m_allaer(:,:,2) = 0.                ! ASSSM
  m_allaer(:,:,3) = bcsol(:,:)        ! ASBCM
  m_allaer(:,:,4) = pomsol(:,:)       ! ASPOMM
  m_allaer(:,:,5) = sulfate(:,:)      ! ASSO4M || CSSO4M   
  m_allaer(:,:,6) = 0.                ! CIDUSTM
  m_allaer(:,:,7) = bcins(:,:)        ! AIBCM
  m_allaer(:,:,8) = pomins(:,:)       ! AIPOMM

!
! Calculate the optical properties for the aerosols
!
  DO k = 1, klev
     pdel(:,k) = paprs(:,k) - paprs (:,k+1)
  END DO

  IF (.NOT. new_aod) THEN 
     CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
          tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex)
  ELSE
     fractnat_allaer(:,:) = 0.
     CALL aeropt_2bands( &
          pdel, m_allaer, pdtphys, rhcl, & 
          tau_aero, piz_aero, cg_aero,   &
          fractnat_allaer, flag_aerosol, &
          pplay, t_seri) 

     CALL aeropt_5wv( &
          pdel, m_allaer, &
          pdtphys, rhcl, aerindex, & 
          flag_aerosol, pplay, t_seri)
  ENDIF

END SUBROUTINE aerosol_optic