      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
     &    pq,ccn,tauref,tau,aerosol,reffrad,nueffrad,
     &    QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,zdqnorm,dsodust)
                                                   
! to use  'getin'
      USE ioipsl_getincom 
       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Computing aerosol optical depth in each gridbox.
c
c   author: F.Forget 
c   ------
c   update F. Montmessin (water ice scheme) 
c      and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
c   update J.-B. Madeleine 2008-2009:
c       - added 3D scattering by aerosols;
c       - dustopacity transferred from physiq.F to callradite.F,
c           and renamed into aeropacity.F;
c   update J FAURE & A SPIGA
c      -Interactive dust for regional dust storms
c
c
c   
c   input:
c   ----- 
c   ngrid             Number of gridpoint of horizontal grid
c   nlayer            Number of layer
c   nq                Number of tracer
c   zday                  Date (time since Ls=0, in martian days)
c   ls                Solar longitude (Ls) , radian
c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
c   pq                Dust mixing ratio (used if tracer =T and active=T).
c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
c   QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
c   QREFir3d(ngridmx,nlayermx,naerkind)  / at reference wavelengths;
c   omegaREFvis3d(ngridmx,nlayermx,naerkind) \ 3d single scat. albedo
c   omegaREFir3d(ngridmx,nlayermx,naerkind)  / at reference wavelengths;
c
c   output:
c   -------
c   tauref       Prescribed mean column optical depth at 700 Pa 
c   tau          Column total visible dust optical depth at each point
c   aerosol      aerosol(ig,l,1) is the dust optical
c                depth in layer l, grid point ig

c
c=======================================================================
#include "dimensions.h"
#include "dimphys.h"
#include "callkeys.h"
#include "comcstfi.h"
#include "comgeomfi.h"
#include "dimradmars.h"
#include "yomaer.h"
#include "tracer.h"
#include "planete.h"
#include "aerkind.h"

c-----------------------------------------------------------------------
c
c    Declarations :
c    --------------
c
c    Input/Output
c    ------------
      INTEGER ngrid,nlayer,nq

      REAL ls,zday,expfactor    
      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
      REAL pq(ngrid,nlayer,nq)
      REAL tauref(ngrid), tau(ngrid,naerkind)
      REAL aerosol(ngrid,nlayer,naerkind)
      REAL dsodust(ngridmx,nlayermx)                       !optical depth in each layer. Diagnostic.
      REAL reffrad(ngrid,nlayer,naerkind)
      REAL nueffrad(ngrid,nlayer,naerkind)
      REAL QREFvis3d(ngridmx,nlayermx,naerkind)
      REAL QREFir3d(ngridmx,nlayermx,naerkind)
      REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
      REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
      REAL zdqnorm(ngridmx,nlayermx,2)                       !mass mixing ratio perturbation due to the dust storm. Output for meso_physiq.F



c
c    Local variables :
c    -----------------
      INTEGER l,ig,iq,i,j
      INTEGER iaer           ! Aerosol index
      real topdust(ngridmx)
      real zlsconst, zp
      real taueq,tauS,tauN
c     Mean Qext(vis)/Qext(ir) profile
      real msolsir(nlayermx,naerkind)
c     Mean Qext(ir)/Qabs(ir) profile
      real mqextsqabs(nlayermx,naerkind)
c     Variables used when multiple particle sizes are used
c       for dust or water ice particles in the radiative transfer
c       (see callradite.F for more information).
      REAL taudusttmp(ngridmx)! Temporary dust opacity
                               !   used before scaling
      REAL taudustvis(ngridmx) ! Dust opacity after scaling
      REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
                               !   "seen" by the GCM.
      REAL taucloudvis(ngridmx)! Cloud opacity at visible
                               !   reference wavelength
      REAL taucloudtes(ngridmx)! Cloud opacity at infrared
                               !   reference wavelength using
                               !   Qabs instead of Qext
                               !   (direct comparison with TES)
      REAL qdust(ngridmx,nlayermx) ! True dust mass mixing ratio
      REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
                                   !   (particules kg-1)
      REAL qtot(ngridmx)           ! Dust column (kg m-2)

