SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
     &    nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
                                                   
! to use  'getin'
      USE ioipsl_getincom, only: getin
      use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
     &                      igcm_dust_submicron, rho_dust, rho_ice,
     &                      nqdust
      use comgeomfi_h, only: lati, sinlat ! grid point latitudes (rad)
#ifdef DUSTSTORM
      use comgeomfi_h, only: long
      use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
#endif
      use planete_h
      USE comcstfi_h
      use dimradmars_mod, only: naerkind, name_iaer, 
     &            iaerdust,tauvis,
     &            iaer_dust_conrath,iaer_dust_doubleq,
     &            iaer_dust_submicron,iaer_h2o_ice
       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 E. Millour, march 2012:
c         - reference pressure is now set to 610Pa (not 700Pa)
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(ngrid,nlayer,naerkind) \ 3d extinction coefficients
c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
c   omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
c
c   output:
c   -------
c   tauref       Prescribed mean column optical depth at 610 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 "callkeys.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(ngrid,nlayer)
      REAL reffrad(ngrid,nlayer,naerkind)
      REAL nueffrad(ngrid,nlayer,naerkind)
      REAL QREFvis3d(ngrid,nlayer,naerkind)
      REAL QREFir3d(ngrid,nlayer,naerkind)
      REAL omegaREFvis3d(ngrid,nlayer,naerkind)
      REAL omegaREFir3d(ngrid,nlayer,naerkind)
c
c    Local variables :
c    -----------------
      INTEGER l,ig,iq,i,j
      INTEGER iaer           ! Aerosol index
      real topdust(ngrid)
      real zlsconst, zp
      real taueq,tauS,tauN
c     Mean Qext(vis)/Qext(ir) profile
      real msolsir(nlayer,naerkind)
c     Mean Qext(ir)/Qabs(ir) profile
      real mqextsqabs(nlayer,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(ngrid)! Temporary dust opacity
                               !   used before scaling
      REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
      REAL taudustvis(ngrid) ! Dust opacity after scaling
      REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
                               !   "seen" by the GCM.
      REAL taucloudvis(ngrid)! Cloud opacity at visible
                               !   reference wavelength
      REAL taucloudtes(ngrid)! Cloud opacity at infrared
                               !   reference wavelength using
                               !   Qabs instead of Qext
                               !   (direct comparison with TES)
      REAL topdust0(ngrid)

#ifdef DUSTSTORM
!! Local dust storms
      logical localstorm        ! =true to create a local dust storm
      real taulocref,ztoploc,radloc,lonloc,latloc  ! local dust storm parameters
      real reffstorm, yeah
      REAL ray(ngrid) ! distance from dust storm center
      REAL tauuser(ngrid) ! opacity perturbation due to dust storm
      REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
      REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
      real l_top ! layer of the storm's top
      REAL zalt(ngrid, nlayer) ! useful factor to compute l_top 
#endif

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

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 :: cstdustlevel0 = 7
      INTEGER, SAVE      :: cstdustlevel

      LOGICAL,SAVE :: firstcall=.true.

! indexes of water ice and dust tracers:
      INTEGER,SAVE :: i_ice=0  ! water ice
      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
      CHARACTER(LEN=20) :: txt ! to temporarly store text
      CHARACTER(LEN=1) :: txt2 ! to temporarly store text
! indexes of dust scatterers:
      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       typical profile of solsir and (1-w)^(-1):
        msolsir(1:nlayer,1:naerkind)=0
        mqextsqabs(1:nlayer,1:naerkind)=0
        WRITE(*,*) "Typical profiles of Qext(vis)/Qext(IR)"
        WRITE(*,*) "  and Qext(IR)/Qabs(IR):"
        DO iaer = 1, naerkind ! Loop on aerosol kind
          WRITE(*,*) "Aerosol # ",iaer
          DO l=1,nlayer
            DO ig=1,ngrid
              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(ngrid)
            mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
          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)

        IF (freedust) THEN
          cstdustlevel = 1
        ELSE
          cstdustlevel = cstdustlevel0
        ENDIF


#ifndef DUSTSTORM
        firstcall=.false.
#endif

      END IF

