      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 tracer_h, only: noms,rho_co2,rho_ice
       use comcstfi_mod, only: g, pi
       use geometry_mod, only: latitude
       use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl,	&
		CLFvarying,CLFfixval,dusttau,			 	&
		pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,	&
		pres_bottom_strato,pres_top_strato,obs_tau_col_strato,  &
                tau_nh3_cloud, pres_nh3_cloud,                          &
                nlayaero, aeronlay_tauref, aeronlay_choice,             & 
                aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght
                  
       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)
!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
!
!     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
      ! BENJAMIN MODIFS
      real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
      real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
      logical,intent(in) :: clearsky

      real aerosol0, obs_tau_col_aurora, pm
      real pcloud_deck, cloud_slope

      real dp_strato(ngrid)
      real dp_tropo(ngrid)
      real dp_layer(ngrid)

      INTEGER l,ig,iq,iaer,ia

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

      INTEGER,SAVE :: i_co2ice=0      ! co2 ice
      INTEGER,SAVE :: i_h2oice=0      ! water ice
!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
      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
        if (iaero_back2lay.ne.0) then
          print*,'iaero_back2lay= ',iaero_back2lay
        endif
        if (iaero_nh3.ne.0) then
          print*,'iaero_nh3= ',iaero_nh3
        endif
        if (iaero_nlay(1).ne.0) then
          print*,'iaero_nlay= ',iaero_nlay(:)
        endif
        if (iaero_aurora.ne.0) then
          print*,'iaero_aurora= ',iaero_aurora
        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:nlayer,iaer)=0.0
!       2. Opacity calculation
            if (noaero) then ! aerosol set to zero
             aerosol(1:ngrid,1:nlayer,iaer)=0.0
            elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
               aerosol(1:ngrid,1:nlayer,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:nlayer,iaer)=0.0
!       2. Opacity calculation
            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
               aerosol(1:ngrid,1:nlayer,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
                  do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
                                ! same correction should b-probably be done for other aerosol types.
                     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,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
                  do ig=1, ngrid
                     !do l=1,nlayer-1 ! to stop the rad tran bug
                     do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
                        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
                     do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
                        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:nlayer,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=(pplev(ig,1)/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 surface pressure pplev(ig,1)

            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) / pplev(ig,1) &
                       *  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:nlayer,iaer)=0.0


!       2. Opacity calculation

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

               zp=(pplev(ig,1)/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) / pplev(ig,1) &
                       *  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
 
           
!     ---------------------------------------------------------
!==================================================================
!    Two-layer aerosols (unknown composition)
!    S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
!    
!    This scheme is deprecated and left for retrocompatibility
!    You should use the n-layer scheme below !
!
!==================================================================

      if (iaero_back2lay .ne.0) then
           iaer=iaero_back2lay
!       1. Initialization
            aerosol(1:ngrid,1:nlayer,iaer)=0.0
!       2. Opacity calculation


!       JVO 20 : Modif to have each of the layers (strato and tropo) correctly normalized
!                Otherwise we previously had the total optical depth correct but for each
!                separately, so  it didn't match the input values + what's more normalizing
!                to the sum was making them non-independent : eg changing tau_tropo was
!                affecting stratopsheric values of optical depth ...
!
!                Note that the main consequence of the former version bug was (in most cases)
!                to strongly underestimate the stratospheric optical depths compared to the
!                required values, eg, with tau_tropo=10 and tau_strato=0.1, you actually ended
!                with an actual tau_strato of 1E-4 ... !
!
!                NB : Because of the extra transition opacity if the layers are non contiguous,
!                be aware that at the the bottom we have tau > tau_strato + tau_tropo

         DO ig=1,ngrid
          dp_tropo(ig)  = 0.D0
          dp_strato(ig) = 0.D0
          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.D0
             !! 2. tropo layer
             ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
               dp_tropo(ig) = dp_tropo(ig) + aerosol(ig,l,iaer)
             !! 3. linear transition 
             ! JVO 20 : This interpolation needs to be done AFTER we set strato and tropo (see below)
             !! 4. strato layer 
             ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
               dp_strato(ig) = dp_strato(ig) + aerosol(ig,l,iaer)
             !! 5. above strato layer: no aerosols
             ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
               aerosol(ig,l,iaer) = 0.D0
             ENDIF
	  ENDDO
         ENDDO

!       3. Re-normalize to the (input) observed (total) column (for each of the layers)

         DO ig=1,ngrid
          DO l=1,nlayer-1
               IF (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)/dp_tropo(ig)
               ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
                 expfactor=log(pplev(ig,l)/pres_top_tropo)/log(pres_bottom_strato/pres_top_tropo)
                 aerosol(ig,l,iaer) = (obs_tau_col_strato/dp_strato(ig))**expfactor     &
                                    * (obs_tau_col_tropo/dp_tropo(ig))**(1.0-expfactor) &
                                    * aerosol(ig,l,iaer)
               ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
                 aerosol(ig,l,iaer) = obs_tau_col_strato*aerosol(ig,l,iaer)/dp_strato(ig)
               ENDIF
            ENDDO
         ENDDO


      end if ! if Two-layer aerosols  

!==================================================================
!    Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
!    S. Guerlet (2013)
!    JVO 20 : You should now use the generic n-layer scheme below
!==================================================================

      if (iaero_nh3 .ne.0) then
           iaer=iaero_nh3
!       1. Initialization
            aerosol(1:ngrid,1:nlayer,iaer)=0.D0
