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 tracer_h, only: noms
       use comcstfi_mod, only: g
       use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,  &
		pres_bottom_strato,pres_top_strato,obs_tau_col_strato
                  
       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(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

      ! for fixed dust profiles
      real expfactor, zp

      ! identify tracers
      IF (firstcall) THEN

        write(*,*) "Tracers found in aeropacity:"

        if (noaero) then
          print*, "No active aerosols found in aeropacity"
        else
          print*, "If you would like to use aerosols, make sure any old"
          print*, "start files are updated in newstart using the option"
          print*, "q=0"
          write(*,*) "Active aerosols found in aeropacity:"
        endif

        if (iaero_back2lay.ne.0) then
          print*,'iaero_back2lay= ',iaero_back2lay
        endif

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


!     ---------------------------------------------------------   

!==================================================================
!    Two-layer aerosols (unknown composition)
!    S. Guerlet (2013)
!==================================================================

      if (iaero_back2lay .ne.0) then
           iaer=iaero_back2lay
!       1. Initialization
            aerosol(1:ngrid,1:nlayer,iaer)=0.0
!       2. Opacity calculation
          DO ig=1,ngrid
           DO l=1,nlayer-1
             aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
             !! 1. below tropospheric layer: no aerosols
             IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
             !! 2. tropo layer
             ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
               aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)
             !! 3. linear transition
             ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
               expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
               aerosol(ig,l,iaer)= obs_tau_col_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)*aerosol(ig,l,iaer)/1.5
             !! 4. strato layer
             ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
               aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
             !! 5. above strato layer: no aerosols
             ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
             ENDIF
	   ENDDO
          ENDDO

 !       3. Re-normalize to observed total column
         tau_col(:)=0.0
         DO l=1,nlayer
          DO ig=1,ngrid
               tau_col(ig) = tau_col(ig) &
                     + aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
            ENDDO
         ENDDO

         DO ig=1,ngrid
           DO l=1,nlayer-1
                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
           ENDDO
         ENDDO


      end if ! if Two-layer 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