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

       use radinc_h, only : L_TAUMAX,naerkind
       use aerosol_mod
       USE comgeomfi_h
       USE tracer_h, only: noms,rho_co2,rho_ice
                  
       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,nlayermx,naerkind) \ 3d extinction coefficients
!     QREFir3d(ngrid,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 "comvert.h"

      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,nlayermx,naerkind) ! extinction coefficient in the visible
      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayermx,naerkind)
      REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
      ! BENJAMIN MODIFS
      real,intent(in) :: cloudfrac(ngrid,nlayermx) ! cloud fraction
      real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
      logical,intent(in) :: clearsky

      real aerosol0

      INTEGER l,ig,iq,iaer

      LOGICAL,SAVE :: firstcall=.true.

      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

      ! 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

      ! identify tracers
      IF (firstcall) THEN

        write(*,*) "Tracers found in aeropacity:"
        do iq=1,nq
          tracername=noms(iq)
          if (tracername.eq."co2_ice") then
            i_co2ice=iq
          write(*,*) "i_co2ice=",i_co2ice

          endif
          if (tracername.eq."h2o_ice") then
            i_h2oice=iq
            write(*,*) "i_h2oice=",i_h2oice
          endif
        enddo

        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_co2.ne.0).and.(.not.noaero)) then
          print*, 'iaero_co2=  ',iaero_co2
        endif
        if (iaero_h2o.ne.0) then
          print*,'iaero_h2o=  ',iaero_h2o    
        endif
        if (iaero_dust.ne.0) then
          print*,'iaero_dust= ',iaero_dust
        endif
        if (iaero_h2so4.ne.0) then
          print*,'iaero_h2so4= ',iaero_h2so4
        endif

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


!     ---------------------------------------------------------
!==================================================================
!    CO2 ice aerosols
!==================================================================

      if (iaero_co2.ne.0) then
           iaer=iaero_co2
!       1. Initialization
            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
!       2. Opacity calculation
            if (noaero) then ! aerosol set to zero
             aerosol(1:ngrid,1:nlayermx,iaer)=0.0
            elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
               aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
               !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
            else
               DO ig=1, ngrid
                  DO l=1,nlayer-1 ! to stop the rad tran bug

                     aerosol0 =                         &
                          (  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
                     aerosol0           = max(aerosol0,1.e-9)
                     aerosol0           = min(aerosol0,L_TAUMAX)
                     aerosol(ig,l,iaer) = aerosol0
!                     aerosol(ig,l,iaer) = 0.0
!                     print*, aerosol(ig,l,iaer)
!        using cloud fraction
!                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
!                     aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)


                  ENDDO
               ENDDO
            end if ! if fixed or varying
      end if ! if CO2 aerosols   
!==================================================================
!     Water ice / liquid 
!==================================================================

      if (iaero_h2o.ne.0) then 
           iaer=iaero_h2o
!       1. Initialization
            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
!       2. Opacity calculation
            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
               aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9

               ! put cloud at cloudlvl
               if(kastprof.and.(cloudlvl.ne.0.0))then
                  ig=1
                  do l=1,nlayer
                     if(int(cloudlvl).eq.l)then
                     !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
                        print*,'Inserting cloud at level ',l
                        !aerosol(ig,l,iaer)=10.0

                        rho_ice=920.0

                        ! the Kasting approximation
                        aerosol(ig,l,iaer) =                      &
                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
                          !( pq(ig,l,i_h2oice) + 1.E-9 ) *         &
                          ( 4.0e-4 + 1.E-9 ) *         &
                          ( pplev(ig,l) - pplev(ig,l+1) ) / g


                        open(115,file='clouds.out',form='formatted')
                        write(115,*) l,aerosol(ig,l,iaer)
                        close(115)

                        return
                     endif
                  end do

                  call abort
               endif

            else

               do ig=1, ngrid
                  do l=1,nlayer-1 ! to stop the rad tran bug

                     aerosol(ig,l,iaer) =                                    & !modification by BC
                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
                          !( pplev(ig,l) - pplev(ig,l+1) ) / g       !   opacity in the clearsky=true and the 
                                                                     !   clear=false/pq=0 case
                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
                          ( pplev(ig,l) - pplev(ig,l+1) ) / g 

                  enddo
               enddo

               if(CLFvarying)then
                  call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer))
                  do ig=1, ngrid
                     do l=1,nlayer-1 ! to stop the rad tran bug
                        CLFtot  = max(totcloudfrac(ig),0.01)
                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
                     enddo
                  enddo
               else
                  do ig=1, ngrid
                     do l=1,nlayer-1 ! to stop the rad tran bug
                        CLFtot  = CLFfixval
                        aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
                        aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
                     enddo
                  enddo
              end if!(CLFvarying)
            endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
	      
      end if ! End if h2o aerosol