c     CCN reduction factor
      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant

c
c   local saved variables
c   ---------------------

      REAL topdust0(ngridmx) 
      SAVE topdust0
c     Level under which the dust mixing ratio is held constant
c       when computing the dust opacity in each layer
c       (this applies when doubleq and active are true)
      INTEGER, PARAMETER :: cstdustlevel = 7

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall
! Local dust storms

         logical justbackground    !to switch on/off dust absorption
         logical localstorm        ! =true to create a local dust storm
         real taulocref,ztoploc,radloc,lonloc,latloc  !local dust storm parameters
         REAL ray(ngridmx)                            !distance from dust storm center
         REAL tauuser(ngridmx)                        !opacity perturbation due to dust storm
         REAL more_dust(ngridmx,nlayermx,2)           !Mass mixing ratio perturbation due to the dust storm
         REAL int_factor(ngridmx)                     !usefull factor to compute mmr perturbation
         real l_top                                   !layer of the dust storm's top
         REAL zalt(ngridmx, nlayer)                   !usefull factor to compute l_top



! indexes of water ice and dust tracers:
      INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers
      INTEGER,SAVE :: i_ice=0  ! water ice
      CHARACTER(LEN=20) :: txt ! to temporarly store text
      CHARACTER(LEN=1) :: txt2 ! to temporarly store text
! indexes of dust scatterers:
      INTEGER,SAVE :: iaerdust(naerkind)
      INTEGER,SAVE :: naerdust ! number of dust scatterers

      tau(1:ngrid,1:naerkind)=0

! identify tracers

      IF (firstcall) THEN
        ! identify scatterers that are dust
        naerdust=0
        DO iaer=1,naerkind
          txt=name_iaer(iaer)
          IF (txt(1:4).eq."dust") THEN
            naerdust=naerdust+1
            iaerdust(naerdust)=iaer
          ENDIF
        ENDDO
        ! identify tracers which are dust
        i=0
        DO iq=1,nq
          txt=noms(iq)
          IF (txt(1:4).eq."dust") THEN
          i=i+1
          nqdust(i)=iq
          ENDIF
        ENDDO

        IF (water.AND.activice) THEN
          i_ice=igcm_h2o_ice
          write(*,*) "aeropacity: i_ice=",i_ice
        ENDIF

c       altitude of the top of the aerosol layer (km) at Ls=2.76rad:
c       in the Viking year scenario
        DO ig=1,ngrid
            topdust0(ig)=60. -22.*SIN(lati(ig))**2
        END DO

