Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col)

       use radinc_h, only : L_TAUMAX,naerkind
       use aerosol_mod
       use radii_mod
                  
       implicit none

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"

!==================================================================
!     
!     Purpose
!     -------
!     Compute aerosol optical depth in each gridbox.
!     
!     Authors
!     ------- 
!     F. Forget
!     F. Montmessin (water ice scheme) 
!     update J.-B. Madeleine (2008)
!     Adaptation Pluto : Tanguy Bertrand (2017)
!
!     Input
!     ----- 
!     ngrid             Number of horizontal gridpoints
!     nlayer            Number of layers
!     nq                Number of tracers
!     pplev             Pressure (Pa) at each layer boundary
!     pq                Aerosol mixing ratio
!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
!     QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
!     QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths
!
!     Output
!     ------
!     aerosol            Aerosol optical depth in layer l, grid point ig
!     tau_col            Total column optical depth at grid point ig
!
!=======================================================================

      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
      INTEGER,INTENT(IN) :: nq     ! number of tracers
      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth

      real aerosol0

      INTEGER l,ig,iq,iaer

      LOGICAL,SAVE :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
      REAL CBRT
      EXTERNAL CBRT

      !INTEGER,SAVE :: i_haze=0      ! 
      CHARACTER(LEN=20) :: tracername ! to temporarily store text

      ! for fixed dust profiles
      real topdust, expfactor, zp
      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling

      real CLFtot
!      integer block

      ! identify tracers
      IF (firstcall) THEN

        write(*,*) "Tracers found in aeropacity:"
        do iq=1,nq
          tracername=noms(iq)
          write(*,*) iq,tracername
        enddo

        if (iaero_haze.eq.0) then
          print*, "No active aerosols found in aeropacity"
          print*, "Check traceur.def or haze option in callphys.def"
          stop
        else
          write(*,*) "Active aerosols found in aeropacity:"
        endif

        if (iaero_haze.ne.0) then
          print*, 'iaero_haze=  ',iaero_haze
        endif
        print*, 'naerkind aeropacity=',naerkind
        firstcall=.false.
      ENDIF ! of IF (firstcall)


!     ---------------------------------------------------------
!==================================================================
!    Haze aerosols
!==================================================================

      if (iaero_haze.ne.0) then
           iaer=iaero_haze
!       1. Initialization
            aerosol(1:ngrid,1:nlayer,iaer)=0.0
!       2. Opacity calculation
            DO ig=1, ngrid
                  DO l=1,nlayer-1 ! to stop the rad tran bug
                     ! if fractal, radius doit etre equivalent sphere radius
                     aerosol0 =                         &
                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
                          ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
                          ( pq(ig,l,i_haze) + 1.E-10 ) *         &
                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
                     aerosol0           = max(aerosol0,1.e-10)
                     aerosol0           = min(aerosol0,L_TAUMAX)
                     aerosol(ig,l,iaer) = aerosol0
                  ENDDO
            ENDDO
            !QREF est le meme dans toute la colonne par def si size uniforme
            !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)            
            !print*, 'TB17: rho_q=',rho_q(i_haze)            
            !print*, 'TB17: reffrad=',reffrad(1,:,1)       
            !print*, 'TB17: pq=',pq(1,:,i_haze)
            !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
      end if ! if haze aerosols   
           
! --------------------------------------------------------------------------
! Column integrated visible optical depth in each point (used for diagnostic)

      tau_col(:)=0.0
      do iaer = 1, naerkind
         do l=1,nlayer
            do ig=1,ngrid
               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
            end do
         end do
      end do

      do ig=1,ngrid
         do l=1,nlayer
            do iaer = 1, naerkind
               if(aerosol(ig,l,iaer).gt.1.e3)then
                  print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
                  print*,'at ig=',ig,',  l=',l,', iaer=',iaer
                  print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
                  print*,'reffrad=',reffrad(ig,l,iaer)
               endif
            end do
         end do
      end do

      do ig=1,ngrid
         if(tau_col(ig).gt.1.e3)then
            print*,'WARNING: tau_col=',tau_col(ig)
            print*,'at ig=',ig
            print*,'aerosol=',aerosol(ig,:,:)
            print*,'QREFvis3d=',QREFvis3d(ig,:,:)
            print*,'reffrad=',reffrad(ig,:,:)
         endif
      end do
      return
    end subroutine aeropacity