!==================================================================
!             Dust 
!==================================================================
      if (iaero_dust.ne.0) then
          iaer=iaero_dust
!         1. Initialization 
          aerosol(1:ngrid,1:nlayermx,iaer)=0.0
          
          topdust=30.0 ! km  (used to be 10.0 km) LK

!       2. Opacity calculation

!           expfactor=0.
           DO l=1,nlayer-1
             DO ig=1,ngrid
!             Typical mixing ratio profile

                 zp=(preff/pplay(ig,l))**(70./topdust)
                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)

!             Vertical scaling function
              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
               *expfactor


             ENDDO
           ENDDO

!          Rescaling each layer to reproduce the choosen (or assimilated)
!          dust extinction opacity at visible reference wavelength, which
!          is scaled to the "preff" reference surface pressure available in comvert.h
!          and stored in startfi.nc

            taudusttmp(1:ngrid)=0.
              DO l=1,nlayer
                DO ig=1,ngrid
                   taudusttmp(ig) = taudusttmp(ig) &
                          +  aerosol(ig,l,iaer)
                ENDDO
              ENDDO
            DO l=1,nlayer-1
               DO ig=1,ngrid
                  aerosol(ig,l,iaer) = max(1E-20, &
                          dusttau &
                       *  pplev(ig,1) / preff &
                       *  aerosol(ig,l,iaer) &
                       /  taudusttmp(ig))

              ENDDO
            ENDDO
      end if ! If dust aerosol   

!==================================================================
!           H2SO4 
!==================================================================
! added by LK
      if (iaero_h2so4.ne.0) then
         iaer=iaero_h2so4

!       1. Initialization
         aerosol(1:ngrid,1:nlayermx,iaer)=0.0


!       2. Opacity calculation

!           expfactor=0.
         DO l=1,nlayer-1
            DO ig=1,ngrid
!              Typical mixing ratio profile

               zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
               expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)

!             Vertical scaling function
               aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor

            ENDDO
         ENDDO
         tauh2so4tmp(1:ngrid)=0.
         DO l=1,nlayer
            DO ig=1,ngrid
               tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
            ENDDO
         ENDDO
         DO l=1,nlayer-1
            DO ig=1,ngrid
               aerosol(ig,l,iaer) = max(1E-20, &
                          1 &
                       *  pplev(ig,1) / preff &
                       *  aerosol(ig,l,iaer) &
                       /  tauh2so4tmp(ig))

            ENDDO
         ENDDO

! 1/700. is assuming a "sulfurtau" of 1
! Sulfur aerosol routine to be improved.
!                     aerosol0 =                         &
!                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
!                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
!                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
!                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
!                     aerosol0           = max(aerosol0,1.e-9)
!                     aerosol0           = min(aerosol0,L_TAUMAX)
!                     aerosol(ig,l,iaer) = aerosol0

!                  ENDDO
!               ENDDO
      end if
 
           

! --------------------------------------------------------------------------
! 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
      
