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

       use radinc_h, only : naerkind
                  
       implicit none

!==================================================================
!     
!     Purpose
!     -------
!     Compute aerosol optical depth in each gridbox.
!     
!     Authors
!     ------- 
!     F. Forget
!     F. Montmessin (water ice scheme) 
!     update J.-B. Madeleine (2008)
!     dust removal, simplification by Robin Wordsworth (2009)
!
!     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(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
!     QREFir3d(ngridmx,nlayermx,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
!
!=======================================================================

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

      INTEGER ngrid,nlayer,nq

      REAL pplev(ngrid,nlayer+1)
      REAL pq(ngrid,nlayer,nq)
      REAL aerosol(ngrid,nlayer,naerkind)
      REAL reffrad(ngrid,nlayer,naerkind)
      REAL QREFvis3d(ngridmx,nlayermx,naerkind)
      REAL QREFir3d(ngridmx,nlayermx,naerkind)

      REAL tauref(ngrid), tau_col(ngrid)

      INTEGER l,ig,iq,iaer

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall

      REAL CBRT
      EXTERNAL CBRT

      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
      INTEGER,SAVE :: i_h2oice=0      ! water ice
      CHARACTER(LEN=20) :: tracername ! to temporarily store text


! identify tracers
      IF (firstcall) THEN

         ! are these tests of any real use?
        IF(ngrid.NE.ngridmx) THEN
            PRINT*,'STOP in aeropacity!'
            PRINT*,'problem of dimensions:'
            PRINT*,'ngrid  =',ngrid
            PRINT*,'ngridmx  =',ngridmx
            STOP
        ENDIF

        if (nq.gt.nqmx) then
           write(*,*) 'STOP in aeropacity: (nq .gt. nqmx)!'
           write(*,*) 'nq=',nq,' nqmx=',nqmx
           stop
        endif


        do iq=1,nqmx
          tracername=noms(iq)
          if (tracername.eq."co2_ice") then
            i_co2ice=iq
          endif
          if (tracername.eq."h2o_ice") then
            i_h2oice=iq
          endif
        enddo
        
        write(*,*) "aeropacity: i_co2ice=",i_co2ice
        write(*,*) "            i_h2oice=",i_h2oice
        
        !if((i_h2oice.lt.1) .and. (naerkind.gt.1))then
        !   print*,'naerkind > 1 but no h2o ice, this will not work.'
        !   stop
        !endif

        firstcall=.false.
      ENDIF ! of IF (firstcall)


      aerosol(:,:,:)=0.0

      DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer)
!     ---------------------------------------------------------
        aerkind: SELECT CASE (iaer)
!==================================================================
        CASE(1) aerkind                          ! CO2 ice (iaer=1)
!==================================================================

!       1. Initialization
        CALL zerophys(ngrid*nlayer,aerosol(1,1,iaer))

!       2. Opacity calculation
        if (aerofixed.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
           do ig=1, ngrid
              do l=1,nlayer
                 aerosol(ig,l,iaer) =1.e-9
              end do
              aerosol(ig,14,iaer)=1.e-9 ! single cloud layer option
           end do
        else
           DO ig=1, ngrid
              DO l=1,nlayer
                 aerosol(ig,l,iaer) =                   &
                      (  0.75 * QREFvis3d(ig,l,iaer) /  &
                      ( rho_co2 * reffrad(ig,l,iaer) )  ) * &
                      ( pq(ig,l,i_co2ice) + 1.E-9 ) *   &
                      ( pplev(ig,l) - pplev(ig,l+1) ) / g
              ENDDO
           ENDDO
        end if

!==================================================================
        CASE(2) aerkind               ! Water ice crystals (iaer=2)
!==================================================================


!       1. Initialization
        CALL zerophys(ngrid*nlayer,aerosol(1,1,iaer))

!       2. Opacity calculation
!     if (1.eq.1) then
     if (aerofixed.or.(i_h2oice.eq.0)) then
!        print*,'No H2O clouds in the radiative transfer!'
        do ig=1, ngrid
           do l=1,nlayer
              aerosol(ig,l,iaer) =1.e-9
           end do
           aerosol(ig,5,iaer)=1.e-9 ! single cloud layer option
        end do
     else
        DO ig=1, ngrid
          DO l=1,nlayer
             aerosol(ig,l,iaer) =                   &
            (  0.75 * QREFvis3d(ig,l,iaer) /        &
              ( rho_ice * reffrad(ig,l,iaer) )  ) * &
              ( pq(ig,l,i_h2oice) + 1.E-9 ) *   &
              ( pplev(ig,l) - pplev(ig,l+1) ) / g
          ENDDO
        ENDDO

     end if


!==================================================================
        END SELECT aerkind

      ENDDO ! iaer (loop on aerosol kind)


! --------------------------------------------------------------------------
! Column integrated visible optical depth in each point (used for diagnostic)

      call zerophys(ngrid,tau_col)
      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

      return
    end subroutine aeropacity