c     Vertical column optical depth at "odpref" Pa 
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      IF(freedust) THEN
         tauref(:) = 0. ! tauref is computed after, instead of being forced

      ELSE IF(iaervar.eq.1) THEN 
         do ig=1, ngrid
          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
          if (lati(ig).ge.0) then
          ! Northern hemisphere
            tauref(ig)= tauN +
     &      (taueq-tauN)*0.5*(1+tanh((45-lati(ig)*180./pi)*6/60))
          else
          ! Southern hemisphere
            tauref(ig)= tauS +
     &      (taueq-tauS)*0.5*(1+tanh((45+lati(ig)*180./pi)*6/60))
          endif
        enddo ! of do ig=1,ngrid
      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 IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
      ! clim, cold or warm synthetic scenarios
        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
      ELSE IF ((iaervar.ge.24).and.(iaervar.le.32))
     &     THEN  ! << MY... dust scenarios >>
        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
      ELSE IF ((iaervar.eq.4).or.
     &         ((iaervar.ge.124).and.(iaervar.le.126))) THEN
       ! "old" TES assimation dust scenario (values at 700Pa in files!)
        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
      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
            ! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
            ! in the Viking year scenario
            topdust0(ig)=60. -22.*sinlat(ig)**2
            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.odpref
     $                        /(988.**(topdust(ig)/70.))) then
                zp=(odpref/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=odpref/pplay(ig,l)
                aerosol(ig,l,1)= tauref(ig)/odpref *
     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
                ! DENSITY SCALED OPACITY IN INFRARED:
                dsodust(ig,l) =
     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
     &          pq(ig,cstdustlevel,igcm_dust_mass)
              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
                ! DENSITY SCALED OPACITY IN INFRARED:
                dsodust(ig,l) =
     &          (  0.75 * QREFir3d(ig,l,iaer) /
     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
     &          pq(ig,l,igcm_dust_mass)
              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 -- Now done in physiq.F
!        IF (ngrid.NE.1) THEN
!          CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
!     &      ' ',2,taucloudvis)
!          CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
!     &      ' ',2,taucloudtes)
!          IF (callstats) THEN
!            CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
!     &        ' ',2,taucloudvis)
!            CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
!     &        ' ',2,taucloudtes)
!          ENDIF
!        ELSE
!         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 odpref Pa pressure surface.
c -----------------------------------------------------------------

#ifdef DUSTSTORM
c -----------------------------------------------------------------
! Calculate reference opacity without perturbation
c -----------------------------------------------------------------
      IF (firstcall) THEN
        DO iaer=1,naerdust
          DO l=1,nlayer
            DO ig=1,ngrid
              tauref(ig) = tauref(ig) +
     &                    aerosol(ig,l,iaerdust(iaer))
            ENDDO
          ENDDO
        ENDDO
        tauref(:) = tauref(:) * odpref / pplev(:,1)

c--------------------------------------------------
c Get parameters of the opacity perturbation
c--------------------------------------------------
        iaer=1  ! just change dust

        write(*,*) "Add a local storm ?"
        localstorm=.true. ! default value
        call getin("localstorm",localstorm)
        write(*,*) " localstorm = ",localstorm

        IF (localstorm) THEN
          WRITE(*,*) "********************"
          WRITE(*,*) "ADDING A LOCAL STORM"
          WRITE(*,*) "********************"

          write(*,*) "ref opacity of local dust storm"
              taulocref = 4.25 ! default value
              call getin("taulocref",taulocref)
              write(*,*) " taulocref = ",taulocref

          write(*,*) "target altitude of local storm (km)"
              ztoploc = 10.0 ! default value
              call getin("ztoploc",ztoploc)
              write(*,*) " ztoploc = ",ztoploc

          write(*,*) "radius of dust storm (degree)"
              radloc = 0.5 ! default value
              call getin("radloc",radloc)
              write(*,*) " radloc = ",radloc

          write(*,*) "center longitude of storm (deg)"
              lonloc = 25.0 ! default value
              call getin("lonloc",lonloc)
              write(*,*) " lonloc = ",lonloc

          write(*,*) "center latitude of storm (deg)"
              latloc = -2.5 ! default value
              call getin("latloc",latloc)
              write(*,*) " latloc = ",latloc
        
          write(*,*) "reff storm (mic) 0. for background"
              reffstorm = 0.0 ! default value
              call getin("reffstorm",reffstorm)
              write(*,*) " reffstorm = ",reffstorm