c       typical profile of solsir and (1-w)^(-1):
        msolsir(1:nlayer,1:naerkind)=0
        mqextsqabs(1:nlayer,1:naerkind)=0
        WRITE(*,*) "Typical profiles of solsir and Qext/Qabs(IR):"
        DO iaer = 1, naerkind ! Loop on aerosol kind
          WRITE(*,*) "Aerosol # ",iaer
          DO l=1,nlayer
            DO ig=1,ngridmx
              msolsir(l,iaer)=msolsir(l,iaer)+
     &              QREFvis3d(ig,l,iaer)/
     &              QREFir3d(ig,l,iaer)
              mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
     &              (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
            ENDDO
            msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
            mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
          ENDDO
          WRITE(*,*) "solsir: ",msolsir(:,iaer)
          WRITE(*,*) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
        ENDDO

!       load value of tauvis from callphys.def (if given there,
!       otherwise default value read from starfi.nc file will be used)
        call getin("tauvis",tauvis)

c        firstcall=.false.

      END IF

c     Vertical column optical depth at 700.Pa 
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF(iaervar.eq.1) THEN 
         do ig=1, ngridmx
          tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
                                       ! or read in starfi
        end do
      ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>

        tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
        do ig=2,ngrid
          tauref(ig) = tauref(1)
        end do

      ELSE IF (iaervar.eq.3) THEN  ! << "MGS" scenario >>

        taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14
        tauS= 0.1 +(0.5-0.1)  *(cos(0.5*(ls-4.363)))**14
        tauN = 0.1
c	   if (peri_day.eq.150) then
c	     tauS=0.1
c	     tauN=0.1 +(0.5-0.1)  *(cos(0.5*(ls+pi-4.363)))**14
c	     taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14
c           endif
        do ig=1,ngrid/2  ! Northern hemisphere
          tauref(ig)= tauN +
     &    (taueq-tauN)*0.5*(1+tanh((45-lati(ig)*180./pi)*6/60))
        end do
        do ig=ngrid/2+1, ngridmx  ! Southern hemisphere
          tauref(ig)= tauS +
     &    (taueq-tauS)*0.5*(1+tanh((45+lati(ig)*180./pi)*6/60))
        end do
      ELSE IF ((iaervar.eq.4).or.
     &        ((iaervar.ge.24).and.(iaervar.le.26)))
     &     THEN  ! << "TES assimilated dust scenarios >>
        call readtesassim(ngrid,nlayer,zday,pplev,tauref)

      ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
c         tauref(1) = 0.2
c         if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.))
c    &                              tauref(1) = 2.5
        tauref(1) = 2.5
        if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
     &                              tauref(1) = .2

        do ig=2,ngrid
          tauref(ig) = tauref(1)
        end do
      ELSE
        stop 'problem with iaervar in aeropacity.F'
      ENDIF



c -----------------------------------------------------------------
c Computing the opacity in each layer
c -----------------------------------------------------------------

      DO iaer = 1, naerkind ! Loop on aerosol kind
c     --------------------------------------------
        aerkind: SELECT CASE (name_iaer(iaer))
c==================================================================
        CASE("dust_conrath") aerkind      ! Typical dust profile
c==================================================================

c       Altitude of the top of the dust layer
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        zlsconst=SIN(ls-2.76)
        if (iddist.eq.1) then
          do ig=1,ngrid
             topdust(ig)=topdustref         ! constant dust layer top
          end do

        else if (iddist.eq.2) then          ! "Viking" scenario
          do ig=1,ngrid
            topdust(ig)=topdust0(ig)+18.*zlsconst
          end do

        else if(iddist.eq.3) then         !"MGS" scenario
          do ig=1,ngrid
            topdust(ig)=60.+18.*zlsconst
     &                -(32+18*zlsconst)*sin(lati(ig))**4
     &                 - 8*zlsconst*(sin(lati(ig)))**5
          end do
        endif

c       Optical depth in each layer :
c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if(iddist.ge.1) then

          expfactor=0.
          DO l=1,nlayer
            DO ig=1,ngrid
c             Typical mixing ratio profile 
              if(pplay(ig,l).gt.700.
     $                        /(988.**(topdust(ig)/70.))) then
                zp=(700./pplay(ig,l))**(70./topdust(ig))
                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
              else    
                expfactor=1.e-3
              endif
c             Vertical scaling function
              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
     &          expfactor *
     &          QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
            ENDDO
          ENDDO

        else if(iddist.eq.0) then   
c         old dust vertical distribution function (pollack90)
          DO l=1,nlayer
             DO ig=1,ngrid
                zp=700./pplay(ig,l)
                aerosol(ig,l,1)= tauref(ig)/700. *
     s           (pplev(ig,l)-pplev(ig,l+1))
     s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
             ENDDO
          ENDDO
        end if

c==================================================================
        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
c        (transport of mass and number mixing ratio)
c==================================================================
             
          DO l=1,nlayer
            IF (l.LE.cstdustlevel) THEN
