SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,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
      use comgeomfi_h, only: lati ! grid point latitudes (rad)
      use yomaer_h, only: tauvis
       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 "dimensions.h"
!#include "dimphys.h"
#include "callkeys.h"
#include "comcstfi.h"
!#include "comgeomfi.h"
!#include "dimradmars.h"
!#include "yomaer.h"
!#include "tracer.h"
! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
#include"scatterers.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(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)

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

      REAL,SAVE,ALLOCATABLE :: 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,SAVE :: firstcall=.true.

! indexes of water ice and dust tracers:
      INTEGER,ALLOCATABLE,SAVE :: nqdust(:) ! to store the indexes of 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 :: iaerdust(naerkind)
      INTEGER,SAVE :: naerdust ! number of dust scatterers

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

! identify tracers

      IF (firstcall) THEN
        ! allocate local saved arrays
        allocate(nqdust(nq))
        allocate(topdust0(ngrid))
        
        ! 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 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)

        firstcall=.false.

      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.7)) THEN
      ! cold or warm synthetic scenarios
        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
      ELSE IF ((iaervar.ge.24).and.(iaervar.le.30))
     &     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
            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
              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 -- 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 -----------------------------------------------------------------

      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)
c            tauscaling(ig) = 1.e-4
        ENDDO

      ENDIF

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

      IF (freedust) THEN
       ! tauref has been initialized to 0 before.
       DO iaer=1,naerdust
        DO l=1,nlayer
         DO ig=1,ngrid
          tauref(ig) = tauref(ig) + 
     &                aerosol(ig,l,iaerdust(iaer))
          ENDDO
        ENDDO
       ENDDO
      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
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