!       2. Opacity calculation
          DO ig=1,ngrid

           DO l=1,nlayer-1
            !! 1. below cloud layer: no opacity
	    
            IF (pplev(ig,l) .gt. pres_nh3_cloud ) THEN
            aerosol(ig,l,iaer) = 0.D0            

             ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN 
	     cloud_slope=5. !!(hard-coded, correspond to scale height 0.2)
             aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)**(cloud_slope))*tau_nh3_cloud 

             ENDIF
            ENDDO

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

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

     end if ! if NH3 cloud  

!=========================================================================================================
!    Generic N-layers aerosols/clouds
!    Author : J. Vatant d'Ollone (2020)
!    
!    Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud
!    
!    + Each layer can have different optical properties, size of particle ...
!    + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...)
!    + You have different choices for vertical profile of the aerosol layers :
!           * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
!           * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).
!                                   In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.
!           * aeronlay_choice = ... feel free to add more cases  !
!    + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it)
!
!=========================================================================================================

      if (iaero_nlay(1) .ne.0) then

        DO ia=1,nlayaero
           iaer=iaero_nlay(ia)

!          a. Initialization
           aerosol(1:ngrid,1:nlayer,iaer)=0.D0

!          b. Opacity calculation
           
           ! Case 1 : Follows atmospheric scale height between boundaries pressures
           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           IF (aeronlay_choice(ia).eq.1) THEN

             dp_layer(:)=0.D0
             DO ig=1,ngrid
               DO l=1,nlayer-1
                 !! i. Opacity follows scale height
                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
                   aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
                 !! ii. Outside aerosol layer boundaries: no aerosols
                 ELSE
                   aerosol(ig,l,iaer) = 0.D0
                 ENDIF
               ENDDO
             ENDDO
             ! iii. Re-normalize to required total opacity
             DO ig=1,ngrid
               DO l=1,nlayer-1
                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
                                     * aeronlay_tauref(ia)
                 ENDIF
               ENDDO
             ENDDO

           ! Case 2 : Follows input scale height
           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
           ELSE IF (aeronlay_choice(ia).eq.2) THEN
           
             cloud_slope  = 1.D0/aeronlay_sclhght(ia)
             pcloud_deck  = aeronlay_pbot(ia)
             dp_layer(:)  = 0.D0

             DO ig=1,ngrid
               DO l=1,nlayer-1
                 !! i. Below cloud layer: no opacity
                 IF (pplev(ig,l) .gt. pcloud_deck) THEN
                   aerosol(ig,l,iaer) = 0.D0            
                 !! ii. Follows scale height above cloud deck
                 ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN 
                   aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope))
                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
                 ENDIF
               ENDDO
             ENDDO
             ! iii. Re-normalize to required total opacity
             DO ig=1,ngrid
               DO l=1,nlayer-1
                 IF (pplev(ig,l) .le. pcloud_deck) THEN 
                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
                                     * aeronlay_tauref(ia)
                 ENDIF
               ENDDO
             ENDDO

           ENDIF ! aeronlay_choice

          ENDDO ! loop on n aerosol layers

      end if ! if N-layer aerosols
  
!==================================================================
!    Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
!    S. Guerlet (2015)
!==================================================================


      if (iaero_aurora .ne.0) then
           iaer=iaero_aurora
!       1. Initialization
            aerosol(1:ngrid,1:nlayer,iaer)=0.D0 
	 pm = 2000. !!case study: maxi aerosols at 20 hPa
!       2. Opacity calculation
          DO ig=1,ngrid

	  !! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015
              DO l=1,nlayer
       		aerosol(ig,l,iaer) = (pplev(ig,l)/pm)**2 * exp(-(pplev(ig,l)/pm)**2)
              ENDDO
          ENDDO
	 
 !       3. Meridional distribution, and re-normalize to observed total column
         dp_layer(:)=0.D0
         DO ig=1,ngrid
	  !!Jupiter
	  !!Hem sud:
          IF (latitude(ig)*180.D0/pi .lt. -45.D0 .and. latitude(ig)*180.D0/pi .gt. -70.) THEN
 	  obs_tau_col_aurora= 10.D0**(-0.06D0*latitude(ig)*180.D0/pi-3.4D0) 
          ELSEIF (latitude(ig)*180.D0/pi .lt. -37.D0 .and. latitude(ig)*180.D0/pi .ge. -45.) THEN
 	  obs_tau_col_aurora= 10.D0**(-0.3D0*latitude(ig)*180.D0/pi-14.3D0) 
           ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN
 	  obs_tau_col_aurora= 10**(0.06*70.-3.4D0) 
	  !!Hem Nord:  
          ELSEIF (latitude(ig)*180.D0/pi .gt. 30.D0 .and. latitude(ig)*180.D0/pi .lt. 70.) THEN
	  obs_tau_col_aurora= 10.D0**(0.03D0*latitude(ig)*180.D0/pi-1.17D0)  
          ELSEIF (latitude(ig)*180.D0/pi .gt. 22.D0 .and. latitude(ig)*180.D0/pi .le. 30.) THEN
	  obs_tau_col_aurora= 10.D0**(0.3D0*latitude(ig)*180.D0/pi-9.4D0)  
          ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN
	  obs_tau_col_aurora= 10**(0.03*70.-1.17D0)  
          ELSEIF (latitude(ig)*180.D0/pi .ge. -37. .and. latitude(ig)*180.D0/pi .le. 22.) THEN
	 obs_tau_col_aurora = 0.001D0    !!Jupiter: mini pas a zero
	  ENDIF

 	  DO l=1,nlayer  
               dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
          ENDDO
         ENDDO

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


      end if ! if Auroral 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
      