c           Opacity in the first levels is held constant to 
c             avoid unrealistic values due to constant lifting:
              DO ig=1,ngrid
                aerosol(ig,l,iaer) = 
     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
     &          pq(ig,cstdustlevel,igcm_dust_mass) *
     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
              ENDDO
            ELSE
              DO ig=1,ngrid
                aerosol(ig,l,iaer) =
     &          (  0.75 * QREFvis3d(ig,l,iaer) /
     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
     &          pq(ig,l,igcm_dust_mass) *
     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
              ENDDO
            ENDIF
          ENDDO

c==================================================================
        CASE("dust_submicron") aerkind   ! Small dust population
c==================================================================

          DO l=1,nlayer
            IF (l.LE.cstdustlevel) THEN
c           Opacity in the first levels is held constant to 
c             avoid unrealistic values due to constant lifting:
              DO ig=1,ngrid
                aerosol(ig,l,iaer) = 
     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
     &          pq(ig,cstdustlevel,igcm_dust_submicron) *
     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
              ENDDO
            ELSE
              DO ig=1,ngrid
                aerosol(ig,l,iaer) = 
     &          (  0.75 * QREFvis3d(ig,l,iaer) /
     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
     &          pq(ig,l,igcm_dust_submicron) *
     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
              ENDDO
            ENDIF
          ENDDO

c==================================================================
        CASE("h2o_ice") aerkind             ! Water ice crystals
c==================================================================

c       1. Initialization
        aerosol(1:ngrid,1:nlayer,iaer) = 0.
        taucloudvis(1:ngrid) = 0.
        taucloudtes(1:ngrid) = 0.
c       2. Opacity calculation
        DO ig=1, ngrid
          DO l=1,nlayer
            aerosol(ig,l,iaer) = max(1E-20,
     &        (  0.75 * QREFvis3d(ig,l,iaer) /
     &        ( rho_ice * reffrad(ig,l,iaer) )  ) *
     &        pq(ig,l,i_ice) *
     &        ( pplev(ig,l) - pplev(ig,l+1) ) / g
     &                              )
            taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
            taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
     &        QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
     &        ( 1.E0 - omegaREFir3d(ig,l,iaer) )
          ENDDO
        ENDDO
c       3. Outputs
        IF (ngrid.NE.1) THEN
          CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
     &      ' ',2,taucloudvis)
          CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
     &      ' ',2,taucloudtes)
          IF (callstats) THEN
            CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
     &        ' ',2,taucloudvis)
            CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
     &        ' ',2,taucloudtes)
          ENDIF
        ELSE
c         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
        ENDIF
c==================================================================
        END SELECT aerkind
c     -----------------------------------
      ENDDO ! iaer (loop on aerosol kind)
c -----------------------------------------------------------------
c Rescaling each layer to reproduce the choosen (or assimilated)
c   dust extinction opacity at visible reference wavelength, which
c   is originally scaled to an equivalent 700Pa pressure surface.
c -----------------------------------------------------------------

c-----------------------------------------------------------------
c Chose justbakground=false to create an interactive local dust storm 
c Switch justbackground to false to enable storm dust absoption
c 
c Note that if justbackground=false, dust background has a conrath
c repartition and absorb ligth.
c--------------------------------------------------------------------

      justbackground=.false.

      IF (justbackground .eq. .true.)  THEN

      taudusttmp(1:ngrid)=0.
      DO iaer=1,naerdust
        DO l=1,nlayer
          DO ig=1,ngrid
c           Scaling factor
            taudusttmp(ig) = taudusttmp(ig) +
     &                       aerosol(ig,l,iaerdust(iaer))
          ENDDO
        ENDDO
      ENDDO
      DO iaer=1,naerdust
        DO l=1,nlayer
          DO ig=1,ngrid
            aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
     &                   tauref(ig) *
     &                   pplev(ig,1) / 700.E0 *
     &                   aerosol(ig,l,iaerdust(iaer)) /
     &                   taudusttmp(ig)
     &                                        )
          ENDDO
        ENDDO
      ENDDO

      ENDIF

c -----------------------------------------------------------------
c the quantity of dust to add at the first time step is calculated to match
c a tunable opacity perturbation.
c -----------------------------------------------------------------

      IF (firstcall) THEN
      WRITE(*,*) " RENORMALISATION !!! "
