! $Id$
!
SUBROUTINE readaerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
     pplay, paprs, t_seri, rhcl, presnivs, &
     mass_solu_aero, mass_solu_aero_pi, &
     tau_aero, piz_aero, cg_aero, &
     tausum_aero, tau3d_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
  USE aero_mod
  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
  REAL, DIMENSION(klev), INTENT(IN)        :: presnivs

! Output arguments
!****************************************************************************************
  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero    ! Total mass for all soluble aerosols
  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_solu_aero_pi !     -"-     preindustrial values
  REAL, DIMENSION(klon,klev,naero_grp,nbands), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
  REAL, DIMENSION(klon,klev,naero_grp,nbands), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
  REAL, DIMENSION(klon,klev,naero_grp,nbands), INTENT(OUT) :: cg_aero     ! asymmetry parameter aerosol
  REAL, DIMENSION(klon,nwave,naero_spc), INTENT(OUT)       :: tausum_aero
  REAL, DIMENSION(klon,klev,nwave,naero_spc), INTENT(OUT)  :: tau3d_aero

! Local variables
!****************************************************************************************
  REAL, DIMENSION(klon)        :: aerindex ! POLDER aerosol index 
  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)   :: cidust    ! DUST aerosol concentration  [ug/m3]
  REAL, DIMENSION(klon,klev)   :: sscoarse  ! SS Coarse concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: sssupco   ! SS Super Coarse concentration [ug/m3]
  REAL, DIMENSION(klon,klev)   :: ssacu     ! SS Acumulation 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)   :: cidust_pi
  REAL, DIMENSION(klon,klev)   :: sscoarse_pi
  REAL, DIMENSION(klon,klev)   :: sssupco_pi
  REAL, DIMENSION(klon,klev)   :: ssacu_pi
  REAL, DIMENSION(klon,klev)   :: pdel
  REAL, DIMENSION(klon,klev,naero_spc) :: m_allaer
  REAL, DIMENSION(klon,naero_tot)      :: fractnat_allaer

  INTEGER :: k, i
  
!****************************************************************************************
! 1) Get aerosol mass
!    
!****************************************************************************************
! Read and interpolate sulfate
  IF ( flag_aerosol .EQ. 1 .OR. &
       flag_aerosol .EQ. 6 ) THEN 

     CALL readaerosol_interp(id_ASSO4M, rjourvrai, debut, pplay, paprs, t_seri, 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. 6 ) THEN 

     ! Get bc aerosol distribution 
     CALL readaerosol_interp(id_ASBCM, rjourvrai, debut, pplay, paprs, t_seri, bcsol, bcsol_pi )
     CALL readaerosol_interp(id_AIBCM, rjourvrai, debut, pplay, paprs, t_seri, 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. 6 ) THEN

     CALL readaerosol_interp(id_ASPOMM, rjourvrai, debut, pplay, paprs, t_seri, pomsol, pomsol_pi)
     CALL readaerosol_interp(id_AIPOMM, rjourvrai, debut, pplay, paprs, t_seri, pomins, pomins_pi)
  ELSE
     pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
     pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
  END IF


! Read and interpolate csssm, ssssm, assssm
  IF (flag_aerosol .EQ. 4 .OR. &
      flag_aerosol .EQ. 6 ) THEN 

      CALL readaerosol_interp(id_SSSSM ,rjourvrai, debut, pplay, paprs, t_seri, sssupco, sssupco_pi) 
      CALL readaerosol_interp(id_CSSSM ,rjourvrai, debut, pplay, paprs, t_seri, sscoarse,sscoarse_pi) 
      CALL readaerosol_interp(id_ASSSM ,rjourvrai, debut, pplay, paprs, t_seri, ssacu, ssacu_pi) 

  ELSE
     sscoarse(:,:) = 0. ; sscoarse_pi(:,:) = 0. 
     ssacu(:,:)    = 0. ; ssacu_pi(:,:) = 0. 
     sssupco(:,:)  = 0. ; sssupco_pi = 0. 
  ENDIF

! Read and interpolate cidustm
  IF (flag_aerosol .EQ. 5 .OR.  &
      flag_aerosol .EQ. 6 ) THEN 

      CALL readaerosol_interp(id_CIDUSTM, rjourvrai, debut, pplay, paprs, t_seri, cidust, cidust_pi) 

  ELSE
      cidust(:,:) = 0. ; cidust_pi(:,:) = 0. 
  ENDIF

!
! Store all aerosols in one variable
!
  m_allaer(:,:,id_ASBCM)  = bcsol(:,:)        ! ASBCM
  m_allaer(:,:,id_ASPOMM) = pomsol(:,:)       ! ASPOMM
  m_allaer(:,:,id_ASSO4M) = sulfate(:,:)      ! ASSO4M (= SO4) 
  m_allaer(:,:,id_CSSO4M) = 0.                ! CSSO4M 
  m_allaer(:,:,id_SSSSM)  = sssupco(:,:)      ! SSSSM
  m_allaer(:,:,id_CSSSM)  = sscoarse(:,:)     ! CSSSM
  m_allaer(:,:,id_ASSSM)  = ssacu(:,:)        ! ASSSM
  m_allaer(:,:,id_CIDUSTM)= cidust(:,:)       ! CIDUSTM
  m_allaer(:,:,id_AIBCM)  = bcins(:,:)        ! AIBCM
  m_allaer(:,:,id_AIPOMM) = pomins(:,:)       ! AIPOMM

!
! Calculate the total mass of all soluble aersosols
!
  mass_solu_aero(:,:)    = sulfate(:,:)    + bcsol(:,:)    + pomsol(:,:)    + &
       sscoarse(:,:)    + ssacu(:,:)    + sssupco(:,:) 
  mass_solu_aero_pi(:,:) = sulfate_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:) + &
       sscoarse_pi(:,:) + ssacu_pi(:,:) + sssupco_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, presnivs) 
     
     ! aeropt_5wv only for validation and diagnostics.
     CALL aeropt_5wv( &
          pdel, m_allaer, &
          pdtphys, rhcl, aerindex, & 
          flag_aerosol, pplay, t_seri, &
          tausum_aero, tau3d_aero, presnivs)
  ELSE

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

END SUBROUTINE readaerosol_optic