! $Id$
!
SUBROUTINE readaerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
     pplay, paprs, t_seri, rhcl, &
     mass_ins_aero, mass_ins_aero_pi, &
     tau_aero, piz_aero, cg_aero )

! This routine will :
! 1) recevie the aerosols(already read and interpolated) corresponding to flag_aerosol
! 2) calculate the optical properties for the aerosols
!
  
  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)     :: mass_ins_aero    ! Total mass for all indissoluble aerosols
  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_ins_aero_pi !     -"-     preindustrial values
  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero     ! asymmetry parameter aerosol

! 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, i
  
!****************************************************************************************
! 1) Get aerosol mass
!    
!****************************************************************************************
! Read and interpolate sulfate
  IF ( flag_aerosol .EQ. 1 .OR. &
       flag_aerosol .EQ. 4 .OR. &
       flag_aerosol .EQ. 6 ) THEN 

     CALL readaerosol_interp(5, rjourvrai, debut, pplay, paprs, sulfate, sulfate_pi)
  ELSE
     sulfate(:,:) = 0. ; sulfate_pi(:,:) = 0.
  END IF

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

     ! Get bc aerosol distribution 
     CALL readaerosol_interp(3, rjourvrai, debut, pplay, paprs, bcsol, bcsol_pi )
     CALL readaerosol_interp(7, rjourvrai, debut, pplay, paprs, bcins, bcins_pi )
  ELSE
     bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
     bcins(:,:) = 0. ; bcins_pi(:,:) = 0.
  END IF


! Read and interpolate 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_interp(4, rjourvrai, debut, pplay, paprs, pomsol, pomsol_pi)
     CALL readaerosol_interp(8, rjourvrai, debut, pplay, paprs, pomins, pomins_pi)
  ELSE
     pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
     pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
  END IF


! Store 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 ! Coarse Soluble Sea Salt Mass
  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        ! Coarse Insoluble DUST Mass
  m_allaer(:,:,7) = bcins(:,:)        ! AIBCM
  m_allaer(:,:,8) = pomins(:,:)       ! AIPOMM

!
! Calculate the total mass of all indissoluble aersosols
!
  mass_ins_aero(:,:)    = sulfate(:,:)    + bcsol(:,:)    + pomsol(:,:) 
  mass_ins_aero_pi(:,:) = sulfate_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:) 

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

  IF (new_aod) THEN 

     fractnat_allaer(:,:) = 0.
     CALL aeropt_2bands( &
          pdel, m_allaer, pdtphys, rhcl, & 
          tau_aero, piz_aero, cg_aero,   &
          fractnat_allaer, flag_aerosol, &
          pplay, t_seri) 
     
     ! aeropt_5wv only for validation and diagnostics.
     ! In this version no diagnostics are set. 
     ! jg : may be desactivated if no diagnostics added.
     CALL aeropt_5wv( &
          pdel, m_allaer, &
          pdtphys, rhcl, aerindex, & 
          flag_aerosol, pplay, t_seri)
  ELSE

     CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
          tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex)
     
  END IF

END SUBROUTINE readaerosol_optic