!==================================================================
module aerosol_mod
implicit none

!==================================================================

!  aerosol indexes: these are initialized to be 0 if the
!                 corresponding aerosol was not activated in callphys.def
!                 -- otherwise a value is set via iniaerosol
      integer, save, protected :: iaero_n2 = 0 
      integer, save, protected :: iaero_dust = 0
      integer, save, protected :: iaero_h2so4 = 0
      logical, save, protected :: noaero = .false.
!$OMP THREADPRIVATE(iaero_n2,iaero_h2o,iaero_dust,iaero_h2so4,noaero)

! two-layer simple aerosol model
      integer, save, protected :: iaero_back2lay = 0
 ! NH3 cloud
      integer, save, protected :: iaero_nh3 = 0
! N-layer aerosol model (replaces the 2-layer and hard-coded clouds)
      integer,dimension(:), allocatable, save, protected :: iaero_nlay
! Auroral aerosols
      integer, save, protected :: iaero_aurora = 0
!$OMP THREADPRIVATE(iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora)

! Generic aerosols
      integer, dimension(:), allocatable, save, protected :: iaero_generic
      integer, dimension(:), allocatable, save, protected :: i_rgcs_ice
!$OMP THREADPRIVATE(iaero_generic,i_rgcs_ice)

!==================================================================

contains

  SUBROUTINE iniaerosol

  use mod_phys_lmdz_para, only : is_master
  use radinc_h, only: naerkind
  use tracer_h, only: n_rgcs, nqtot, is_rgcs
  use callkeys_mod, only: aeron2, dusttau, aeroh2so4, &
                          aeroback2lay, aeronh3, nlayaero, aeronlay, &
                          aeroaurora, aerogeneric, &
                          aerovenus1, aerovenus2, aerovenus2p, &
                          aerovenus3, aerovenusUV

  IMPLICIT NONE
!=======================================================================
!   subject:
!   --------
!   Initialization related to aerosols 
!   (N2 aerosols, dust, water, chemical species, ice...)   
!
!   author: Laura Kerber, S. Guerlet
!   ------
!        
!=======================================================================

  integer :: i, ia, iq

  ! Special case, dyn. allocation for n-layer depending on callphys.def
  IF(.NOT.ALLOCATED(iaero_nlay)) ALLOCATE(iaero_nlay(nlayaero))
  iaero_nlay(:) = 0
  ! Do the same for iaero_generic and i_rgcs_ice
  IF (.not. allocated(iaero_generic)) allocate(iaero_generic(aerogeneric))
  if (.not. allocated(i_rgcs_ice)) allocate(i_rgcs_ice(aerogeneric)) 

  ! Init of i_rgcs_ice
  i_rgcs_ice(:) =0
  ia = 1
  do iq=1,nqtot
    if (is_rgcs(iq) .eq. 1) then 
        i_rgcs_ice(ia)=iq
        ia = ia+1
     endif
  enddo

  iaero_generic(:)=0
  ia=0
  if (aeron2) then
     ia=ia+1
     iaero_n2=ia
  endif
  if (is_master) write(*,*) '--- N2 aerosol = ', iaero_n2

  if (dusttau.gt.0) then
     ia=ia+1
     iaero_dust=ia
  endif
  if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust

  if (aeroh2so4) then
     ia=ia+1
     iaero_h2so4=ia
  endif
  if (is_master) write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
      
  if (aeroback2lay) then
     ia=ia+1
     iaero_back2lay=ia
  endif
  if (is_master) write(*,*) '--- Two-layer aerosol = ', iaero_back2lay

  if (aeronh3) then
     ia=ia+1
     iaero_nh3=ia
  endif
  if (is_master) write(*,*) '--- NH3 Cloud = ', iaero_nh3

  if (aeronlay) then
     do i=1,nlayaero
       ia=ia+1
       iaero_nlay(i)=ia
     enddo
  endif
  if (is_master) write(*,*) '--- N-layer aerosol = ', iaero_nlay

  if (aeroaurora) then
     ia=ia+1
     iaero_aurora=ia
  endif
  if (is_master) write(*,*) '--- Auroral aerosols = ', iaero_aurora

  if (aerogeneric .ne. 0) then 
     do i=1,aerogeneric
        ia = ia+1
        iaero_generic(i) = ia
     enddo
  endif 
      
  if (is_master) then
    write(*,*)'--- Radiative Generic Condensable Species = ',iaero_generic

    write(*,*) '=== Number of aerosols= ', ia
  endif ! of is_master

! For the zero aerosol case, we currently make a dummy n2 aerosol which is zero everywhere.
! (See aeropacity.F90 for how this works). A better solution would be to turn off the 
! aerosol machinery in the no aerosol case, but this would be complicated. LK

  if (ia.eq.0) then  !For the zero aerosol case. 
     ia = 0
     noaero = .true.
  endif

  if (.not.noaero .and. ia.ne.naerkind) then
    if (is_master) then
      print*, 'Aerosols counted not equal to naerkind'
      print*, 'set correct value for nearkind in callphys.def'
      print*, 'which should be ',ia
      print*, 'according to current options in callphys.def'
      print*, 'or change/correct incompatible options there'
      print*, 'Abort in iniaerosol'
    endif
    call abort_physic("iniaerosl",'wrong number of aerosols',1)
  endif ! of if (ia.ne.naerkind)

  END SUBROUTINE iniaerosol

end module aerosol_mod
!==================================================================