c--------------------------------------------------
c  Parameters of the opacity perturbation
c--------------------------------------------------

      iaer=1  !!!! PROVISOIRE !!!!
      taulocref = 4.25 !10  ! ref optical depth of the local dust storm
      ztoploc = 10      ! target pseudo-altitude of local storm (km)
      radloc = 0.5      ! radius of dust storm (degree)
      lonloc = 85       ! center longitude of storm (deg)
      latloc = 25       ! center latitude of storm (deg)

      DO ig=1,ngrid

c---------------------------------------
c        distance to the center:
c-----------------------------------------

      ray(ig)=SQRT((lati(ig)*180./pi-latloc)**2 +
     &          (long(ig)*180./pi -lonloc)**2)

c-------------------------------------------------
c           Tau's new map:
c------------------------------------------

      tauuser(ig)=max(tauref(ig) * pplev(ig,1) / 700.E0 , taulocref
     &          * (TANH(2.+(radloc-ray(ig))*2)+1.)/2.)

c---------------------------------------------------------
c           compute l_top
c----------------------------------------------------------

          DO l=1,nlayer

            zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
     &                      / g / 44.01
     &                    * 8.31 * 210.
                IF (     (ztoploc .lt. zalt(ig,l)  )
     &          .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1

          ENDDO
      ENDDO

c----------------------------------------------------------------------------
c    compute int_factor
c---------------------------------------------------------------------------


      DO ig=1,ngrid
          int_factor(ig)=0.
          DO l=1,nlayer
             IF (l .lt. l_top+1) THEN
                      int_factor(ig) =
     &                int_factor(ig) +
     &          (  0.75 * QREFvis3d(ig,l,iaer) /
     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
             ENDIF
          ENDDO
          DO l=1, nlayer

c-------------------------------------------------------------------------
c Mass mixing ratio perturbation due to the local dust storm in each layer
c-------------------------------------------------------------------------
          more_dust(ig,l,1)=
     &                     (tauuser(ig)-(tauref(ig)
     &                      * pplev(ig,1) / 700.E0)) /
     &                      int_factor(ig)
          more_dust(ig,l,2)=
     &                     (tauuser(ig)-(tauref(ig) *
     &                     pplev(ig,1) / 700.E0))
     &                      / int_factor(ig) *
     &                     ((ref_r0/reffrad(ig,l,iaer))**3)
     &                      * r3n_q 
          ENDDO
      ENDDO

c--------------------------------------------------------------------------------------
c   quantity of dust whiwh is added at the first time step in dynamical core.
c--------------------------------------------------------------------------------------
      DO l=1, l_top
      zdqnorm(:,l,1) = more_dust(:,l,1)
      zdqnorm(:,l,2) = more_dust(:,l,2)
      ENDDO

!        DO ig=1, ngrid
!            DO l=1, l_top
!                aerosol(ig,l,iaer) =
!     &          (  0.75 * QREFvis3d(ig,l,iaer) /
!     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
!     &          more_dust(ig,l,1) *
!     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
!             ENDDO
!        ENDDO


      firstcall = .false.
      ENDIF  !!! firstcall

c -----------------------------------------------------------------
c Computing the number of condensation nuclei
c -----------------------------------------------------------------
      DO iaer = 1, naerkind ! Loop on aerosol kind
c     --------------------------------------------
        aerkind2: SELECT CASE (name_iaer(iaer))
c==================================================================
        CASE("dust_conrath") aerkind2     ! Typical dust profile
c==================================================================
          DO l=1,nlayer
            DO ig=1,ngrid
              ccn(ig,l) = max(aerosol(ig,l,iaer) /
     &                  pi / QREFvis3d(ig,l,iaer) *
     &                  (1.+nueffrad(ig,l,iaer))**3. /
     &                  reffrad(ig,l,iaer)**2. * g /
     &                  (pplev(ig,l)-pplev(ig,l+1)),1e-30)
            ENDDO
          ENDDO
c==================================================================
        CASE("dust_doubleq") aerkind2!Two-moment scheme for dust
c        (transport of mass and number mixing ratio)
c==================================================================
          qtot(1:ngridmx) = 0.
          DO l=1,nlayer
            DO ig=1,ngrid

c--------------------------------------------------------------
c commented useless because (mass mixing ration have a real physical sense now)
c---------------------------------------------------------------------
c       if (localstorm. NE. .true.) then
c              qdust(ig,l) = pq(ig,l,igcm_dust_mass) * tauref(ig) *
c     &                      pplev(ig,1) / 700.E0 / taudusttmp(ig)
c       else
c              qdust(ig,l) = pq(ig,l,igcm_dust_mass) * taureftache(ig) *
c     &                      pplev(ig,1) / 700.E0 / taudusttmp(ig)
c       endif
c              qtot(ig) = qtot(ig) + qdust(ig,l) *
c     &                   (pplev(ig,l)-pplev(ig,l+1)) / g
              ccn(ig,l) = max( ( ref_r0 /
     &                    reffrad(ig,l,iaer) )**3. *
     &                    r3n_q * pq(ig,l,igcm_dust_mass) ,1e-30)
            ENDDO
          ENDDO
c==================================================================
        END SELECT aerkind2
c     -----------------------------------
      ENDDO ! iaer (loop on aerosol kind)


c -----------------------------------------------------------------
c -----------------------------------------------------------------
c  Reduce number of nuclei
!         TEMPORAIRE : r�duction du nombre de nuclei FF 04/200
!         reduction facteur 3
!         ccn(ig,l) = ccn(ig,l) / 27.
!         reduction facteur 2
!         ccn(ig,l) = ccn(ig,l) / 8.
c -----------------------------------------------------------------
       write(*,*) "water_param CCN reduc. fac. ", ccn_factor
       DO l=1,nlayer
         DO ig=1,ngrid
            ccn(ig,l) = ccn(ig,l) / ccn_factor
         ENDDO
       ENDDO
c -----------------------------------------------------------------
c -----------------------------------------------------------------


c -----------------------------------------------------------------
c Column integrated visible optical depth in each point
c -----------------------------------------------------------------
      DO iaer=1,naerkind
        do l=1,nlayer
           do ig=1,ngrid
             tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
           end do
        end do
      ENDDO

                
c -----------------------------------------------------------------
c Density scaled opacity and column opacity output
c -----------------------------------------------------------------
      dsodust(1:ngrid,1:nlayer) = 0.
      DO iaer=1,naerdust
        DO l=1,nlayer
          DO ig=1,ngrid
            dsodust(ig,l) = dsodust(ig,l) +
     &                      aerosol(ig,l,iaerdust(iaer)) * g /
     &                      (pplev(ig,l) - pplev(ig,l+1))
          ENDDO
        ENDDO
        IF (ngrid.NE.1) THEN
          write(txt2,'(i1.1)') iaer
          call WRITEDIAGFI(ngridmx,'taudust'//txt2,
     &                    'Dust col opacity',
     &                    ' ',2,tau(1,iaerdust(iaer)))
          IF (callstats) THEN
            CALL wstats(ngridmx,'taudust'//txt2,
     &                 'Dust col opacity',
     &                 ' ',2,tau(1,iaerdust(iaer)))
          ENDIF
        ENDIF
      ENDDO

      IF (ngrid.NE.1) THEN
c       CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp',
c    &                    'm2.kg-1',3,dsodust)
        IF (callstats) THEN
          CALL wstats(ngridmx,'dsodust',
     &               'tau*g/dp',
     &               'm2.kg-1',3,dsodust)
        ENDIF
c       CALL WRITEDIAGFI(ngridmx,'ccn','Cond. nuclei',
c    &                    'part kg-1',3,ccn)
      ELSE
        CALL writeg1d(ngrid,nlayer,dsodust,'dsodust','m2.kg-1')
      ENDIF
c -----------------------------------------------------------------
      return
      end 