!! LOOP: modify opacity
      DO ig=1,ngrid

      !! distance to the center:
      ray(ig)=SQRT((lati(ig)*180./pi-latloc)**2 +
     &          (long(ig)*180./pi -lonloc)**2)

      !! transition factor for storm
      !! factor is hardcoded -- increase it to steepen
      yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.

      !! new opacity field
      !! -- add an opacity set to taulocref
      !! -- the additional reference opacity will
      !!      thus be taulocref*odpref/pplev
      tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref , 
     &          taulocref * yeah )

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

     !! change reffrad if ever needed
      IF (reffstorm .gt. 0.) THEN
          DO l=1,nlayer
             IF (l .lt. l_top+1) THEN
                reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
     &          * 1.e-6 * yeah )
             ENDIF
          ENDDO
      ENDIF

      ENDDO
!! END LOOP

      !! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
      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
          !! Mass mixing ratio perturbation due to local dust storm in each layer
          more_dust(ig,l,1)=
     &                     (tauuser(ig)-(tauref(ig)
     &                      * pplev(ig,1) /odpref)) /
     &                      int_factor(ig)
          more_dust(ig,l,2)=
     &                     (tauuser(ig)-(tauref(ig) *
     &                     pplev(ig,1) /odpref))
     &                      / int_factor(ig) *
     &                     ((ref_r0/reffrad(ig,l,iaer))**3)
     &                      * r3n_q 
          ENDDO
      ENDDO

      !! quantity of dust for each layer with the addition of the perturbation
      DO l=1, l_top
          pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
     .            + more_dust(:,l,1)
          pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
     .            + more_dust(:,l,2)
      ENDDO
      ENDIF !! IF (localstorm)
      tauref(:)=0.
      ENDIF !! IF (firstcall)
#endif

      IF (freedust) THEN
          tauscaling(:) = 1.

      ELSE
c       Temporary scaling factor
        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

c       Saved scaling factor
        DO ig=1,ngrid
            tauscaling(ig) = tauref(ig) * 
     &                       pplev(ig,1) / odpref / taudusttmp(ig)
        ENDDO

      ENDIF ! IF (freedust)

c     Opacity computation
      DO iaer=1,naerdust
        DO l=1,nlayer
          DO ig=1,ngrid
            aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
     &                aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
          ENDDO
        ENDDO
      ENDDO

      DO l=1,nlayer
        DO ig=1,ngrid
          dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
        ENDDO
      ENDDO

      IF (freedust) THEN
        ! tauref has been initialized to 0 before.
        DO iaer=1,naerdust
          DO l=1,nlayer
            DO ig=1,ngrid
#ifdef DUSTSTORM
      !! recalculate opacity because storm perturbation has been added
      IF (firstcall) THEN
              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
      ENDIF
#endif
              tauref(ig) = tauref(ig) + 
     &                    aerosol(ig,l,iaerdust(iaer))
            ENDDO
          ENDDO
        ENDDO
        tauref(:) = tauref(:) * odpref / pplev(:,1)
      ENDIF
      
c output for debug
c        IF (ngrid.NE.1) THEN
c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
c     &      '#',2,taudusttmp)
c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
c     &      '#',2,tauscaling)
c        ELSE
c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
c     &      '#',0,taudusttmp)
c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
c     &      '#',0,tauscaling)
c        ENDIF
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

#ifdef DUSTSTORM
      IF (firstcall) THEN
        firstcall=.false.
      ENDIF
#endif

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

c      IF (ngrid.NE.1) THEN
c       CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
c    &                    'm2.kg-1',3,dsodust)
c        IF (callstats) THEN
c          CALL wstats(ngrid,'dsodust',
c     &               'tau*g/dp',
c     &               'm2.kg-1',3,dsodust)
c        ENDIF
c      ELSE
c        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
c     &                   dsodust)
c      ENDIF ! of IF (ngrid.NE.1)
c -----------------------------------------------------------------
      return
      end