Index: trunk/LMDZ.TITAN/libf/phytitan/aerave.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/aerave.F	(revision 1787)
+++ 	(revision )
@@ -1,182 +1,0 @@
-      subroutine aerave( ndata,
-     &     longdata,epdata,omegdata,gdata,          
-     &     longref,epref,temp,nir,longir
-     &     ,epir,omegir,gir,qref )
-
-
-      implicit none
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Calculate mean values of aerosol quantities in each band.
-!     
-!     Authors
-!     ------- 
-!     R. Fournier (1996) 
-!     F. Forget (1996)
-!     
-!     Called by
-!     ---------
-!     suaer_corrk.F90
-!     
-!     Calls
-!     -----
-!     blackl.F
-!     
-!==================================================================
-
-!     
-!     R.Fournier 02/1996 
-!     (modif F.Forget 02/1996)
-!     le spectre est decoupe en "nir" bandes et cette routine calcule
-!     les donnees radiatives moyenne sur chaque bande : l'optimisation
-!     est faite pour une temperature au sol "temp" et une epaisseur
-!     optique de l'atmosphere "epref" a la longueur d'onde "longref"
-!     
-!     dans la version actuelle, les ponderations sont independantes de
-!     l'epaisseur optique : c'est a dire que "omegir", "gir"
-!     et "epir/epre" sont independants de "epref".
-!     en effet les ponderations sont choisies pour une solution exacte
-!     en couche mince et milieu isotherme. 
-!     
-!     entree
-!     
-!     ndata : taille des champs data
-!     longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol
-!     (longdata longueur d'onde en METRES)
-!     * longref : longueur d'onde a laquelle l'epaisseur optique
-!     est connue
-!     * epref : epaisseur optique a longref
-!     * temp : temperature choisie pour la ponderation (Planck)
-!     * nir : nombre d'intervals dans la discretisation spectrale
-!     du GCM
-!     * longir : longueurs d'onde definissant ces intervals
-!     
-!     sortie
-!     
-!     * epir : epaisseur optique moyenne pour chaque interval
-!     * omegir : "scattering albedo" moyen pour chaque interval
-!     * gir : "assymetry factor" moyen pour chaque interval
-!     * qref : extinction coefficient at reference wavelength
-
-!=======================================================================
-! output
-
-      REAL longref
-      REAL epref
-      REAL temp
-      INTEGER nir
-      REAL*8 longir(nir+1)
-      REAL epir(nir)
-      REAL omegir(nir)
-      REAL gir(nir)
-
-!=======================================================================
-
-      INTEGER iir,nirmx
-      PARAMETER (nirmx=100)
-      INTEGER idata,ndata
-
-      DOUBLE PRECISION tmp1
-      REAL tmp2,tmp3
-
-!=======================================================================
-! input
-
-      REAL emit
-      REAL totalemit(nirmx)
-      REAL longdata(ndata),epdata(ndata)
-     &     ,omegdata(ndata),gdata(ndata)
-      REAL qextcorrdata(ndata)
-      INTEGER ibande,nbande
-      PARAMETER (nbande=1000)
-
-      REAL long,deltalong
-      INTEGER ilong
-      INTEGER i1,i2
-      REAL c1,c2
-      REAL factep,qextcorr,omeg,g,qref
-
-      long=longref
-
-
-!=======================================================================
-!     pre-interpolation
-      ilong=1
-      DO idata=2,ndata
-         IF (long.gt.longdata(idata)) ilong=idata
-      ENDDO
-      i1=ilong
-      i2=ilong+1
-      IF (i2.gt.ndata) i2=ndata
-      IF (long.lt.longdata(1)) i2=1
-      IF (i1.eq.i2) THEN
-         c1=1.E+0
-         c2=0.E+0
-      ELSE
-         c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
-         c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
-      ENDIF
-
-      qref=c1*epdata(i1)+c2*epdata(i2)
-      factep=qref/epref
-      DO idata=1,ndata
-         qextcorrdata(idata)=epdata(idata)/factep
-      ENDDO
-!=======================================================================
-
-      DO iir=1,nir
-
-         deltalong=(longir(iir+1)-longir(iir)) / nbande
-         totalemit(iir)=0.E+0
-         epir(iir)=0.E+0
-         omegir(iir)=0.E+0
-         gir(iir)=0.E+0
-
-         DO ibande=1,nbande
-
-            long=longir(iir) + (ibande-0.5E+0) * deltalong
-            CALL blackl(DBLE(long),DBLE(temp),tmp1)
-            emit=REAL(tmp1)
-
-!=======================================================================
-!     interpolation
-            ilong=1
-            DO idata=2,ndata
-               IF (long.gt.longdata(idata)) ilong=idata
-            ENDDO
-            i1=ilong
-            i2=ilong+1
-
-            IF (i2.gt.ndata) i2=ndata
-            IF (long.lt.longdata(1)) i2=1
-            IF (i1.eq.i2) THEN
-               c1=1.E+0
-               c2=0.E+0
-            ELSE
-               c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
-               c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
-            ENDIF
-!=======================================================================
-
-            qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2)
-            omeg=c1*omegdata(i1)+c2*omegdata(i2)
-            g=c1*gdata(i1)+c2*gdata(i2)
-
-            totalemit(iir)=totalemit(iir)+deltalong*emit
-            epir(iir)=epir(iir)+deltalong*emit*qextcorr
-            omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr
-            gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g
-
-         ENDDO
-
-         gir(iir)=gir(iir)/omegir(iir)
-         omegir(iir)=omegir(iir)/epir(iir)
-         epir(iir)=epir(iir)/totalemit(iir)
-
-      ENDDO
-
-      return
-      end
Index: trunk/LMDZ.TITAN/libf/phytitan/aerave_new.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/aerave_new.F	(revision 1787)
+++ 	(revision )
@@ -1,223 +1,0 @@
-      SUBROUTINE aerave_new ( ndata,
-     & longdata,epdata,omegdata,gdata,          
-     &            longref,epref,temp,nir,longir
-     &            ,epir,omegir,gir,qref,omegaref        )
-
-
-      IMPLICIT NONE
-c.......................................................................
-c
-c R.Fournier 02/1996 
-c (modif F.Forget 02/1996)
-c le spectre est decoupe en "nir" bandes et cette routine calcule
-c les donnees radiatives moyenne sur chaque bande : l'optimisation
-c est faite pour une temperature au sol "temp" et une epaisseur
-c optique de l'atmosphere "epref" a la longueur d'onde "longref"
-c
-c dans la version actuelle, les ponderations sont independantes de
-c l'epaisseur optique : c'est a dire que "omegir", "gir"
-c et "epir/epre" sont independants de "epref".
-c en effet les ponderations sont choisies pour une solution exacte
-c en couche mince et milieu isotherme. 
-c
-c entree
-c
-c    ndata : taille des champs data
-c    longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol
-c                  (longdata longueur d'onde en METRES)
-c  * longref : longueur d'onde a laquelle l'epaisseur optique
-c              est connue
-c  * epref : epaisseur optique a longref
-c  * temp : temperature choisie pour la ponderation (Planck)
-c  * nir : nombre d'intervals dans la discretisation spectrale
-c           du GCM
-c  * longir : longueurs d'onde definissant ces intervals
-c
-c sortie
-c
-c  * epir : epaisseur optique moyenne pour chaque interval
-c  * omegir : "scattering albedo" moyen pour chaque interval
-c  * gir : "assymetry factor" moyen pour chaque interval
-c  * qref : extinction coefficient at reference wavelength
-c  * omegaref : single scat. albedo at reference wavelength
-c
-c.......................................................................
-c
-      REAL longref
-      REAL epref
-      REAL temp
-      INTEGER nir
-      REAL*8 longir(nir+1)
-      REAL epir(nir)
-      REAL omegir(nir)
-      REAL gir(nir)
-c
-c.......................................................................
-c
-      INTEGER iir,nirmx
-      PARAMETER (nirmx=100)
-      INTEGER idata,ndata
-c
-c.......................................................................
-c
-      REAL emit
-      REAL totalemit(nirmx)
-      REAL longdata(ndata),epdata(ndata)
-     &    ,omegdata(ndata),gdata(ndata)
-      REAL qextcorrdata(ndata)
-      INTEGER ibande,nbande
-      PARAMETER (nbande=1000)
-      REAL long,deltalong
-      INTEGER ilong
-      INTEGER i1,i2
-      REAL c1,c2
-      REAL factep,qextcorr,omeg,g
-      REAL qref,omegaref
-c
-c.......................................................................
-c
-      DOUBLE PRECISION tmp1
-      REAL tmp2,tmp3
-c
-c
-      long=longref
-
-
-      !if(nir.eq.27)then
-      !print*,'long',long
-      !print*,'longdata',longdata
-      !print*,'epdata',epdata
-      !print*,'omegdata',omegdata
-      !print*,'gdata',gdata
-      !print*,'data looks aok!'
-
-      !print*,'ndata=',ndata
-      !print*,'longdata',shape(longdata)
-      !print*,'epdata',shape(epdata)
-      !print*,'omegdata',shape(omegdata)
-      !print*,'gdata',shape(gdata)
-      ! print*,'longref',longref
-      !print*,'epref',epref
-      !print*,'temp',temp
-      !print*,'nir',nir
-      !print*,'longir',longir
-      !print*,'epir',epir
-      !print*,'omegir',gir
-      !print*,'qref',qref
-      !print*,'longir=',longir
-      !stop
-      !endif
-
-
-c********************************************************
-c interpolation
-      ilong=1
-      DO idata=2,ndata
-        IF (long.gt.longdata(idata)) ilong=idata
-      ENDDO
-      i1=ilong
-      i2=ilong+1
-      IF (i2.gt.ndata) i2=ndata
-      IF (long.lt.longdata(1)) i2=1
-      IF (i1.eq.i2) THEN
-        c1=1.E+0
-        c2=0.E+0
-      ELSE
-        c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
-        c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
-      ENDIF
-c********************************************************
-c
-      qref=c1*epdata(i1)+c2*epdata(i2)
-      omegaref=c1*omegdata(i1)+c2*omegdata(i2)
-      factep=qref/epref
-      DO idata=1,ndata
-        qextcorrdata(idata)=epdata(idata)/factep
-      ENDDO
-c
-c.......................................................................
-c
-      DO iir=1,nir
-c
-c.......................................................................
-c
-        deltalong=(longir(iir+1)-longir(iir)) / nbande
-        totalemit(iir)=0.E+0
-        epir(iir)=0.E+0
-        omegir(iir)=0.E+0
-        gir(iir)=0.E+0
-c
-c.......................................................................
-c
-        DO ibande=1,nbande
-c
-c.......................................................................
-c
-          long=longir(iir) + (ibande-0.5E+0) * deltalong
-          CALL blackl(DBLE(long),DBLE(temp),tmp1)
-          emit=REAL(tmp1)
-c
-c.......................................................................
-c
-c********************************************************
-c interpolation
-      ilong=1
-      DO idata=2,ndata
-        IF (long.gt.longdata(idata)) ilong=idata
-      ENDDO
-      i1=ilong
-      i2=ilong+1
-      IF (i2.gt.ndata) i2=ndata
-      IF (long.lt.longdata(1)) i2=1
-      IF (i1.eq.i2) THEN
-        c1=1.E+0
-        c2=0.E+0
-      ELSE
-        c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1))
-        c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2))
-      ENDIF
-c********************************************************
-c
-          qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2)
-          omeg=c1*omegdata(i1)+c2*omegdata(i2)
-          g=c1*gdata(i1)+c2*gdata(i2)
-c
-c.......................................................................
-c
-          totalemit(iir)=totalemit(iir)+deltalong*emit
-          epir(iir)=epir(iir)+deltalong*emit*qextcorr
-          omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr
-          gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g
-c
-c.......................................................................
-c
-        ENDDO
-c
-c.......................................................................
-c
-        gir(iir)=gir(iir)/omegir(iir)
-        omegir(iir)=omegir(iir)/epir(iir)
-        epir(iir)=epir(iir)/totalemit(iir)
-c
-c.......................................................................
-c
-      ENDDO
-c
-c......................................................................
-c
-c     Diagnostic de controle si on moyenne sur tout le spectre vis ou IR :
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c     tmp2=0.E+0
-c     DO iir=1,nir
-c       tmp2=tmp2+totalemit(iir)
-c     ENDDO
-c     tmp3=5.67E-8 * temp**4
-c     IF (abs((tmp2-tmp3)/tmp3).gt.0.05E+0) THEN
-c       PRINT *,'!!!! <---> il manque du Planck (voir moyenne.F)'
-c       PRINT *,'somme des bandes :',tmp2,'--- Planck:',tmp3
-c     ENDIF
-c
-c......................................................................
-c
-      RETURN
-      END
Index: trunk/LMDZ.TITAN/libf/phytitan/aeropacity.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/aeropacity.F90	(revision 1787)
+++ 	(revision )
@@ -1,184 +1,0 @@
-      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
-         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col)
-
-       use radinc_h, only : L_TAUMAX,naerkind
-       use aerosol_mod
-       USE tracer_h, only: noms
-       use comcstfi_mod, only: g
-       use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,  &
-		pres_bottom_strato,pres_top_strato,obs_tau_col_strato
-                  
-       implicit none
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Compute aerosol optical depth in each gridbox.
-!     
-!     Authors
-!     ------- 
-!     F. Forget
-!     F. Montmessin (water ice scheme) 
-!     update J.-B. Madeleine (2008)
-!     dust removal, simplification by Robin Wordsworth (2009)
-!
-!     Input
-!     ----- 
-!     ngrid             Number of horizontal gridpoints
-!     nlayer            Number of layers
-!     nq                Number of tracers
-!     pplev             Pressure (Pa) at each layer boundary
-!     pq                Aerosol mixing ratio
-!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
-!     QREFvis3d(ngrid,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
-
-      real aerosol0
-
-      INTEGER l,ig,iq,iaer
-
-      LOGICAL,SAVE :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-      REAL CBRT
-      EXTERNAL CBRT
-
-      ! for fixed dust profiles
-      real expfactor, zp
-
-      ! identify tracers
-      IF (firstcall) THEN
-
-        write(*,*) "Tracers found in aeropacity:"
-
-        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_back2lay.ne.0) then
-          print*,'iaero_back2lay= ',iaero_back2lay
-        endif
-
-        firstcall=.false.
-      ENDIF ! of IF (firstcall)
-
-
-!     ---------------------------------------------------------   
-
-      if (noaero) then
-        aerosol(1:ngrid,1:nlayer,1)=0.0 ! JVO 2017 : Now iaer = 1 is always dummy co2 for noaero case, since we don't use aerosols
-      endif
-
-!==================================================================
-!    Two-layer aerosols (unknown composition)
-!    S. Guerlet (2013)
-!==================================================================
-
-      if (iaero_back2lay .ne.0) then
-           iaer=iaero_back2lay
-!       1. Initialization
-            aerosol(1:ngrid,1:nlayer,iaer)=0.0
-!       2. Opacity calculation
-          DO ig=1,ngrid
-           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.*aerosol(ig,l,iaer)
-             !! 2. tropo layer
-             ELSEIF (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)
-             !! 3. linear transition
-             ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
-               expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
-               aerosol(ig,l,iaer)= obs_tau_col_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)*aerosol(ig,l,iaer)/1.5
-             !! 4. strato layer
-             ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
-               aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
-             !! 5. above strato layer: no aerosols
-             ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
-               aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
-             ENDIF
-	   ENDDO
-          ENDDO
-
- !       3. Re-normalize to observed total column
-         tau_col(:)=0.0
-         DO l=1,nlayer
-          DO ig=1,ngrid
-               tau_col(ig) = tau_col(ig) &
-                     + aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
-            ENDDO
-         ENDDO
-
-         DO ig=1,ngrid
-           DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
-           ENDDO
-         ENDDO
-
-
-      end if ! if Two-layer 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
-      
Index: trunk/LMDZ.TITAN/libf/phytitan/aeroptproperties.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/aeroptproperties.F90	(revision 1787)
+++ 	(revision )
@@ -1,754 +1,0 @@
-      SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad,   &
-                                  QVISsQREF3d,omegaVIS3d,gVIS3d,   &
-                                  QIRsQREF3d,omegaIR3d,gIR3d,      &
-                                  QREFvis3d,QREFir3d)!,		   &
-!                                  omegaREFvis3d,omegaREFir3d)
-
-      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
-      use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
-      use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
-      use radcommon_h, only: radiustab,nsize
-
-      implicit none
-
-!     =============================================================
-!     Aerosol Optical Properties
-!
-!     Description:
-!       Compute the scattering parameters in each grid
-!       box, depending on aerosol grain sizes. Log-normal size
-!       distribution and Gauss-Legendre integration are used.
-
-!     Parameters:
-!       Don't forget to set the value of varyingnueff below; If
-!       the effective variance of the distribution for the given
-!       aerosol is considered homogeneous in the atmosphere, please
-!       set varyingnueff(iaer) to .false. Resulting computational
-!       time will be much better.
-
-!     Authors: J.-B. Madeleine, F. Forget, F. Montmessin
-!     Slightly modified and converted to F90 by R. Wordsworth (2009)
-!     Varying nueff section removed by R. Wordsworth for simplicity
-!     ==============================================================
-
-!     Local variables 
-!     ---------------
-
-
-
-!     =============================================================
-      LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
-!     =============================================================
-
-!     Min. and max radius of the interpolation grid (in METERS)
-      REAL, PARAMETER :: refftabmin = 2e-8 !2e-8
-!      REAL, PARAMETER :: refftabmax = 35e-6
-      REAL, PARAMETER :: refftabmax = 1e-3
-!     Log of the min and max variance of the interpolation grid
-      REAL, PARAMETER :: nuefftabmin = -4.6
-      REAL, PARAMETER :: nuefftabmax = 0.
-!     Number of effective radius of the interpolation grid
-      INTEGER, PARAMETER :: refftabsize = 200
-!     Number of effective variances of the interpolation grid
-!      INTEGER, PARAMETER :: nuefftabsize = 100
-      INTEGER, PARAMETER :: nuefftabsize = 1
-!     Interpolation grid indices (reff,nueff)
-      INTEGER :: grid_i,grid_j
-!     Intermediate variable
-      REAL :: var_tmp,var3d_tmp(ngrid,nlayer)
-!     Bilinear interpolation factors
-      REAL :: kx,ky,k1,k2,k3,k4
-!     Size distribution parameters
-      REAL :: sizedistk1,sizedistk2
-!     Pi!
-      REAL,SAVE :: pi
-!$OMP THREADPRIVATE(pi)
-!     Variables used by the Gauss-Legendre integration:
-      INTEGER radius_id,gausind
-      REAL kint
-      REAL drad
-      INTEGER, PARAMETER :: ngau = 10
-      REAL weightgaus(ngau),radgaus(ngau)
-      SAVE weightgaus,radgaus
-!      DATA weightgaus/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/
-!      DATA radgaus/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/
-      DATA    radgaus/0.07652652113350,0.22778585114165, &
-                      0.37370608871528,0.51086700195146, &
-                      0.63605368072468,0.74633190646476, &
-                      0.83911697181213,0.91223442826796, &
-                      0.96397192726078,0.99312859919241/
-
-      DATA weightgaus/0.15275338723120,0.14917298659407, &
-                      0.14209610937519,0.13168863843930, &
-                      0.11819453196154,0.10193011980823, &
-                      0.08327674160932,0.06267204829828, &
-                      0.04060142982019,0.01761400714091/
-!$OMP THREADPRIVATE(radgaus,weightgaus)
-!     Indices
-      INTEGER :: i,j,k,l,m,iaer,idomain
-      INTEGER :: ig,lg,chg
-
-!     Local saved variables
-!     ---------------------
-
-!     Radius axis of the interpolation grid
-      REAL,SAVE :: refftab(refftabsize)
-!     Variance axis of the interpolation grid
-      REAL,SAVE :: nuefftab(nuefftabsize)
-!     Volume ratio of the grid
-      REAL,SAVE :: logvratgrid,vratgrid
-!     Grid used to remember which calculation is done
-      LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2) = .false.
-!$OMP THREADPRIVATE(refftab,nuefftab,logvratgrid,vratgrid,checkgrid)
-!     Optical properties of the grid (VISIBLE)
-      REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
-      REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
-      REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
-      REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
-      REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
-!$OMP THREADPRIVATE(qsqrefVISgrid,qextVISgrid,qscatVISgrid,omegVISgrid,gVISgrid)
-!     Optical properties of the grid (INFRARED)
-      REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
-      REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
-      REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
-      REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
-      REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
-!$OMP THREADPRIVATE(qsqrefIRgrid,qextIRgrid,qscatIRgrid,omegIRgrid,gIRgrid)
-!     Optical properties of the grid (REFERENCE WAVELENGTHS)
-      REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
-      REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
-      REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
-      REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
-      REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
-      REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
-!$OMP THREADPRIVATE(qrefVISgrid,qscatrefVISgrid,qrefIRgrid,qscatrefIRgrid,omegrefVISgrid,&
-	!$OMP omegrefIRgrid)
-!     Firstcall
-      LOGICAL,SAVE :: firstcall = .true.
-!$OMP THREADPRIVATE(firstcall)
-!     Variables used by the Gauss-Legendre integration:
-      REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
-      REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
-      REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
-!$OMP THREADPRIVATE(normd,dista,distb)
-
-      REAL,SAVE :: radGAUSa(ngau,naerkind,2)
-      REAL,SAVE :: radGAUSb(ngau,naerkind,2)
-!$OMP THREADPRIVATE(radGAUSa,radGAUSb)
-
-      REAL,SAVE :: qsqrefVISa(L_NSPECTV,ngau,naerkind)
-      REAL,SAVE :: qrefVISa(ngau,naerkind)
-      REAL,SAVE :: qsqrefVISb(L_NSPECTV,ngau,naerkind)
-      REAL,SAVE :: qrefVISb(ngau,naerkind)
-      REAL,SAVE :: omegVISa(L_NSPECTV,ngau,naerkind)
-      REAL,SAVE :: omegrefVISa(ngau,naerkind)
-      REAL,SAVE :: omegVISb(L_NSPECTV,ngau,naerkind)
-      REAL,SAVE :: omegrefVISb(ngau,naerkind)
-      REAL,SAVE :: gVISa(L_NSPECTV,ngau,naerkind)
-      REAL,SAVE :: gVISb(L_NSPECTV,ngau,naerkind)
-!$OMP THREADPRIVATE(qsqrefVISa,qrefVISa,qsqrefVISb,qrefVISb,omegVISa, &
-	!$OMP omegrefVISa,omegVISb,omegrefVISb,gVISa,gVISb)
-
-      REAL,SAVE :: qsqrefIRa(L_NSPECTI,ngau,naerkind)
-      REAL,SAVE :: qrefIRa(ngau,naerkind)
-      REAL,SAVE :: qsqrefIRb(L_NSPECTI,ngau,naerkind)
-      REAL,SAVE :: qrefIRb(ngau,naerkind)
-      REAL,SAVE :: omegIRa(L_NSPECTI,ngau,naerkind)
-      REAL,SAVE :: omegrefIRa(ngau,naerkind)
-      REAL,SAVE :: omegIRb(L_NSPECTI,ngau,naerkind)
-      REAL,SAVE :: omegrefIRb(ngau,naerkind)
-      REAL,SAVE :: gIRa(L_NSPECTI,ngau,naerkind)
-      REAL,SAVE :: gIRb(L_NSPECTI,ngau,naerkind)
-!$OMP THREADPRIVATE(qsqrefIRa,qrefIRa,qsqrefIRb,qrefIRb,omegIRa,omegrefIRa,&
-	!$OMP omegIRb,omegrefIRb,gIRa,gIRb)
-
-      REAL :: radiusm
-      REAL :: radiusr
-
-!     Inputs
-!     ------
-
-      INTEGER :: ngrid,nlayer
-!     Aerosol effective radius used for radiative transfer (meter)
-      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind)
-!     Aerosol effective variance used for radiative transfer (n.u.)
-      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind)
-
-!     Outputs
-!     -------
-
-      REAL,INTENT(OUT) :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
-      REAL,INTENT(OUT) :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
-      REAL,INTENT(OUT) :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
-
-      REAL,INTENT(OUT) :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
-      REAL,INTENT(OUT) :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
-      REAL,INTENT(OUT) :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
-
-      REAL,INTENT(OUT) :: QREFvis3d(ngrid,nlayer,naerkind)
-      REAL,INTENT(OUT) :: QREFir3d(ngrid,nlayer,naerkind)
-
-      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
-      REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
-
-      DO iaer = 1, naerkind ! Loop on aerosol kind
-        IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
-!==================================================================
-!       If there is one single particle size, optical
-!         properties of the considered aerosol are homogeneous
-          DO lg = 1, nlayer
-            DO ig = 1, ngrid
-              DO chg = 1, L_NSPECTV
-                QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
-                omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
-                gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
-              ENDDO
-              DO chg = 1, L_NSPECTI
-                QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
-                omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
-                gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
-              ENDDO
-              QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
-              QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
-              omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
-              omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
-            ENDDO
-          ENDDO
-
-
-          if (firstcall) then
-             print*,'Optical prop. of the aerosol are homogenous for:'
-             print*,'iaer = ',iaer
-          endif
-
-        ELSE ! Varying effective radius and variance
-      DO idomain = 1, 2 ! Loop on visible or infrared channel
-!==================================================================
-!     1. Creating the effective radius and variance grid
-!     --------------------------------------------------
-      IF (firstcall) THEN
-
-!       1.1 Pi!
-        pi = 2. * asin(1.e0)
-
-!       1.2 Effective radius
-        refftab(1)    = refftabmin
-        refftab(refftabsize) = refftabmax
-
-        logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3.
-        vratgrid = exp(logvratgrid)
-
-        do i = 2, refftabsize-1
-          refftab(i) = refftab(i-1)*vratgrid**(1./3.)
-        enddo
-
-!       1.3 Effective variance
-        if(nuefftabsize.eq.1)then ! addded by RDW
-           print*,'Warning: no variance range in aeroptproperties'
-           nuefftab(1)=0.2
-        else
-           do i = 0, nuefftabsize-1
-              nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
-           enddo
-        endif
-
-        firstcall = .false.
-      ENDIF
-
-!       1.4 Radius middle point and range for Gauss integration
-        radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1))
-        radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1))
-
-!       1.5 Interpolating data at the Gauss quadrature points:
-        DO gausind=1,ngau
-          drad=radiusr*radgaus(gausind)
-          radGAUSa(gausind,iaer,idomain)=radiusm-drad
-
-          radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1)
-          IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN
-            radius_id=radius_id-1
-          ENDIF
-          IF (radius_id.GE.nsize(iaer,idomain)) THEN
-            radius_id=nsize(iaer,idomain)-1
-            kint = 1.
-          ELSEIF (radius_id.LT.1) THEN
-            radius_id=1
-            kint = 0.
-          ELSE
-          kint = ( (radiusm-drad) -				&
-                   radiustab(iaer,idomain,radius_id) ) /	&
-                 ( radiustab(iaer,idomain,radius_id+1) -	&
-                   radiustab(iaer,idomain,radius_id) )
-          ENDIF
-          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
-            DO m=1,L_NSPECTV
-               qsqrefVISa(m,gausind,iaer)=                      &
-                    (1-kint)*QVISsQREF(m,iaer,radius_id) +      &
-                    kint*QVISsQREF(m,iaer,radius_id+1)
-            omegVISa(m,gausind,iaer)=                           &
-                    (1-kint)*omegaVIS(m,iaer,radius_id) +       &
-                    kint*omegaVIS(m,iaer,radius_id+1)
-            gVISa(m,gausind,iaer)=                              &
-                    (1-kint)*gVIS(m,iaer,radius_id) +           &
-                    kint*gVIS(m,iaer,radius_id+1)
-            ENDDO
-            qrefVISa(gausind,iaer)=                             &
-                    (1-kint)*QREFvis(iaer,radius_id) +          &
-                    kint*QREFvis(iaer,radius_id+1)
-            omegrefVISa(gausind,iaer)=                          &
-                    (1-kint)*omegaREFvis(iaer,radius_id) +      &
-                    kint*omegaREFvis(iaer,radius_id+1)
-          ELSE ! INFRARED DOMAIN ----------------------------------
-            DO m=1,L_NSPECTI
-            qsqrefIRa(m,gausind,iaer)=                          &
-                    (1-kint)*QIRsQREF(m,iaer,radius_id) +       &
-                    kint*QIRsQREF(m,iaer,radius_id+1)
-            omegIRa(m,gausind,iaer)=                            &
-                    (1-kint)*omegaIR(m,iaer,radius_id) +        &
-                    kint*omegaIR(m,iaer,radius_id+1)
-            gIRa(m,gausind,iaer)=                               &
-                    (1-kint)*gIR(m,iaer,radius_id) +            &
-                    kint*gIR(m,iaer,radius_id+1)
-            ENDDO
-            qrefIRa(gausind,iaer)=                              &
-                    (1-kint)*QREFir(iaer,radius_id) +           &
-                    kint*QREFir(iaer,radius_id+1)
-            omegrefIRa(gausind,iaer)=                           &
-                    (1-kint)*omegaREFir(iaer,radius_id) +       &
-                    kint*omegaREFir(iaer,radius_id+1)
-          ENDIF
-        ENDDO
-
-        DO gausind=1,ngau
-          drad=radiusr*radgaus(gausind)
-          radGAUSb(gausind,iaer,idomain)=radiusm+drad
-
-          radius_id=minloc(abs(radiustab(iaer,idomain,:) -      &
-                               (radiusm+drad)),DIM=1)
-          IF ((radiustab(iaer,idomain,radius_id) -              &
-               (radiusm+drad)).GT.0) THEN
-            radius_id=radius_id-1
-          ENDIF
-          IF (radius_id.GE.nsize(iaer,idomain)) THEN
-            radius_id=nsize(iaer,idomain)-1
-            kint = 1.
-          ELSEIF (radius_id.LT.1) THEN
-            radius_id=1
-            kint = 0.
-          ELSE
-            kint = ( (radiusm+drad) -                           &
-                     radiustab(iaer,idomain,radius_id) ) /      &
-                   ( radiustab(iaer,idomain,radius_id+1) -      &
-                     radiustab(iaer,idomain,radius_id) )
-          ENDIF
-          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
-            DO m=1,L_NSPECTV
-            qsqrefVISb(m,gausind,iaer)=                         &
-                    (1-kint)*QVISsQREF(m,iaer,radius_id) +      &
-                    kint*QVISsQREF(m,iaer,radius_id+1)    
-            omegVISb(m,gausind,iaer)=                           &
-                    (1-kint)*omegaVIS(m,iaer,radius_id) +       &
-                    kint*omegaVIS(m,iaer,radius_id+1)
-            gVISb(m,gausind,iaer)=                              &
-                    (1-kint)*gVIS(m,iaer,radius_id) +           &
-                    kint*gVIS(m,iaer,radius_id+1)
-            ENDDO
-            qrefVISb(gausind,iaer)=                             &
-                    (1-kint)*QREFvis(iaer,radius_id) +          &
-                    kint*QREFvis(iaer,radius_id+1)
-            omegrefVISb(gausind,iaer)=                          &
-                    (1-kint)*omegaREFvis(iaer,radius_id) +      &
-                    kint*omegaREFvis(iaer,radius_id+1)
-          ELSE ! INFRARED DOMAIN ----------------------------------
-            DO m=1,L_NSPECTI
-            qsqrefIRb(m,gausind,iaer)=                          &
-                    (1-kint)*QIRsQREF(m,iaer,radius_id) +       &
-                    kint*QIRsQREF(m,iaer,radius_id+1)
-            omegIRb(m,gausind,iaer)=                            &
-                    (1-kint)*omegaIR(m,iaer,radius_id) +        &
-                    kint*omegaIR(m,iaer,radius_id+1)
-            gIRb(m,gausind,iaer)=                               &
-                    (1-kint)*gIR(m,iaer,radius_id) +            &
-                    kint*gIR(m,iaer,radius_id+1)
-            ENDDO
-            qrefIRb(gausind,iaer)=                              &
-                    (1-kint)*QREFir(iaer,radius_id) +           &
-                    kint*QREFir(iaer,radius_id+1)
-            omegrefIRb(gausind,iaer)=                           &
-                    (1-kint)*omegaREFir(iaer,radius_id) +       &
-                    kint*omegaREFir(iaer,radius_id+1)
-          ENDIF
-        ENDDO
-
-!==================================================================
-! CONSTANT NUEFF FROM HERE
-!==================================================================
-
-!     2. Compute the scattering parameters using linear
-!       interpolation over grain sizes and constant nueff
-!     ---------------------------------------------------
-
-      DO lg = 1,nlayer
-        DO ig = 1, ngrid
-!         2.1 Effective radius index and kx calculation
-          var_tmp=reffrad(ig,lg,iaer)/refftabmin
-          var_tmp=log(var_tmp)*3.
-          var_tmp=var_tmp/logvratgrid+1.
-          grid_i=floor(var_tmp)
-          IF (grid_i.GE.refftabsize) THEN
-!           WRITE(*,*) 'Warning: particle size in grid box #'
-!           WRITE(*,*) ig,' is too large to be used by the '
-!           WRITE(*,*) 'radiative transfer; please extend the '
-!           WRITE(*,*) 'interpolation grid to larger grain sizes.'
-            grid_i=refftabsize-1
-            kx = 1.
-          ELSEIF (grid_i.LT.1) THEN
-!           WRITE(*,*) 'Warning: particle size in grid box #'
-!           WRITE(*,*) ig,' is too small to be used by the '
-!           WRITE(*,*) 'radiative transfer; please extend the '
-!           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
-            grid_i=1
-            kx = 0.
-          ELSE
-            kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /            &
-                 ( refftab(grid_i+1)-refftab(grid_i) )
-          ENDIF
-!         2.3 Integration
-          DO j=grid_i,grid_i+1
-!             2.3.1 Check if the calculation has been done
-              IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
-!               2.3.2 Log-normal dist., r_g and sigma_g are defined
-!                 in [hansen_1974], "Light scattering in planetary
-!                 atmospheres", Space Science Reviews 16 527-610.
-!                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
-                sizedistk2 = log(1.+nueffrad(1,1,iaer))
-                sizedistk1 = exp(2.5*sizedistk2)
-                sizedistk1 = refftab(j) / sizedistk1
-
-                normd(j,1,iaer,idomain) = 1e-30
-                DO gausind=1,ngau
-                  drad=radiusr*radgaus(gausind)
-                  dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
-                  dista(j,1,iaer,idomain,gausind) =                   &
-                    EXP(-dista(j,1,iaer,idomain,gausind) *            &
-                    dista(j,1,iaer,idomain,gausind) *                 &
-                    0.5e0/sizedistk2)/(radiusm-drad)                  
-                  dista(j,1,iaer,idomain,gausind) =                   &
-                    dista(j,1,iaer,idomain,gausind) /                 &
-                    (sqrt(2e0*pi*sizedistk2))
-
-                  distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
-                  distb(j,1,iaer,idomain,gausind) =                   &
-                    EXP(-distb(j,1,iaer,idomain,gausind) *            &
-                    distb(j,1,iaer,idomain,gausind) *                 &
-                    0.5e0/sizedistk2)/(radiusm+drad)
-                  distb(j,1,iaer,idomain,gausind) =                   &
-                    distb(j,1,iaer,idomain,gausind) /                 &
-                    (sqrt(2e0*pi*sizedistk2))
-
-                  normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +   &
-                    weightgaus(gausind) *                             &
-                    (                                                 &
-                    distb(j,1,iaer,idomain,gausind) * pi *            &
-                    radGAUSb(gausind,iaer,idomain) *                  &
-                    radGAUSb(gausind,iaer,idomain) +                  &
-                    dista(j,1,iaer,idomain,gausind) * pi *            &
-                    radGAUSa(gausind,iaer,idomain) *                  &
-                    radGAUSa(gausind,iaer,idomain)                    &
-                    )
-                ENDDO
-                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
-!                 2.3.3.vis Initialization
-                  qsqrefVISgrid(j,1,:,iaer)=0.
-                  qextVISgrid(j,1,:,iaer)=0.
-                  qscatVISgrid(j,1,:,iaer)=0.
-                  omegVISgrid(j,1,:,iaer)=0.
-                  gVISgrid(j,1,:,iaer)=0.
-                  qrefVISgrid(j,1,iaer)=0.
-                  qscatrefVISgrid(j,1,iaer)=0.
-                  omegrefVISgrid(j,1,iaer)=0.
-
-                  DO gausind=1,ngau
-                    DO m=1,L_NSPECTV
-!                     Convolution:
-                      qextVISgrid(j,1,m,iaer) =              &
-                        qextVISgrid(j,1,m,iaer) +            & 
-                        weightgaus(gausind) *                &
-                        (                                    &
-                        qsqrefVISb(m,gausind,iaer) *         &
-                        qrefVISb(gausind,iaer) *             &
-                        pi*radGAUSb(gausind,iaer,idomain) *  &
-                        radGAUSb(gausind,iaer,idomain) *     &
-                        distb(j,1,iaer,idomain,gausind) +    &
-                        qsqrefVISa(m,gausind,iaer) *         &
-                        qrefVISa(gausind,iaer) *             &
-                        pi*radGAUSa(gausind,iaer,idomain) *  &
-                        radGAUSa(gausind,iaer,idomain) *     &
-                        dista(j,1,iaer,idomain,gausind)      &
-                        )
-                      qscatVISgrid(j,1,m,iaer) =             &
-                        qscatVISgrid(j,1,m,iaer) +           &
-                        weightgaus(gausind) *                &
-                        (                                    &
-                        omegVISb(m,gausind,iaer) *           &
-                        qsqrefVISb(m,gausind,iaer) *         &
-                        qrefVISb(gausind,iaer) *             &
-                        pi*radGAUSb(gausind,iaer,idomain) *  &
-                        radGAUSb(gausind,iaer,idomain) *     &
-                        distb(j,1,iaer,idomain,gausind) +    &
-                        omegVISa(m,gausind,iaer) *           &
-                        qsqrefVISa(m,gausind,iaer) *         &
-                        qrefVISa(gausind,iaer) *             &
-                        pi*radGAUSa(gausind,iaer,idomain) *  &
-                        radGAUSa(gausind,iaer,idomain) *     &
-                        dista(j,1,iaer,idomain,gausind)      &
-                        )
-                      gVISgrid(j,1,m,iaer) =                 &
-                        gVISgrid(j,1,m,iaer) +               &
-                        weightgaus(gausind) *                &
-                        (                                    &
-                        omegVISb(m,gausind,iaer) *           &
-                        qsqrefVISb(m,gausind,iaer) *         &
-                        qrefVISb(gausind,iaer) *             &
-                        gVISb(m,gausind,iaer) *              &
-                        pi*radGAUSb(gausind,iaer,idomain) *  &
-                        radGAUSb(gausind,iaer,idomain) *     &
-                        distb(j,1,iaer,idomain,gausind) +    &
-                        omegVISa(m,gausind,iaer) *           &
-                        qsqrefVISa(m,gausind,iaer) *         &
-                        qrefVISa(gausind,iaer) *             &
-                        gVISa(m,gausind,iaer) *              &
-                        pi*radGAUSa(gausind,iaer,idomain) *  &
-                        radGAUSa(gausind,iaer,idomain) *     &
-                        dista(j,1,iaer,idomain,gausind)      &
-                        )
-                    ENDDO
-                    qrefVISgrid(j,1,iaer) =                  &
-                      qrefVISgrid(j,1,iaer) +                &
-                      weightgaus(gausind) *                  &
-                      (                                      &
-                      qrefVISb(gausind,iaer) *               &
-                      pi*radGAUSb(gausind,iaer,idomain) *    &
-                      radGAUSb(gausind,iaer,idomain) *       &
-                      distb(j,1,iaer,idomain,gausind) +      &
-                      qrefVISa(gausind,iaer) *               &
-                      pi*radGAUSa(gausind,iaer,idomain) *    &
-                      radGAUSa(gausind,iaer,idomain) *       &
-                      dista(j,1,iaer,idomain,gausind)        &
-                      )
-                    qscatrefVISgrid(j,1,iaer) =              &
-                      qscatrefVISgrid(j,1,iaer) +            &
-                      weightgaus(gausind) *                  &
-                      (                                      &
-                      omegrefVISb(gausind,iaer) *            &
-                      qrefVISb(gausind,iaer) *               & 
-                      pi*radGAUSb(gausind,iaer,idomain) *    &
-                      radGAUSb(gausind,iaer,idomain) *       &
-                      distb(j,1,iaer,idomain,gausind) +      &
-                      omegrefVISa(gausind,iaer) *            &
-                      qrefVISa(gausind,iaer) *               &
-                      pi*radGAUSa(gausind,iaer,idomain) *    &
-                      radGAUSa(gausind,iaer,idomain) *       &
-                      dista(j,1,iaer,idomain,gausind)        &
-                      )
-                  ENDDO
-
-                  qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /          &
-                                normd(j,1,iaer,idomain)       
-                  qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /  &
-                                normd(j,1,iaer,idomain)
-                  omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /   &
-                               qrefVISgrid(j,1,iaer)
-                  DO m=1,L_NSPECTV
-                    qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /    &
-                                normd(j,1,iaer,idomain)
-                    qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /  &
-                                normd(j,1,iaer,idomain)
-                    gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /          &
-                                qscatVISgrid(j,1,m,iaer) /               &
-                                normd(j,1,iaer,idomain)
-
-                    qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /  &
-                                qrefVISgrid(j,1,iaer)
-                    omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /   &
-                                qextVISgrid(j,1,m,iaer)
-                  ENDDO
-                ELSE                   ! INFRARED DOMAIN ----------
-!                 2.3.3.ir Initialization
-                  qsqrefIRgrid(j,1,:,iaer)=0.
-                  qextIRgrid(j,1,:,iaer)=0.
-                  qscatIRgrid(j,1,:,iaer)=0.
-                  omegIRgrid(j,1,:,iaer)=0.
-                  gIRgrid(j,1,:,iaer)=0.
-                  qrefIRgrid(j,1,iaer)=0.
-                  qscatrefIRgrid(j,1,iaer)=0.
-                  omegrefIRgrid(j,1,iaer)=0.
-
-                  DO gausind=1,ngau
-                    DO m=1,L_NSPECTI
-!                     Convolution:
-                      qextIRgrid(j,1,m,iaer) =                  &
-                        qextIRgrid(j,1,m,iaer) +                &
-                        weightgaus(gausind) *                   &
-                        (                                       &
-                        qsqrefIRb(m,gausind,iaer) *             &
-                        qrefVISb(gausind,iaer) *                &
-                        pi*radGAUSb(gausind,iaer,idomain) *     &
-                        radGAUSb(gausind,iaer,idomain) *        &
-                        distb(j,1,iaer,idomain,gausind) +       &
-                        qsqrefIRa(m,gausind,iaer) *             &
-                        qrefVISa(gausind,iaer) *                &
-                        pi*radGAUSa(gausind,iaer,idomain) *     &
-                        radGAUSa(gausind,iaer,idomain) *        &
-                        dista(j,1,iaer,idomain,gausind)         &
-                        )
-                      qscatIRgrid(j,1,m,iaer) =                 &
-                        qscatIRgrid(j,1,m,iaer) +               &
-                        weightgaus(gausind) *                   &
-                        (                                       &
-                        omegIRb(m,gausind,iaer) *               &
-                        qsqrefIRb(m,gausind,iaer) *             &
-                        qrefVISb(gausind,iaer) *                &
-                        pi*radGAUSb(gausind,iaer,idomain) *     &
-                        radGAUSb(gausind,iaer,idomain) *        &
-                        distb(j,1,iaer,idomain,gausind) +       &
-                        omegIRa(m,gausind,iaer) *               &
-                        qsqrefIRa(m,gausind,iaer) *             &
-                        qrefVISa(gausind,iaer) *                &
-                        pi*radGAUSa(gausind,iaer,idomain) *     &
-                        radGAUSa(gausind,iaer,idomain) *        &
-                        dista(j,1,iaer,idomain,gausind)         &
-                        )
-                      gIRgrid(j,1,m,iaer) =                     &
-                        gIRgrid(j,1,m,iaer) +                   &
-                        weightgaus(gausind) *                   &
-                        (                                       &
-                        omegIRb(m,gausind,iaer) *               &
-                        qsqrefIRb(m,gausind,iaer) *             &
-                        qrefVISb(gausind,iaer) *                &
-                        gIRb(m,gausind,iaer) *                  &
-                        pi*radGAUSb(gausind,iaer,idomain) *     &
-                        radGAUSb(gausind,iaer,idomain) *        &
-                        distb(j,1,iaer,idomain,gausind) +       &
-                        omegIRa(m,gausind,iaer) *               &
-                        qsqrefIRa(m,gausind,iaer) *             &
-                        qrefVISa(gausind,iaer) *                &
-                        gIRa(m,gausind,iaer) *                  &
-                        pi*radGAUSa(gausind,iaer,idomain) *     &
-                        radGAUSa(gausind,iaer,idomain) *        &
-                        dista(j,1,iaer,idomain,gausind)         &
-                        )
-                    ENDDO
-                    qrefIRgrid(j,1,iaer) =                      &
-                      qrefIRgrid(j,1,iaer) +                    &
-                      weightgaus(gausind) *                     &
-                      (                                         &
-                      qrefIRb(gausind,iaer) *                   &
-                      pi*radGAUSb(gausind,iaer,idomain) *       &
-                      radGAUSb(gausind,iaer,idomain) *          &
-                      distb(j,1,iaer,idomain,gausind) +         &
-                      qrefIRa(gausind,iaer) *                   &
-                      pi*radGAUSa(gausind,iaer,idomain) *       &
-                      radGAUSa(gausind,iaer,idomain) *          &
-                      dista(j,1,iaer,idomain,gausind)           &
-                      )
-                    qscatrefIRgrid(j,1,iaer) =                  &
-                      qscatrefIRgrid(j,1,iaer) +                &
-                      weightgaus(gausind) *                     &
-                      (                                         &
-                      omegrefIRb(gausind,iaer) *                &
-                      qrefIRb(gausind,iaer) *                   &
-                      pi*radGAUSb(gausind,iaer,idomain) *       &
-                      radGAUSb(gausind,iaer,idomain) *          &
-                      distb(j,1,iaer,idomain,gausind) +         &
-                      omegrefIRa(gausind,iaer) *                &
-                      qrefIRa(gausind,iaer) *                   &
-                      pi*radGAUSa(gausind,iaer,idomain) *       &
-                      radGAUSa(gausind,iaer,idomain) *          &
-                      dista(j,1,iaer,idomain,gausind)           &
-                      )
-                  ENDDO
- 
-                  qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /          &
-                                normd(j,1,iaer,idomain)
-                  qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /  &
-                                normd(j,1,iaer,idomain)
-                  omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /   &
-                               qrefIRgrid(j,1,iaer)
-                  DO m=1,L_NSPECTI
-                    qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /    &
-                                normd(j,1,iaer,idomain)
-                    qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /  &
-                                normd(j,1,iaer,idomain)
-                    gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /          &
-                                qscatIRgrid(j,1,m,iaer) /              &
-                                normd(j,1,iaer,idomain)
-
-                    qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /  &
-                                qrefVISgrid(j,1,iaer)
-                    omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /   &
-                                qextIRgrid(j,1,m,iaer)
-                  ENDDO
-                ENDIF                  ! --------------------------
-                checkgrid(j,1,iaer,idomain) = .true.
-              ENDIF !checkgrid
-          ENDDO !grid_i
-!         2.4 Linear interpolation
-          k1 = (1-kx)
-          k2 = kx
-          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
-          DO m=1,L_NSPECTV
-             QVISsQREF3d(ig,lg,m,iaer) =                           &
-                        k1*qsqrefVISgrid(grid_i,1,m,iaer) +        &
-                        k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
-            omegaVIS3d(ig,lg,m,iaer) =                             &
-                        k1*omegVISgrid(grid_i,1,m,iaer) +          &
-                        k2*omegVISgrid(grid_i+1,1,m,iaer)
-            gVIS3d(ig,lg,m,iaer) =                                 &
-                        k1*gVISgrid(grid_i,1,m,iaer) +             &
-                        k2*gVISgrid(grid_i+1,1,m,iaer)
-          ENDDO !L_NSPECTV
-          QREFvis3d(ig,lg,iaer) =                                  &
-                        k1*qrefVISgrid(grid_i,1,iaer) +            &
-                        k2*qrefVISgrid(grid_i+1,1,iaer)
-          omegaREFvis3d(ig,lg,iaer) =                              &
-                        k1*omegrefVISgrid(grid_i,1,iaer) +         &
-                        k2*omegrefVISgrid(grid_i+1,1,iaer)
-          ELSE                   ! INFRARED -----------------------
-          DO m=1,L_NSPECTI
-            QIRsQREF3d(ig,lg,m,iaer) =                             &
-                        k1*qsqrefIRgrid(grid_i,1,m,iaer) +         &
-                        k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
-            omegaIR3d(ig,lg,m,iaer) =                              &
-                        k1*omegIRgrid(grid_i,1,m,iaer) +           &
-                        k2*omegIRgrid(grid_i+1,1,m,iaer) 
-            gIR3d(ig,lg,m,iaer) =                                  & 
-                        k1*gIRgrid(grid_i,1,m,iaer) +              &
-                        k2*gIRgrid(grid_i+1,1,m,iaer)
-          ENDDO !L_NSPECTI
-          QREFir3d(ig,lg,iaer) =                                   &
-                        k1*qrefIRgrid(grid_i,1,iaer) +             &
-                        k2*qrefIRgrid(grid_i+1,1,iaer)
-          omegaREFir3d(ig,lg,iaer) =                               &
-                        k1*omegrefIRgrid(grid_i,1,iaer) +          &
-                        k2*omegrefIRgrid(grid_i+1,1,iaer)
-          ENDIF                  ! --------------------------------
-        ENDDO !nlayer
-      ENDDO !ngrid
-
-!==================================================================
-
-
-
-      ENDDO ! idomain
-
-      ENDIF ! nsize = 1
-
-      ENDDO ! iaer (loop on aerosol kind)
-
-      RETURN
-    END SUBROUTINE aeroptproperties
-
-
-
-      
Index: trunk/LMDZ.TITAN/libf/phytitan/aerosol_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/aerosol_mod.F90	(revision 1787)
+++ 	(revision )
@@ -1,18 +1,0 @@
-!==================================================================
-module aerosol_mod
-implicit none
-save
-!==================================================================
-
-!  aerosol indexes: these are initialized to be 0 if the
-!                 corresponding aerosol was not activated in callphys.def
-!                 -- otherwise a value is given in iniaerosol
-      logical :: noaero = .false.
-
-! two-layer simple aerosol model
-      integer :: iaero_back2lay = 0
-!$OMP THREADPRIVATE(noaero,iaero_back2lay)
-      
-!==================================================================
-end module aerosol_mod
-!==================================================================
Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1788)
@@ -1,10 +1,10 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    & 
-          tsurf,fract,dist_star,aerosol,                       &
+          tsurf,fract,dist_star,                               &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
           fluxabs_sw,fluxtop_dn,                               &
           OLR_nu,OSR_nu,                                       &
-          tau_col,firstcall,lastcall)
+          firstcall,lastcall)
 
       use mod_phys_lmdz_para, only : is_master
@@ -14,6 +14,4 @@
       use ioipsl_getin_p_mod, only: getin_p
       use gases_h
-      use radii_mod, only : su_aer_radii,back2lay_reffrad
-      use aerosol_mod, only : iaero_back2lay
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
@@ -64,5 +62,4 @@
       
       ! OUTPUT
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
       REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)             ! Heating rate (K/s) due to LW radiation.
       REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)             ! Heating rate (K/s) due to SW radiation.
@@ -75,23 +72,7 @@
       REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)        ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
       REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
-      REAL,INTENT(OUT) :: tau_col(ngrid)                 ! Diagnostic from aeropacity.
       REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
       
       
-      ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.   
-      REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
-      REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
-      REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
-      REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
-      REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
-      REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
-
-!      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
-!      REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
-
-      REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:)  ! aerosol effective radius (m)
-      REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
-!$OMP THREADPRIVATE(reffrad,nueffrad)
-
 !-----------------------------------------------------------------------
 !     Declaration of the variables required by correlated-k subroutines
@@ -114,5 +95,4 @@
       REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
 
-      REAL*8 tauaero(L_LEVELS,naerkind)
       REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
       REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
@@ -125,5 +105,5 @@
       REAL*8 albv(L_NSPECTV)                         ! Spectral Visible Albedo.
 
-      INTEGER ig,l,k,nw,iaer
+      INTEGER ig,l,k,nw
 
       real szangle
@@ -134,19 +114,6 @@
       real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
 
-      ! Local aerosol optical properties for each column on RADIATIVE grid.
-      real*8,save,allocatable ::  QXVAER(:,:,:)
-      real*8,save,allocatable ::  QSVAER(:,:,:)
-      real*8,save,allocatable ::  GVAER(:,:,:)
-      real*8,save,allocatable ::  QXIAER(:,:,:)
-      real*8,save,allocatable ::  QSIAER(:,:,:)
-      real*8,save,allocatable ::  GIAER(:,:,:)
-
-      real, dimension(:,:,:), save, allocatable :: QREFvis3d
-      real, dimension(:,:,:), save, allocatable :: QREFir3d
-!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d)
-
 
       ! Miscellaneous :
-      real*8  temp,temp1,temp2,pweight
       character(len=10) :: tmp1
       character(len=10) :: tmp2
@@ -167,28 +134,6 @@
       if(firstcall) then
 
-        ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
-        if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind))
-        if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind))
-        if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind))
-        if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind))
-        if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind))
-        if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind))
-
-         !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
-         IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
-         IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
-         ! Effective radius and variance of the aerosols
-         IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
-         IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind))
-
          call system('rm -f surf_vals_long.out')
 
-         if(naerkind.gt.4)then
-            print*,'Code not general enough to deal with naerkind > 4 yet.'
-            call abort
-         endif
-         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-         
-         
 !--------------------------------------------------
 !             Set up correlated k
@@ -207,6 +152,4 @@
          call setspv            ! Basic visible properties.
          call sugas_corrk       ! Set up gaseous absorption properties.
-         call suaer_corrk       ! Set up aerosol optical properties.
-        
 
          OLR_nu(:,:) = 0.
@@ -240,24 +183,4 @@
 !=======================================================================
  
-      qxvaer(:,:,:)=0.0
-      qsvaer(:,:,:)=0.0
-      gvaer(:,:,:) =0.0
-
-      qxiaer(:,:,:)=0.0
-      qsiaer(:,:,:)=0.0
-      giaer(:,:,:) =0.0
-
-!--------------------------------------------------
-!     Effective radius and variance of the aerosols
-!--------------------------------------------------
-
-      do iaer=1,naerkind
-     
-          if(iaer.eq.iaero_back2lay)then
-            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
-         endif
-         
-     end do !iaer=1,naerkind.
-
 
       ! How much light do we get ?
@@ -266,16 +189,4 @@
       end do
 
-      ! Get 3D aerosol optical properties.
-      call aeroptproperties(ngrid,nlayer,reffrad,nueffrad,         &
-           QVISsQREF3d,omegaVIS3d,gVIS3d,                          &
-           QIRsQREF3d,omegaIR3d,gIR3d,                             &
-           QREFvis3d,QREFir3d)                                     
-
-      ! Get aerosol optical depths.
-      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
-           reffrad,QREFvis3d,QREFir3d,                             & 
-           tau_col)                
-          
-
 
 !-----------------------------------------------------------------------    
@@ -288,139 +199,4 @@
 !=======================================================================
 
-
-!-----------------------------------------------------------------------
-!    Aerosol optical properties Qext, Qscat and g.
-!    The transformation in the vertical is the same as for temperature.
-!-----------------------------------------------------------------------
-           
-           
-            do iaer=1,naerkind
-               ! Shortwave.
-               do nw=1,L_NSPECTV 
-               
-                  do l=1,nlayer
-
-                     temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer)         &
-                         *QREFvis3d(ig,nlayer+1-l,iaer)
-
-                     temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
-                         *QREFvis3d(ig,max(nlayer-l,1),iaer)
-
-                     qxvaer(2*l,nw,iaer)  = temp1
-                     qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
-
-                     temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
-                     temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
-
-                     qsvaer(2*l,nw,iaer)  = temp1
-                     qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
-
-                     temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
-                     temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
-
-                     gvaer(2*l,nw,iaer)  = temp1
-                     gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
-
-                  end do ! nlayer
-
-                  qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
-                  qxvaer(2*nlayer+1,nw,iaer)=0.
-
-                  qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
-                  qsvaer(2*nlayer+1,nw,iaer)=0.
-
-                  gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
-                  gvaer(2*nlayer+1,nw,iaer)=0.
-
-               end do ! L_NSPECTV
-             
-               do nw=1,L_NSPECTI
-                  ! Longwave
-                  do l=1,nlayer
-
-                     temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer)         &
-                          *QREFir3d(ig,nlayer+1-l,iaer)
-
-                     temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer)    &
-                          *QREFir3d(ig,max(nlayer-l,1),iaer)
-
-                     qxiaer(2*l,nw,iaer)  = temp1
-                     qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
-
-                     temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
-                     temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
-
-                     qsiaer(2*l,nw,iaer)  = temp1
-                     qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
-
-                     temp1=gir3d(ig,nlayer+1-l,nw,iaer)
-                     temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
-
-                     giaer(2*l,nw,iaer)  = temp1
-                     giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
-
-                  end do ! nlayer
-
-                  qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
-                  qxiaer(2*nlayer+1,nw,iaer)=0.
-
-                  qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
-                  qsiaer(2*nlayer+1,nw,iaer)=0.
-
-                  giaer(1,nw,iaer)=giaer(2,nw,iaer)
-                  giaer(2*nlayer+1,nw,iaer)=0.
-
-               end do ! L_NSPECTI
-               
-            end do ! naerkind
-
-            ! Test / Correct for freaky s. s. albedo values.
-            do iaer=1,naerkind
-               do k=1,L_LEVELS
-
-                  do nw=1,L_NSPECTV
-                     if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
-                        print*,'Serious problems with qsvaer values' 
-                        print*,'in callcorrk'
-                        call abort
-                     endif
-                     if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
-                        qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
-                     endif
-                  end do
-
-                  do nw=1,L_NSPECTI 
-                     if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
-                        print*,'Serious problems with qsiaer values'
-                        print*,'in callcorrk'
-                        call abort
-                     endif
-                     if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
-                        qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
-                     endif
-                  end do
-
-               end do ! L_LEVELS
-            end do ! naerkind
-
-!-----------------------------------------------------------------------
-!     Aerosol optical depths
-!-----------------------------------------------------------------------
-            
-         do iaer=1,naerkind     ! a bug was here           
-            do k=0,nlayer-1
-               
-               pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
-                       (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
-               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
-               tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
-               tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
-
-            end do
-            ! boundary conditions
-            tauaero(1,iaer)          = tauaero(2,iaer)
-            !tauaero(1,iaer)          = 0.
-            
-         end do ! naerkind
 
          ! Albedo and Emissivity.
@@ -586,6 +362,5 @@
             
             call optcv(dtauv,tauv,taucumv,plevrad,                 &
-                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
-                 tmid,pmid,taugsurf)
+                 wbarv,cosbv,tauray,tmid,pmid,taugsurf)
 
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
@@ -629,6 +404,5 @@
 
          call optci(plevrad,tlevrad,dtaui,taucumi,                  &
-              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
-              taugsurfi)
+              cosbi,wbari,tmid,pmid,taugsurfi)
 
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
@@ -655,5 +429,4 @@
             print*,'fluxtop_dn=',fluxtop_dn(ig)
             print*,'acosz=',acosz
-            print*,'aerosol=',aerosol(ig,:,:)
             print*,'temp=   ',pt(ig,:)
             print*,'pplay=  ',pplay(ig,:)
@@ -754,6 +527,4 @@
 !$OMP END MASTER
 !$OMP BARRIER
-        IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
-        IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
       endif
 
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1788)
@@ -31,8 +31,4 @@
       logical,save :: mass_redistrib
 !$OMP THREADPRIVATE(stelbbody,ozone,tracer,mass_redistrib)
-      logical,save :: sedimentation
-!$OMP THREADPRIVATE(sedimentation)
-      logical,save :: aeroback2lay
-!$OMP THREADPRIVATE(aeroback2lay)
       logical,save :: nosurf
       logical,save :: oblate
@@ -41,24 +37,11 @@
       integer,save :: ichim
       integer,save :: iddist
-      integer,save :: iaervar
       integer,save :: iradia
       integer,save :: startype
-!$OMP THREADPRIVATE(ichim,iddist,iaervar,iradia,startype)
+!$OMP THREADPRIVATE(ichim,iddist,iradia,startype)
 
       real,save :: Fat1AU
       real,save :: stelTbb
 !$OMP THREADPRIVATE(Fat1AU,stelTbb)
-      real,save :: tplanet
-      real,save :: obs_tau_col_tropo
-      real,save :: obs_tau_col_strato
-!$OMP THREADPRIVATE(tplanet,obs_tau_col_tropo,obs_tau_col_strato)
-      real,save :: pres_bottom_tropo
-      real,save :: pres_top_tropo
-      real,save :: pres_bottom_strato
-      real,save :: pres_top_strato
-!$OMP THREADPRIVATE(pres_bottom_tropo,pres_top_tropo,pres_bottom_strato,pres_top_strato)
-      real,save :: size_tropo
-      real,save :: size_strato
-!$OMP THREADPRIVATE(size_tropo,size_strato)
       real,save :: pceil
 !$OMP THREADPRIVATE(pceil)
Index: trunk/LMDZ.TITAN/libf/phytitan/callsedim.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callsedim.F	(revision 1787)
+++ 	(revision )
@@ -1,128 +1,0 @@
-      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
-     &                pplev,zlev, pt, pdt,
-     &                pq, pdqfi, pdqsed,pdqs_sed,nq)
-
-      use radinc_h, only : naerkind
-      USE tracer_h, only : radius, rho_q
-      use comcstfi_mod, only: g
-
-      IMPLICIT NONE
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Calculates sedimentation of aerosols depending on their
-!     density and radius.
-!     
-!     Authors
-!     -------
-!     F. Forget (1999)
-!     Tracer generalisation by E. Millour (2009)
-!     
-!==================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-c   arguments:
-c   ----------
-
-      integer,intent(in):: ngrid ! number of horizontal grid points
-      integer,intent(in):: nlay  ! number of atmospheric layers
-      real,intent(in):: ptimestep  ! physics time step (s)
-      real,intent(in):: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
-      real,intent(in):: pt(ngrid,nlay)      ! temperature at mid-layer (K)
-      real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature
-      real,intent(in):: zlev(ngrid,nlay+1)  ! altitude at layer boundaries
-      integer,intent(in) :: nq ! number of tracers
-      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
-      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency on tracers before
-                                               ! sedimentation (kg/kg.s-1)
-      
-      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
-      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
-      
-c   local:
-c   ------
-
-      INTEGER l,ig, iq
-
-      ! for particles with varying radii:
-      real reffrad(ngrid,nlay,naerkind) ! particle radius (m) 
-      real nueffrad(ngrid,nlay,naerkind) ! aerosol effective radius variance
-
-      real zqi(ngrid,nlay,nq) ! to locally store tracers
-      real zt(ngrid,nlay) ! to locally store temperature (K)
-      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
-      real epaisseur (ngrid,nlay) ! Layer thickness (m)
-      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
-
-
-      LOGICAL,SAVE :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-
-c    ** un petit test de coherence
-c       --------------------------
-
-      IF (firstcall) THEN
-        firstcall=.false.
-      ENDIF ! of IF (firstcall)
-      
-!=======================================================================
-!     Preliminary calculation of the layer characteristics
-!     (mass (kg.m-2), thickness (m), etc.
-
-      do  l=1,nlay
-        do ig=1, ngrid
-          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
-          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
-          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
-        end do
-      end do
-
-!======================================================================
-! Calculate the transport due to sedimentation for each tracer
-! [This has been rearranged by L. Kerber to allow the sedimentation
-!  of general tracers.]
- 
-      do iq=1,nq
-       if(radius(iq).gt.1.e-9) then   
-!         (no sedim for gases)      
-
-! store locally updated tracers
-
-          do l=1,nlay 
-            do ig=1, ngrid
-              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
-            enddo
-          enddo ! of do l=1,nlay
-        
-!======================================================================
-! Sedimentation 
-!======================================================================
-! General Case
-             call newsedim(ngrid,nlay,1,ptimestep,
-     &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
-     &            zqi(1,1,iq),wq)
-
-!=======================================================================
-!     Calculate the tendencies
-!======================================================================
-
-          do ig=1,ngrid
-!     Ehouarn: with new way of tracking tracers by name, this is simply
-            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
-          end do
-
-          DO l = 1, nlay
-            DO ig=1,ngrid
-              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
-     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
-            ENDDO
-          ENDDO
-       endif ! of no gases
-      enddo ! of do iq=1,nq
-      return
-      end
Index: trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 1788)
@@ -18,7 +18,4 @@
       character(len=12),parameter :: surfdir="surface_data"
       
-      ! aerdir stores aerosol properties files (optprop_*dat files)
-      character(LEN=18),parameter :: aerdir="aerosol_properties" 
-
       end module datafile_mod
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.TITAN/libf/phytitan/iniaerosol.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/iniaerosol.F	(revision 1787)
+++ 	(revision )
@@ -1,49 +1,0 @@
-      SUBROUTINE iniaerosol()
-
-
-      use radinc_h, only: naerkind
-      use aerosol_mod
-      use callkeys_mod, only: aeroback2lay
-
-      IMPLICIT NONE
-c=======================================================================
-c   subject:
-c   --------
-c   Initialization related to aerosols 
-c   (Chemical species, ice...)   
-c
-c   author: Laura Kerber, S. Guerlet
-c   ------
-c        
-c=======================================================================
-
-      integer ia
-
-      ia=0
-       
-      if (aeroback2lay) then
-         ia=ia+1
-         iaero_back2lay=ia
-      endif
-      write(*,*) '--- Two-layer aerosol = ', iaero_back2lay
-
-      write(*,*) '=== Number of aerosols= ', ia
-      
-! For the zero aerosol case, we currently make a dummy co2 aerosol which is zero everywhere.
-! (See aeropacity.F90 for how this works). A better solution would be to turn off the 
-! aerosol machinery in the no aerosol case, but this would be complicated. LK
-
-      if (ia.eq.0) then  !For the zero aerosol case. 
-         noaero = .true.
-         ia = 1
-      endif
-
-      if (ia.ne.naerkind) then
-          print*, 'Aerosols counted not equal to naerkind'
-          print*, 'Compile with tag -s',ia,'to run'
-          print*, 'or change options in callphys.def'
-          print*, 'Abort in iniaerosol.F'
-          call abort
-      endif
-
-      end
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1788)
@@ -9,5 +9,5 @@
              prad,pg,pr,pcpp)
 
-  use radinc_h, only: ini_radinc_h, naerkind
+  use radinc_h, only: ini_radinc_h
   use datafile_mod, only: datadir
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
@@ -385,9 +385,4 @@
      endif
 
-     write(*,*)"Default planetary temperature?"
-     tplanet=215.0
-     call getin_p("tplanet",tplanet)
-     write(*,*)" tplanet = ",tplanet
-
      write(*,*)"Which star?"
      startype=1 ! default value = Sol
@@ -399,71 +394,4 @@
      call getin_p("Fat1AU",Fat1AU)
      write(*,*)" Fat1AU = ",Fat1AU
-
-
-! TRACERS:
-
-!         write(*,*)"Number of radiatively active aerosols:"
-!         naerkind=0. ! default value
-!         call getin_p("naerkind",naerkind)
-!         write(*,*)" naerkind = ",naerkind
-
- 
-!=================================
-
-     write(*,*)"Radiatively active two-layer aersols?"
-     aeroback2lay=.false.     ! default value
-     call getin_p("aeroback2lay",aeroback2lay)
-     write(*,*)" aeroback2lay = ",aeroback2lay
-
-     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
-                    "in the tropospheric layer (visible)"
-     obs_tau_col_tropo=8.D0
-     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
-     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
-
-     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
-                    "in the stratospheric layer (visible)"
-     obs_tau_col_strato=0.08D0
-     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
-     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
-
-     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
-     pres_bottom_tropo=66000.0
-     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
-     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
-
-     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
-     pres_top_tropo=18000.0
-     call getin_p("pres_top_tropo",pres_top_tropo)
-     write(*,*)" pres_top_tropo = ",pres_top_tropo
-
-     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
-     pres_bottom_strato=2000.0
-     call getin_p("pres_bottom_strato",pres_bottom_strato)
-     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
-
-     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
-     pres_top_strato=100.0
-     call getin_p("pres_top_strato",pres_top_strato)
-     write(*,*)" pres_top_strato = ",pres_top_strato
-
-     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
-                    "tropospheric layer, in meters"
-     size_tropo=2.e-6
-     call getin_p("size_tropo",size_tropo)
-     write(*,*)" size_tropo = ",size_tropo
-
-     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
-                    "stratospheric layer, in meters"
-     size_strato=1.e-7
-     call getin_p("size_strato",size_strato)
-     write(*,*)" size_strato = ",size_strato
-
-!=================================
-
-     write(*,*) "Gravitationnal sedimentation ?"
-     sedimentation=.false. ! default value
-     call getin_p("sedimentation",sedimentation)
-     write(*,*) " sedimentation = ",sedimentation       
 
      write(*,*) "Does user want to force cpp and mugaz?"
Index: trunk/LMDZ.TITAN/libf/phytitan/initracer.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/initracer.F	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/initracer.F	(revision 1788)
@@ -32,12 +32,7 @@
 
 
-c-----------------------------------------------------------------------
-c  radius(nq)      ! aerosol particle radius (m)
+c------------------------------------------------
 c  rho_q(nq)       ! tracer densities (kg.m-3)
-c  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
-c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
-c  alpha_devil(nq) ! lifting coeeficient by dust devil
-c-----------------------------------------------------------------------
-
+c------------------------------------------------
        nqtot=nq
        !! we allocate once for all arrays in common in tracer_h.F90
@@ -45,18 +40,5 @@
        IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
        ALLOCATE(mmol(nq))
-       ALLOCATE(radius(nq))
        ALLOCATE(rho_q(nq))
-       ALLOCATE(qext(nq))
-       ALLOCATE(alpha_lift(nq))
-       ALLOCATE(alpha_devil(nq))
-       ALLOCATE(qextrhor(nq))
-       !! initialization
-       alpha_lift(:)=0.
-       alpha_devil(:)=0.
-       
-       ! Added by JVO 2017 : these arrays are handled later
-       ! -> initialization is the least we can do, please !!!
-       radius(:)=0.
-       qext(:)=0.
 
 ! Initialization: get tracer names from the dynamics and check if we are
@@ -373,14 +355,3 @@
 
 
-c     Output for records:
-c     ~~~~~~~~~~~~~~~~~~
-      write(*,*)
-      Write(*,*) '******** initracer : dust transport parameters :'
-      write(*,*) 'alpha_lift = ', alpha_lift
-      write(*,*) 'alpha_devil = ', alpha_devil
-      write(*,*) 'radius  = ', radius
-      
-      write(*,*) 'Qext  = ', qext 
-      write(*,*)
-
       end
Index: trunk/LMDZ.TITAN/libf/phytitan/newsedim.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/newsedim.F	(revision 1787)
+++ 	(revision )
@@ -1,197 +1,0 @@
-      SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep,
-     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq)
-      
-      use comcstfi_mod, only: r, g
-      IMPLICIT NONE
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!      Calculates sedimentation of 1 tracer 
-!      of radius rd (m) and density rho (kg.m-3) 
-!     
-!==================================================================
-
-!-----------------------------------------------------------------------
-!   declarations
-!   ------------
-
-!   arguments
-!   ---------
-
-      integer,intent(in) :: ngrid  ! number of atmospheric columns
-      integer,intent(in) :: nlay  ! number of atmospheric layers
-      integer,intent(in) :: naersize   ! number of particle sizes (1 or number
-                                       ! of grid boxes)
-      real,intent(in) :: ptimestep ! physics time step (s)
-      real,intent(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressures (Pa)
-      real,intent(in) :: pt(ngrid,nlay)    ! mid-layer temperatures (K)
-      real,intent(in) :: masse (ngrid,nlay) ! atmospheric mass (kg)
-      real,intent(in) :: epaisseur (ngrid,nlay)  ! thickness of atm. layers (m)
-      real,intent(in) :: rd(naersize) ! particle radius (m)
-      real,intent(in) :: rho ! particle density (kg.m-3)
-      real,intent(inout) :: pqi(ngrid,nlay)  ! tracer   (e.g. ?/kg)
-      real,intent(out) :: wq(ngrid,nlay+1)  ! flux of tracer during timestep (?/m-2)
-      
-c   local:
-c   ------
-
-      INTEGER l,ig, k, i
-      REAL rfall
-
-      LOGICAL,SAVE :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-
-c    Traceurs :
-c    ~~~~~~~~ 
-      real traversee (ngrid,nlay)
-      real vstokes(ngrid,nlay)
-      real w(ngrid,nlay)
-      real ptop, dztop, Ep, Stra
-
-
-c    Physical constant
-c    ~~~~~~~~~~~~~~~~~
-c     Gas molecular viscosity (N.s.m-2)
-      real,parameter :: visc=1.e-5       ! CO2
-c     Effective gas molecular radius (m)
-      real,parameter :: molrad=2.2e-10   ! CO2
-
-c     local and saved variable
-      real,save :: a,b
-!$OMP THREADPRIVATE(a,b)
-
-c    ** un petit test de coherence
-c       --------------------------
-
-
-      !print*,'temporary fixed particle rad in newsedim!!'
-
-      IF (firstcall) THEN
-         firstcall=.false.
-
-
-
-!=======================================================================
-!     Preliminary calculations for sedimenation velocity
-
-!       - Constant to compute stokes speed simple formulae
-!        (Vstokes =  b * rho* r**2   avec   b= (2/9) * rho * g / visc
-         b = 2./9. * g / visc
-
-!       - Constant  to compute gas mean free path
-!        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
-         a = 0.707*8.31/(4*3.1416* molrad**2  * 6.023e23)
-
-c       - Correction to account for non-spherical shape (Murphy et al.  1990)
-c   (correction = 0.85 for irregular particles, 0.5 for disk shaped particles)
-c        a = a    * 0.85
-
-
-      ENDIF
-
-c-----------------------------------------------------------------------
-c    1. initialisation
-c    -----------------
-
-c     Sedimentation velocity (m/s)
-c     ~~~~~~~~~~~~~~~~~~~~~~
-c     (stokes law corrected for low pressure by the Cunningham
-c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
-       
-        do  l=1,nlay
-          do ig=1, ngrid
-            if (naersize.eq.1) then 
-              rfall=rd(1)
-            else
-              i=ngrid*(l-1)+ig
-              rfall=rd(i) ! how can this be correct!!?
-            endif  
-
-            vstokes(ig,l) = b * rho * rfall**2 *
-     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
-
-c           Layer crossing time (s) :
-            traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l)
-          end do
-        end do
-
-
-c     Calcul de la masse d'atmosphere correspondant a q transferee
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c     (e.g. on recherche le niveau  en dessous de laquelle le traceur
-c      va traverser le niveau intercouche l : "dztop" est sa hauteur
-c      au dessus de l (m), "ptop" est sa pression (Pa))
-      do  l=1,nlay
-        do ig=1, ngrid
-             
-             dztop = vstokes(ig,l)*  ptimestep 
-             Ep=0
-             k=0
-           w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS)
-c **************************************************************
-c            Simple Method
-cc             w(ig,l) =
-cc     &       (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
-cc             write(*,*) 'OK simple method l,w =', l, w(ig,l)
-cc            write(*,*) 'OK simple method dztop =', dztop
-           w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
-             !!! Diagnostic: JF. Fix: AS. Date: 05/11
-             !!! Probleme arrondi avec la quantite ci-dessus
-             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
-             !!! ---> dans ce cas on utilise le developpement limite !
-             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2  
-
-             IF ( w(ig,l) .eq. 0. ) THEN
-                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
-             ELSE
-                w(ig,l) = w(ig,l) * pplev(ig,l) / g
-             ENDIF
-! LK borrowed simple method from Mars model (AS/JF)
-
-!**************************************************************
-cccc         Complex method :
-           if (dztop.gt.epaisseur(ig,l)) then
-cccc            Cas ou on "epuise" la couche l : On calcule le flux
-cccc            Venant de dessus en tenant compte de la variation de Vstokes
-
-               Ep= epaisseur(ig,l)
-               Stra= traversee(ig,l)
-               do while(dztop.gt.Ep.and.l+k+1.le.nlay)
-                 k=k+1
-                 dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra)
-                 Ep = Ep + epaisseur(ig,l+k)
-                 Stra = Stra + traversee(ig,l+k)
-               enddo 
-               Ep = Ep - epaisseur(ig,l+k)
-!             ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
-             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
-	     IF ( ptop .eq. 1. ) THEN
-                !PRINT*, 'newsedim: exposant trop petit ', ig, l
-                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
-             ELSE
-                ptop=pplev(ig,l+k) * ptop
-             ENDIF
-
-             w(ig,l) = (pplev(ig,l) - ptop)/g
-
-            endif   !!! complex method
-c
-cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
-cc           write(*,*) 'OK new    method dztop =', dztop
-cc       if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop
-cc       if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep
-cc            if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l)
-cc            if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l)
-cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
-c **************************************************************
-
-        end do
-      end do
-
-      call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)
-c         write(*,*) ' newsed: wq(6), wq(7), q(6)',
-c    &                wq(1,6),wq(1,7),pqi(1,6)
-
-      END
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1788)
@@ -1,5 +1,4 @@
 subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
-     QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
-     TMID,PMID,TAUGSURF)
+     COSBI,WBARI,TMID,PMID,TAUGSURF)
 
   use radinc_h
@@ -41,12 +40,4 @@
   real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
 
-  ! for aerosols
-  real*8  QXIAER(L_LEVELS,L_NSPECTI,NAERKIND)
-  real*8  QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
-  real*8  GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
-  real*8  TAUAERO(L_LEVELS,NAERKIND)
-  real*8  TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND)
-  real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
-
   ! Titan customisation
   ! J. Vatant d'Ollone (2016)
@@ -68,5 +59,5 @@
 
   real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
-  real*8 DCONT,DAERO
+  real*8 DCONT
   double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
   double precision p_cross
@@ -117,13 +108,4 @@
   end do                    ! levels
 
-  ! Spectral dependance of aerosol absorption
-  do iaer=1,naerkind
-     DO NW=1,L_NSPECTI
-        do K=2,L_LEVELS
-           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
-        end do                    ! levels
-     END DO
-  end do
-
   do NW=1,L_NSPECTI
 
@@ -132,6 +114,4 @@
         ilay = k / 2 ! int. arithmetic => gives the gcm layer index
         
-        DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
-
         !================= Titan customisation ========================================
         call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
@@ -240,5 +220,4 @@
            DTAUKI(K,nw,ng) = TAUGAS    & 
                              + DCONT   & ! For parameterized continuum absorption
-			     + DAERO   & ! For aerosol absorption
 			     + DHAZE_T(K,NW)  ! For Titan haze
 
@@ -251,5 +230,4 @@
         DTAUKI(K,nw,ng) = 0.d0      & 
                           + DCONT   & ! For parameterized continuum absorption
-	                  + DAERO   & ! For aerosol absorption
 	                  + DHAZE_T(K,NW)     ! For Titan Haze
 
@@ -262,12 +240,4 @@
   ! ======================================================================
 
-  do iaer=1,naerkind
-    DO NW=1,L_NSPECTI
-     DO K=2,L_LEVELS
-           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
-     ENDDO
-    ENDDO
-  end do
-  
   ! Haze scattering
   DO NW=1,L_NSPECTI
@@ -280,6 +250,5 @@
      DO L=1,L_NLAYRAD-1
         K              = 2*L+1
-        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
-                      + DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
+        btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
      END DO ! L vertical loop
      
@@ -287,5 +256,5 @@
      L           = L_NLAYRAD
      K           = 2*L+1
-     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW)
+     btemp(L,NW) = DHAZES_T(K,NW)
      
   END DO                    ! NW spectral loop
@@ -301,9 +270,4 @@
         atemp = 0.
         if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
-           do iaer=1,naerkind
-              atemp = atemp +                                     &
-                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
-                   GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
-           end do
            atemp = atemp +                   &
                 ASF_T(K,NW)*DHAZES_T(K,NW) + &
@@ -332,7 +296,4 @@
      atemp = 0.
      if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
-        do iaer=1,naerkind
-           atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
-        end do
         atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
         WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1788)
@@ -1,5 +1,4 @@
 SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
-     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
-     TAURAY,TAUAERO,TMID,PMID,TAUGSURF)
+     WBARV,COSBV,TAURAY,TMID,PMID,TAUGSURF)
 
   use radinc_h
@@ -46,12 +45,4 @@
   real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
   real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-
-  ! for aerosols
-  real*8  QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
-  real*8  QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
-  real*8  GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
-  real*8  TAUAERO(L_LEVELS,NAERKIND)
-  real*8  TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND)
-  real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
 
   ! Titan customisation
@@ -76,5 +67,5 @@
 
   real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
-  real*8 DCONT,DAERO
+  real*8 DCONT
   real*8 DRAYAER
   double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
@@ -129,13 +120,4 @@
   end do                    ! levels
 
-  ! Spectral dependance of aerosol absorption
-  do iaer=1,naerkind
-     do NW=1,L_NSPECTV
-        do K=2,L_LEVELS
-           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
-        end do                    ! levels
-     end do
-  end do
-  
   ! Rayleigh scattering
   do NW=1,L_NSPECTV
@@ -158,8 +140,4 @@
         DRAYAER = TRAY(K,NW)
         !     DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
-        do iaer=1,naerkind
-           DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
-        end do
-
         DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol
 
@@ -275,12 +253,4 @@
   !     we need to calculate the scattering albedo and asymmetry factors
 
-  do iaer=1,naerkind
-    DO NW=1,L_NSPECTV
-      DO K=2,L_LEVELS
-           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
-      ENDDO
-    ENDDO
-  end do
-  
   ! Haze scattering
   DO NW=1,L_NSPECTV
@@ -294,8 +264,6 @@
      DO L=1,L_NLAYRAD-1
         K              = 2*L+1
-	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind)) &
-                    + ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
-        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
-                    + DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
+	atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
+        btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
 	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
 	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
@@ -306,8 +274,6 @@
      L           = L_NLAYRAD
      K           = 2*L+1
-     atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
-                 + ASF_T(K,NW)*DHAZES_T(K,NW)
-     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) &
-                 + DHAZES_T(K,NW)
+     atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW)
+     btemp(L,NW) = DHAZES_T(K,NW)
      ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
      btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1788)
@@ -14,5 +14,5 @@
                   pdu,pdv,pdt,pdq,pdpsrf)
  
-      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
+      use radinc_h, only : L_NSPECTI,L_NSPECTV
       use radcommon_h, only: sigma, glat, grav, BWNV
       use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
@@ -22,6 +22,5 @@
       use geometry_mod, only: latitude, longitude, cell_area
       USE comgeomfi_h, only: totarea, totarea_planet
-      USE tracer_h, only: noms, mmol, radius, qext, &
-                          alpha_lift, alpha_devil, qextrhor
+      USE tracer_h, only: noms, mmol
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyetat0_mod, only: phyetat0
@@ -74,5 +73,4 @@
 !      V. Tracers
 !         V.1. Chemistry
-!         V.2. Aerosols and particles.
 !         V.3. Updates (pressure variations, surface budget).
 !         V.4. Surface Tracer Update.
@@ -225,8 +223,4 @@
 ! -----------------
 
-!     Aerosol (dust or ice) extinction optical depth  at reference wavelength 
-!     for the "naerkind" optically active aerosols:
-
-      real aerosol(ngrid,nlayer,naerkind) ! Aerosols.
       real zh(ngrid,nlayer)               ! Potential temperature (K).
       real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
@@ -274,5 +268,4 @@
       real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
       real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
-      real zdqssed(ngrid,nq)                ! Callsedim routine.
       real zdqsurfmr(ngrid,nq)              ! Mass_redistribution routine.
                   
@@ -281,5 +274,4 @@
       real zdqdif(ngrid,nlayer,nq)    ! Turbdiff/vdifc routines.
       real zdqevap(ngrid,nlayer)      ! Turbdiff routine.
-      real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
       real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
       
@@ -327,11 +319,7 @@
       real vmr(ngrid,nlayer)                        ! volume mixing ratio
       real time_phys
-      real,allocatable,dimension(:),save :: tau_col ! Total Aerosol Optical Depth.
-!$OMP THREADPRIVATE(tau_col)
       
       real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
       
-      real qcol(ngrid,nq) ! Tracer Column Mass (kg/m2).
-
 !     to test energy conservation (RW+JL)
       real mass(ngrid,nlayer),massarea(ngrid,nlayer)
@@ -353,16 +341,8 @@
       real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid)                            ! For SW/LW flux.
       real albedo_equivalent1(ngrid)                                         ! For Equivalent albedo calculation.
-      real tau_col1(ngrid)                                                   ! For aerosol optical depth diagnostic.
-      real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV)                ! For Outgoing Radiation diagnostics.
       real tf, ntf    
 
       real,allocatable,dimension(:,:),save :: qsurf_hist
 !$OMP THREADPRIVATE(qsurf_hist)
-
-      real,allocatable,dimension(:,:,:),save :: reffrad  ! Aerosol effective radius (m). By RW
-      real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW
-!$OMP THREADPRIVATE(reffrad,nueffrad)
-
-      real reffcol(ngrid,naerkind)
 
 ! Local variables for Titan chemistry and microphysics (JVO 2017)
@@ -425,6 +405,4 @@
         ALLOCATE(zuprevious(ngrid,nlayer))
         ALLOCATE(qsurf_hist(ngrid,nq))
-        ALLOCATE(reffrad(ngrid,nlayer,naerkind))
-        ALLOCATE(nueffrad(ngrid,nlayer,naerkind))
         ALLOCATE(fluxsurf_lw(ngrid))
         ALLOCATE(fluxsurf_sw(ngrid))
@@ -439,5 +417,4 @@
         ALLOCATE(zdtlw(ngrid,nlayer))
         ALLOCATE(zdtsw(ngrid,nlayer))
-        ALLOCATE(tau_col(ngrid))
         ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
         ALLOCATE(qysat(nlayer,nq-nmicro))
@@ -455,12 +432,6 @@
          dtrad(:,:) = 0.0
          fluxrad(:) = 0.0
-         tau_col(:) = 0.0
          zdtsw(:,:) = 0.0
          zdtlw(:,:) = 0.0
-
-
-!        Initialize aerosol indexes.
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call iniaerosol()
 
       
@@ -781,9 +752,9 @@
                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
-                              tsurf,fract,dist_star,aerosol,                      &
+                              tsurf,fract,dist_star,                              &
                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
                               fluxsurfabs_sw,fluxtop_lw,                          &
                               fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
-                              tau_col,firstcall,lastcall)
+                              firstcall,lastcall)
 
                ! Radiative flux from the sky absorbed by the surface (W.m-2).
@@ -1057,25 +1028,4 @@
          endif
       
-  ! -------------------------
-  !   V.2. Aerosol particles
-  ! -------------------------
-
-         ! Sedimentation.
-         if (sedimentation) then
-        
-            zdqsed(1:ngrid,1:nlayer,1:nq) = 0.0
-            zdqssed(1:ngrid,1:nq)  = 0.0
-           
-            call callsedim(ngrid,nlayer,ptimestep,           &
-                          pplev,zzlev,pt,pdt,pq,pdq,        &
-                          zdqsed,zdqssed,nq)
-
-            ! Whether it falls as rain or snow depends only on the surface temperature
-            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqsed(1:ngrid,1:nlayer,1:nq)
-            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
-
-         end if ! end of 'sedimentation'
-
-
   ! ---------------
   !   V.3 Updates
@@ -1214,6 +1164,4 @@
                print*,'fluxtop_dn has gone crazy'
                print*,'fluxtop_dn=',fluxtop_dn(ig)
-               print*,'tau_col=',tau_col(ig)
-               print*,'aerosol=',aerosol(ig,:,:)
                print*,'temp=   ',pt(ig,:)
                print*,'pplay=  ',pplay(ig,:)
@@ -1268,19 +1216,4 @@
          call abort
       endif
-
-
-      ! Compute column amounts (kg m-2) if tracers are enabled.
-      if(tracer)then
-         qcol(1:ngrid,1:nq)=0.0
-         do iq=1,nq
-            do ig=1,ngrid
-               qcol(ig,iq) = SUM( zq(ig,1:nlayer,iq) * mass(ig,1:nlayer))
-            enddo
-         enddo
-
-         ! Generalised for arbitrary aerosols now. By LK
-         reffcol(1:ngrid,1:naerkind)=0.0
-
-      endif ! end of 'tracer'
 
 
@@ -1359,6 +1292,4 @@
                call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
                            'kg m^-2',2,qsurf(1,iq) )
-               call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
-                          'kg m^-2',2,qcol(1,iq) )
                           
 !              call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
@@ -1472,18 +1403,7 @@
             call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
                              'kg m^-2',2,qsurf_hist(1,iq) )
-            call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
-                            'kg m^-2',2,qcol(1,iq) )
                           
 !          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
 !                          'kg m^-2',2,qsurf(1,iq) )                          
-
-            do ig=1,ngrid
-               if(tau_col(ig).gt.1.e3)then
-                  print*,'WARNING: tau_col=',tau_col(ig)
-                  print*,'at ig=',ig,'in PHYSIQ'
-               endif         
-            end do
-
-            call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
 
          enddo ! end of 'nq' loop
Index: trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 1788)
@@ -1,5 +1,4 @@
 module radcommon_h
-      use radinc_h, only: L_NSPECTI, L_NSPECTV, L_NGAUSS, NTstar, NTstop, &
-                          naerkind, nsizemax
+      use radinc_h, only: L_NSPECTI, L_NSPECTV, L_NGAUSS, NTstar, NTstop
       implicit none
 
@@ -41,21 +40,4 @@
 !                  each temperature, pressure, and spectral interval
 !
-!     AEROSOL RADIATIVE OPTICAL CONSTANTS
-! 
-!   Shortwave
-!   ~~~~~~~~~
-! 
-! For the "naerkind" kind of aerosol radiative properties : 
-! QVISsQREF  :  Qext / Qext("longrefvis")
-! omegavis   :  single scattering albedo
-! gvis       :  assymetry factor
-! 
-!   Longwave
-!   ~~~~~~~~
-! 
-! For the "naerkind" kind of aerosol radiative properties : 
-! QIRsQREF :  Qext / Qext("longrefvis")
-! omegaIR  :  mean single scattering albedo
-! gIR      :  mean assymetry factor
 
       REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
@@ -65,8 +47,4 @@
 	!$OMP WNOV,DWNV,WAVEV,&
 	!$OMP STELLARF,TAURAY)
-
-      REAL*8 blami(L_NSPECTI+1)
-      REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90
-!$OMP THREADPRIVATE(blami,blamv)
 
       !! AS: introduced to avoid doing same computations again for continuum
@@ -85,37 +63,4 @@
 	!$OMP FZEROI,FZEROV)     !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
 
-      real QVISsQREF(L_NSPECTV,naerkind,nsizemax)
-      real omegavis(L_NSPECTV,naerkind,nsizemax)
-      real gvis(L_NSPECTV,naerkind,nsizemax)
-      real QIRsQREF(L_NSPECTI,naerkind,nsizemax)
-      real omegair(L_NSPECTI,naerkind,nsizemax)
-      real gir(L_NSPECTI,naerkind,nsizemax)
-!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir)
-
-
-! Reference wavelengths used to compute reference optical depth (m)
-! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      REAL lamrefir(naerkind),lamrefvis(naerkind)
-
-! Actual number of grain size classes in each domain for a
-!   given aerosol:
-
-      INTEGER          :: nsize(naerkind,2)
-
-! Particle size axis (depend on the kind of aerosol and the
-!   radiation domain)
-
-      DOUBLE PRECISION :: radiustab(naerkind,2,nsizemax)
-!$OMP THREADPRIVATE(lamrefir,lamrefvis,radiustab) !nsize read by suaer_corrk
-
-! Extinction coefficient at reference wavelengths;
-!   These wavelengths are defined in aeroptproperties, and called
-!   longrefvis and longrefir.
-
-      REAL :: QREFvis(naerkind,nsizemax)
-      REAL :: QREFir(naerkind,nsizemax)
-      REAL :: omegaREFvis(naerkind,nsizemax)
-      REAL :: omegaREFir(naerkind,nsizemax)
 
       REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
Index: trunk/LMDZ.TITAN/libf/phytitan/radii_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radii_mod.F90	(revision 1787)
+++ 	(revision )
@@ -1,121 +1,0 @@
-!==================================================================
-module radii_mod
-!==================================================================
-!  module to centralize the radii calculations for aerosols
-!==================================================================
-      
-      use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo, &
-                size_tropo,pres_bottom_strato,size_strato
-
-contains
-
-
-!==================================================================
-   subroutine su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of liquid and icy water particles
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-      use ioipsl_getin_p_mod, only: getin_p
-      use radinc_h, only: naerkind
-      use aerosol_mod, only: iaero_back2lay
-      Implicit none
-
-      integer,intent(in) :: ngrid
-      integer,intent(in) :: nlayer
-
-      real, intent(out) :: reffrad(ngrid,nlayer,naerkind)      !aerosols radii (K)
-      real, intent(out) :: nueffrad(ngrid,nlayer,naerkind)     !variance     
-
-      logical, save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-      integer :: iaer   
-      
-      print*,'enter su_aer_radii'
-          do iaer=1,naerkind
-!     these values will change once the microphysics gets to work
-!     UNLESS tracer=.false., in which case we should be working with
-!     a fixed aerosol layer, and be able to define reffrad in a 
-!     .def file. To be improved!
-
-
-!     WARNING : Titan adapt. (J. Vatant d'Ollone, 2017)
-!            - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY !
-!            - This routine is just here to keep the code running without unplugging all (yet)
-!            - There's only the dummy aerosol case on iaer = 1      
-            if(iaer.eq.1)then
-               reffrad(1:ngrid,1:nlayer,iaer) = 1.e-4
-               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-            endif
-            
-            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
-               reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
-               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
-            endif
-
-            if(iaer.gt.5)then
-               print*,'Error in callcorrk, naerkind is too high (>5).'
-               print*,'The code still needs generalisation to arbitrary'
-               print*,'aerosol kinds and number.'
-               call abort
-            endif
-
-         enddo
-
-      print*,'exit su_aer_radii'
-
-   end subroutine su_aer_radii
-!==================================================================
-
-
-!==================================================================
-   subroutine back2lay_reffrad(ngrid,reffrad,nlayer,pplev)
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the effective radii of particles in a 2-layer model
-!
-!     Authors
-!     -------
-!     Sandrine Guerlet (2013)
-!
-!==================================================================
- 
-      use aerosol_mod   !! Particle sizes and boundaries of aerosol layers defined there
-     Implicit none
-
-      integer,intent(in) :: ngrid
-
-      real, intent(out) :: reffrad(ngrid,nlayer)      ! particle radii (m)
-      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
-      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
-      REAL :: expfactor
-      INTEGER l,ig
-            
-      reffrad(:,:)=1e-6  !!initialization, not important
-          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
-                reffrad(ig,l) = size_tropo
-              ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
-                expfactor=log(size_strato/size_tropo) / log(pres_bottom_strato/pres_top_tropo)
-                reffrad(ig,l)= size_tropo*((pplev(ig,l)/pres_top_tropo)**expfactor)
-              ELSEIF (pplev(ig,l) .le. pres_bottom_strato) then
-                reffrad(ig,l) = size_strato
-              ENDIF
-            ENDDO
-          ENDDO
-
-   end subroutine back2lay_reffrad
-!==================================================================
-
-
-
-end module radii_mod
-!==================================================================
Index: trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/radinc_h.F90	(revision 1788)
@@ -49,9 +49,4 @@
 !               the k-coefficients. Variable component of the mixture
 !		can in princple be anything: currently it's H2O.
-!
-!     NAERKIND  The number of radiatively active aerosol types
-!
-!     NSIZEMAX  The maximum number of aerosol particle sizes
-!
 !----------------------------------------------------------------------
 
@@ -70,5 +65,4 @@
       integer, parameter :: L_NSPECTV = NBvisible
 
-!      integer, parameter :: NAERKIND  = 2 ! set in scatterers.h
       real,    parameter :: L_TAUMAX  = 35
 
@@ -84,9 +78,4 @@
       !integer, parameter :: NTstop = 50000
       !real*8,parameter :: NTfac = 1.0D+2    
-
-      ! Maximum number of grain size classes for aerosol convolution:
-      ! This must correspond to size of largest dataset used for aerosol 
-      ! optical properties in datagcm folder.
-      integer, parameter :: nsizemax = 60
 
       character(len=100),save :: corrkdir
Index: trunk/LMDZ.TITAN/libf/phytitan/setspi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/setspi.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/setspi.F90	(revision 1788)
@@ -23,5 +23,5 @@
 
       use radinc_h,    only: L_NSPECTI,corrkdir,banddir,NTstar,NTstop,NTfac
-      use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma
+      use radcommon_h, only: BWNI,WNOI,DWNI,WAVEI,planckir,sigma
       use datafile_mod, only: datadir
       use comcstfi_mod, only: pi
@@ -136,7 +136,5 @@
          DWNI(M)  = BWNI(M+1)-BWNI(M)
          WAVEI(M) = 1.0D+4/WNOI(M)
-         BLAMI(M) = 0.01D0/BWNI(M)         
-      end do
-      BLAMI(M) = 0.01D0/BWNI(M)
+      end do
 !     note M=L_NSPECTI+1 after loop due to Fortran bizarreness
 
Index: trunk/LMDZ.TITAN/libf/phytitan/setspv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/setspv.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/setspv.F90	(revision 1788)
@@ -24,5 +24,5 @@
 
       use radinc_h,    only: L_NSPECTV, corrkdir, banddir
-      use radcommon_h, only: BWNV,BLAMV,WNOV,DWNV,WAVEV, &
+      use radcommon_h, only: BWNV,WNOV,DWNV,WAVEV, &
                              STELLARF,TAURAY
       use datafile_mod, only: datadir
@@ -113,7 +113,5 @@
          DWNV(M)  = BWNV(M+1)-BWNV(M)
          WAVEV(M) = 1.0E+4/WNOV(M)
-         BLAMV(M) = 0.01/BWNV(M)
       end do
-      BLAMV(M) = 0.01/BWNV(M) ! wavelength in METERS for aerosol stuff
 !     note M=L_NSPECTV+1 after loop due to Fortran bizarreness
 
Index: trunk/LMDZ.TITAN/libf/phytitan/stokes.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/stokes.F90	(revision 1787)
+++ 	(revision )
@@ -1,68 +1,0 @@
-      subroutine stokes(p,t,rd,w,rho_aer)
-
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the sedimentation velocity in a low pressure 
-!     atmosphere.
-!
-!     Authors
-!     ------- 
-!     Francois Forget (1997)
-!
-!==================================================================
-
-!      use radcommon_h, only : Rgas
-      use comcstfi_mod, only: g, pi, avocado
-
-      implicit none
-
-!     input
-!     -----
-!     pressure (Pa), Temperature (K), particle radius (m), density
-      real p, t, rd, rho_aer 
-
-!     output
-!     ------
-!     sedimentation velocity (m/s, >0)
-      real w
-
-!     locally saved variables
-!     ---------------------
-      real a,b,molrad,visc
-      save a,b
-!$OMP THREADPRIVATE(a,b)
-  
-      LOGICAL firstcall
-      SAVE firstcall
-      DATA firstcall/.true./
-!$OMP THREADPRIVATE(firstcall)
-
-      if (firstcall) then
-
-         !print*,'Routine not working: replace Rgas with r'
-	 !stop
-         !a = 0.707*Rgas/(4*pi*molrad**2 * avocado)
-
-         molrad=2.2e-10   ! CO2 (only used in condense_co2cloud at the moment)
-         visc=1.e-5       ! CO2
-
-
-         a = 0.707*8.31/(4*pi*molrad**2 * avocado)
-         b = (2./9.) * rho_aer * g / visc 
-	!print*,'molrad=',molrad
-	!print*,'visc=',visc
-	!print*,'a=',a
-	!print*,'b=',b
-	!print*,'rho_aer=',rho_aer
-	!stop
- 
-         firstcall=.false.
-      end if
-
-!     Sedimentation velocity =
-!     Stokes' Law corrected for low pressures by the Cunningham
-!     slip-flow correction according to Rossow (Icarus 36, 1-50, 1978)
-      w = b * rd*rd * (1 + 1.333* (a*T/P)/rd )  
-      return
-    end subroutine stokes
Index: trunk/LMDZ.TITAN/libf/phytitan/suaer_corrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/suaer_corrk.F90	(revision 1787)
+++ 	(revision )
@@ -1,433 +1,0 @@
-      subroutine suaer_corrk
-
-      ! inputs
-      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
-      use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
-      use datafile_mod, only: datadir, aerdir
-
-      ! outputs
-      use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
-      use radcommon_h, only: radiustab,nsize,tstellar
-      use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
-      use aerosol_mod
-      use callkeys_mod, only: tplanet
-
-      implicit none
-
-!==================================================================
-!     Purpose
-!     -------
-!     Initialize all radiative aerosol properties
-!
-!     Notes
-!     -----
-!     Reads the optical properties -> Mean  -> Variable assignment
-!     (ASCII files)                  (see radcommon_h.F90)
-!     wvl(nwvl)                      longsun
-!     ep(nwvl)                       epav     QVISsQREF(L_NSPECTV)
-!     omeg(nwvl)                     omegav   omegavis(L_NSPECTV)
-!     gfactor(nwvl)                  gav      gvis(L_NSPECTV)
-!     
-!     Authors
-!     ------- 
-!     Richard Fournier (1996) Francois Forget (1996)
-!     Frederic Hourdin
-!     Jean-jacques morcrette *ECMWF*
-!     MODIF Francois Forget (2000)
-!     MODIF Franck Montmessin (add water ice)
-!     MODIF J.-B. Madeleine 2008W27
-!     - Optical properties read in ASCII files
-!     - Add varying radius of the particules
-!     - Add water ice clouds
-!     MODIF R. Wordsworth (2009)
-!     - generalisation to arbitrary spectral bands 
-!     
-!     WARNING : Titan adapt. (J. Vatant d'Ollone, 2017)
-!            - ONLY THE NO AEROSOL CASE FOR NOW SINCE WE COMPUTE THEM ANOTHER WAY !
-!            - This routine is just here to keep the code running without unplugging all (yet)                           
-!================================================================================================
-
-!     Optical properties (read in external ASCII files)
-      INTEGER,SAVE      :: nwvl  ! Number of wavelengths in
-                                ! the domain (VIS or IR), read by master
-
-!      REAL             :: solsir ! visible to infrared ratio
-!                                ! (tau_0.67um/tau_9um)
-
-      REAL, DIMENSION(:),&
-      ALLOCATABLE, SAVE :: wvl  ! Wavelength axis, read by master
-      REAL, DIMENSION(:),&
-      ALLOCATABLE, SAVE :: radiusdyn ! Particle size axis, read by master
-
-      REAL, DIMENSION(:,:),&
-      ALLOCATABLE, SAVE :: ep,& ! Extinction coefficient Qext, read by master
-      omeg,&                    ! Single Scattering Albedo, read by master
-      gfactor                   ! Assymetry Factor, read by master
-
-!     Local variables:
-
-      INTEGER :: iaer           ! Aerosol index
-      INTEGER :: idomain        ! Domain index (1=VIS,2=IR)
-      INTEGER :: ilw            ! longwave index
-      INTEGER :: isw           ! shortwave index
-      INTEGER :: isize          ! Particle size index
-      INTEGER :: jfile          ! ASCII file scan index
-      INTEGER :: file_unit = 60
-      LOGICAL :: file_ok, endwhile
-      CHARACTER(LEN=132) :: scanline ! ASCII file scanning line
-
-!     I/O  of "aerave" (subroutine that spectrally averages
-!     the single scattering parameters)
-
-      REAL lamref                      ! reference wavelengths
-      REAL epref                       ! reference extinction ep
-
-!      REAL epav(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
-!      REAL omegav(L_NSPECTI)          ! Average single scattering albedo
-!      REAL gav(L_NSPECTI)             ! Average assymetry parameter
-
-      REAL epavVI(L_NSPECTV)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
-      REAL omegavVI(L_NSPECTV)          ! Average single scattering albedo
-      REAL gavVI(L_NSPECTV)             ! Average assymetry parameter
-
-      REAL epavIR(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
-      REAL omegavIR(L_NSPECTI)          ! Average single scattering albedo
-      REAL gavIR(L_NSPECTI)             ! Average assymetry parameter
-      
-      logical forgetit                  ! use francois' data?
-      integer iwvl
-
-!     Local saved variables:
-
-      CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id
-!$OMP THREADPRIVATE(file_id)
-!---- Please indicate the names of the optical property files below
-!     Please also choose the reference wavelengths of each aerosol
-
-      if (noaero) then
-        print*, 'naerkind= 0'
-      endif
-      do iaer=1,naerkind
-       if (iaer.eq.1) then ! JVO 2017 : Now iaer = 1 is always dummy co2 for noaero case, since we don't use aerosols
-        forgetit=.true.
-!     visible
-        if(forgetit)then
-           file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
-        else
-           file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
-        endif
-        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
-
-!     infrared
-        if(forgetit)then
-           file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
-        else
-           file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
-        endif
-        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
-       endif ! CO2 aerosols  
-! NOTE: these lamref's are currently for dust, not CO2 ice.
-! JB thinks this shouldn't matter too much, but it needs to be 
-! fixed / decided for the final version.
-
-       if (iaer.eq.iaero_back2lay) then
-         print*, 'naerkind= back2lay', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_saturn_vis_n20.dat'
-         lamrefvis(iaer)=0.8E-6  ! 
-!     infrared
-         file_id(iaer,2) = 'optprop_saturn_ir_n20.dat'
-         lamrefir(iaer)=6.E-6  ! 
-! added by SG
-       endif
-       
-       
-      enddo
-
-!------------------------------------------------------------------
-
-!     Initializations:
-
-      radiustab(:,:,:) = 0.0
-      gvis(:,:,:)      = 0.0
-      omegavis(:,:,:)  = 0.0
-      QVISsQREF(:,:,:) = 0.0
-      gir(:,:,:)       = 0.0
-      omegair(:,:,:)   = 0.0
-      QIRsQREF(:,:,:)  = 0.0
-
-      DO iaer = 1, naerkind     ! Loop on aerosol kind
-         DO idomain = 1, 2      ! Loop on radiation domain (VIS or IR)
-!==================================================================
-!     1. READ OPTICAL PROPERTIES
-!==================================================================
-
-!     1.1 Open the ASCII file
-
-!$OMP MASTER
-            INQUIRE(FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
-                    '/'//TRIM(file_id(iaer,idomain)),&
-                    EXIST=file_ok)
-            IF (file_ok) THEN
-              OPEN(UNIT=file_unit,&
-                   FILE=TRIM(datadir)//'/'//TRIM(aerdir)//&
-                        '/'//TRIM(file_id(iaer,idomain)),&
-                   FORM='formatted')
-            ELSE
-             ! In ye old days these files were stored in datadir;
-             ! let's be retro-compatible
-              INQUIRE(FILE=TRIM(datadir)//&
-                      '/'//TRIM(file_id(iaer,idomain)),&
-                      EXIST=file_ok)
-              IF (file_ok) THEN
-                OPEN(UNIT=file_unit,&
-                     FILE=TRIM(datadir)//&
-                          '/'//TRIM(file_id(iaer,idomain)),&
-                     FORM='formatted')
-              ENDIF              
-            ENDIF
-            IF(.NOT.file_ok) THEN
-               write(*,*)'suaer_corrk: Problem opening ',&
-               TRIM(file_id(iaer,idomain))
-               write(*,*)'It should be in: ',TRIM(datadir)//'/'//TRIM(aerdir)
-               write(*,*)'1) You can set this directory address ',&
-               'in callphys.def with:'
-               write(*,*)' datadir = /absolute/path/to/datagcm'
-               write(*,*)'2) If ',&
-              TRIM(file_id(iaer,idomain)),&
-               ' is a LMD reference datafile, it'
-               write(*,*)' can be obtained online at:'
-               write(*,*)' http://www.lmd.jussieu.fr/',&
-               '~lmdz/planets/LMDZ.GENERIC/datagcm/'
-               CALL ABORT 
-            ENDIF
-
-!     1.2 Allocate the optical property table
-
-            jfile = 1
-            endwhile = .false.
-            DO WHILE (.NOT.endwhile)
-               READ(file_unit,*) scanline
-               IF ((scanline(1:1) .ne. '#').and.&
-               (scanline(1:1) .ne. ' ')) THEN
-               BACKSPACE(file_unit)
-               reading1_seq: SELECT CASE (jfile) ! ====================
-            CASE(1) reading1_seq ! nwvl ----------------------------
-               read(file_unit,*) nwvl
-               jfile = jfile+1
-            CASE(2) reading1_seq ! nsize ---------------------------
-               read(file_unit,*) nsize(iaer,idomain)
-               endwhile = .true.
-            CASE DEFAULT reading1_seq ! ---------------------------- 
-               WRITE(*,*) 'readoptprop: ',& 
-               'Error while loading optical properties.' 
-               CALL ABORT 
-            END SELECT reading1_seq ! ==============================
-         ENDIF
-      ENDDO
-
-      ALLOCATE(wvl(nwvl))       ! wvl
-      ALLOCATE(radiusdyn(nsize(iaer,idomain))) ! radiusdyn
-      ALLOCATE(ep(nwvl,nsize(iaer,idomain))) ! ep
-      ALLOCATE(omeg(nwvl,nsize(iaer,idomain))) ! omeg 
-      ALLOCATE(gfactor(nwvl,nsize(iaer,idomain))) ! g
-
-
-!     1.3 Read the data
-
-      jfile = 1
-      endwhile = .false.
-      DO WHILE (.NOT.endwhile)
-         READ(file_unit,*) scanline
-         IF ((scanline(1:1) .ne. '#').and.&
-         (scanline(1:1) .ne. ' ')) THEN
-         BACKSPACE(file_unit)
-         reading2_seq: SELECT CASE (jfile) ! ====================
-      CASE(1) reading2_seq      ! wvl -----------------------------
-         read(file_unit,*) wvl
-         jfile = jfile+1
-      CASE(2) reading2_seq      ! radiusdyn -----------------------
-         read(file_unit,*) radiusdyn
-         jfile = jfile+1
-      CASE(3) reading2_seq      ! ep ------------------------------
-         isize = 1
-         DO WHILE (isize .le. nsize(iaer,idomain))
-            READ(file_unit,*) scanline
-            IF ((scanline(1:1) .ne. '#').and.&
-            (scanline(1:1) .ne. ' ')) THEN
-            BACKSPACE(file_unit)
-            read(file_unit,*) ep(:,isize)
-            isize = isize + 1
-         ENDIF
-      ENDDO
-
-      jfile = jfile+1
-      CASE(4) reading2_seq      ! omeg ----------------------------
-         isize = 1
-         DO WHILE (isize .le. nsize(iaer,idomain))
-            READ(file_unit,*) scanline
-            IF ((scanline(1:1) .ne. '#').and.&
-            (scanline(1:1) .ne. ' ')) THEN
-            BACKSPACE(file_unit)
-            read(file_unit,*) omeg(:,isize)
-            isize = isize + 1
-         ENDIF
-      ENDDO
-
-      jfile = jfile+1
-      CASE(5) reading2_seq      ! gfactor -------------------------
-         isize = 1
-         DO WHILE (isize .le. nsize(iaer,idomain))
-            READ(file_unit,*) scanline
-            IF ((scanline(1:1) .ne. '#').and.&
-            (scanline(1:1) .ne. ' ')) THEN
-            BACKSPACE(file_unit)
-            read(file_unit,*) gfactor(:,isize)
-            isize = isize + 1
-         ENDIF
-      ENDDO
-
-      jfile = jfile+1
-      IF ((idomain.NE.1).OR.(iaer.NE.1)) THEN
-         endwhile = .true.
-      ENDIF
-      CASE(6) reading2_seq
-         endwhile = .true.
-      CASE DEFAULT reading2_seq ! ----------------------------
-         WRITE(*,*) 'readoptprop: ',&
-         'Error while loading optical properties.'
-         CALL ABORT
-      END SELECT reading2_seq   ! ==============================
-      ENDIF
-      ENDDO
-
-!     1.4 Close the file
-
-      CLOSE(file_unit)
-
-!     1.5 If Francois' data, convert wvl to metres
-       if(iaer.eq.1)then
-         if(forgetit)then
-         !   print*,'please sort out forgetit for naerkind>1'
-            do iwvl=1,nwvl
-               wvl(iwvl)=wvl(iwvl)*1.e-6
-            enddo
-         endif
-      endif
-
-!$OMP END MASTER
-!$OMP BARRIER
-
-
-
-
-
-!==================================================================
-!     2. AVERAGED PROPERTIES AND VARIABLE ASSIGNMENTS
-!==================================================================
-      domain: SELECT CASE (idomain)
-!==================================================================
-      CASE(1) domain            !       VISIBLE DOMAIN (idomain=1)
-!==================================================================
-
-         lamref=lamrefvis(iaer)
-         epref=1.E+0
-
-         DO isize=1,nsize(iaer,idomain)
-
-!     Save the particle sizes
-            radiustab(iaer,idomain,isize)=radiusdyn(isize)
-
-!     Averaged optical properties (GCM channels)
-
-!            CALL aerave ( nwvl,&
-!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
-!            lamref,epref,tstellar,&
-!            L_NSPECTV,blamv,epavVI,&
-!            omegavVI,gavVI,QREFvis(iaer,isize))!,omegaREFir(iaer,isize))
-            CALL aerave_new ( nwvl,&
-            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
-            lamref,epref,tstellar,&
-            L_NSPECTV,blamv,epavVI,&
-            omegavVI,gavVI,QREFvis(iaer,isize),omegaREFir(iaer,isize))
-
-!     Variable assignments (declared in radcommon)
-            DO isw=1,L_NSPECTV
-               QVISsQREF(isw,iaer,isize)=epavVI(isw)
-               gvis(isw,iaer,isize)=gavVI(isw)
-               omegavis(isw,iaer,isize)=omegavVI(isw)
-            END DO
-
-         ENDDO
-
-!==================================================================
-      CASE(2) domain            !      INFRARED DOMAIN (idomain=2)
-!==================================================================
-
-
-         DO isize=1,nsize(iaer,idomain) ! ----------------------------------
-
-            lamref=lamrefir(iaer)
-            epref=1.E+0
-
-!     Save the particle sizes
-            radiustab(iaer,idomain,isize)=radiusdyn(isize)
-
-!     Averaged optical properties (GCM channels)
-
-!     epav is <QIR>/Qext(lamrefir) since epref=1
-!     Note: aerave also computes the extinction coefficient at
-!     the reference wavelength. This is called QREFvis or QREFir
-!     (not epref, which is a different parameter).
-!     Reference wavelengths SHOULD BE defined for each aerosol in
-!     radcommon_h.F90
-
-!            CALL aerave ( nwvl,&
-!            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
-!            lamref,epref,tplanet,&
-!            L_NSPECTI,blami,epavIR,&
-!            omegavIR,gavIR,QREFir(iaer,isize))!,omegaREFir(iaer,isize))
-            CALL aerave_new ( nwvl,&
-            wvl(:),ep(:,isize),omeg(:,isize),gfactor(:,isize),&
-            lamref,epref,tplanet,&
-            L_NSPECTI,blami,epavIR,&
-            omegavIR,gavIR,QREFir(iaer,isize),omegaREFir(iaer,isize))
-
-
-!     Variable assignments (declared in radcommon)
-            DO ilw=1,L_NSPECTI
-               QIRsQREF(ilw,iaer,isize)=epavIR(ilw)
-               gir(ilw,iaer,isize)=gavIR(ilw)
-               omegair(ilw,iaer,isize)=omegavIR(ilw)
-            END DO
-
-
-      ENDDO                     ! isize (particle size) -------------------------------------
-
-      END SELECT domain
-
-!========================================================================
-!     3. Deallocate temporary variables that were read in the ASCII files
-!========================================================================
-
-!$OMP BARRIER
-!$OMP MASTER
-      IF (ALLOCATED(wvl)) DEALLOCATE(wvl)                 ! wvl
-      IF (ALLOCATED(radiusdyn)) DEALLOCATE(radiusdyn)     ! radiusdyn
-      IF (ALLOCATED(ep)) DEALLOCATE(ep)                   ! ep 
-      IF (ALLOCATED(omeg)) DEALLOCATE(omeg)          	  ! omeg 
-      IF (ALLOCATED(gfactor)) DEALLOCATE(gfactor)         ! g
-!$OMP END MASTER
-!$OMP BARRIER
-
-      END DO                    ! Loop on iaer
-      END DO                    ! Loop on idomain
-!========================================================================
-
-      RETURN
-
-
-
-    END subroutine suaer_corrk
-      
Index: trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1787)
+++ trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1788)
@@ -10,12 +10,7 @@
        character*20, DIMENSION(:), ALLOCATABLE :: noms   ! name of the tracer
        real, DIMENSION(:), ALLOCATABLE :: mmol     ! mole mass of tracer (g/mol-1) 
-       real, DIMENSION(:), ALLOCATABLE :: radius   ! dust and ice particle radius (m)
        real, DIMENSION(:), ALLOCATABLE :: rho_q    ! tracer densities (kg.m-3)
-       real, DIMENSION(:), ALLOCATABLE :: qext     ! Single Scat. Extinction coeff at 0.67 um
-       real, DIMENSION(:), ALLOCATABLE :: alpha_lift  ! saltation vertical flux/horiz flux ratio (m-1)
-       real, DIMENSION(:), ALLOCATABLE :: alpha_devil ! lifting coeeficient by dust devil
-       real, DIMENSION(:), ALLOCATABLE :: qextrhor ! Intermediate for computing opt. depth from q
 
-!$OMP THREADPRIVATE(noms,mmol,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor)
+!$OMP THREADPRIVATE(noms,mmol,rho_q)
 
 ! tracer indexes: these are initialized in initracer and should be 0 if the
Index: trunk/LMDZ.TITAN/libf/phytitan/vlz_fi.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/vlz_fi.F	(revision 1787)
+++ 	(revision )
@@ -1,190 +1,0 @@
-      SUBROUTINE vlz_fi(ngrid,nlayer,q,pente_max,masse,w,wq)
-c
-c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
-c
-c    ********************************************************************
-c     Shema  d'advection " pseudo amont " dans la verticale
-c    pour appel dans la physique (sedimentation)
-c    ********************************************************************
-c    q rapport de melange (kg/kg)...
-c    masse : masse de la couche Dp/g
-c    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
-c    pente_max = 2 conseillee
-c
-c
-c   --------------------------------------------------------------------
-      IMPLICIT NONE
-c
-c
-c   Arguments:
-c   ----------
-      integer,intent(in) :: ngrid, nlayer
-      real,intent(in) :: masse(ngrid,nlayer),pente_max
-      REAL,INTENT(INOUT) :: q(ngrid,nlayer)
-      REAL,INTENT(INOUT) :: w(ngrid,nlayer)
-      REAL,INTENT(OUT) :: wq(ngrid,nlayer+1)
-c
-c      Local 
-c   ---------
-c
-      INTEGER i,ij,l,j,ii
-c
-
-      real dzq(ngrid,nlayer),dzqw(ngrid,nlayer),adzqw(ngrid,nlayer)
-      real dzqmax
-      real newmasse
-      real sigw, Mtot, MQtot
-      integer m
-
-      REAL      SSUM,CVMGP,CVMGT
-      integer ismax,ismin
-
-
-c    On oriente tout dans le sens de la pression c'est a dire dans le
-c    sens de W
-
-      do l=2,nlayer
-         do ij=1,ngrid
-            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
-            adzqw(ij,l)=abs(dzqw(ij,l))
-         enddo
-      enddo
-
-      do l=2,nlayer-1
-         do ij=1,ngrid
-#ifdef CRAY
-            dzq(ij,l)=0.5*
-     ,      cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
-#else
-            if(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) then
-                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
-            else
-                dzq(ij,l)=0.
-            endif
-#endif
-            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
-            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
-         enddo
-      enddo
-
-      do ij=1,ngrid
-         dzq(ij,1)=0.
-         dzq(ij,nlayer)=0.
-      enddo
-c ---------------------------------------------------------------
-c   .... calcul des termes d'advection verticale  .......
-c ---------------------------------------------------------------
-
-c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
-c
-c      No flux at the model top:
-       do ij=1,ngrid
-          wq(ij,nlayer+1)=0.
-       enddo
-
-c      1) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)     
-c      ===============================
-
-       do l = 1,nlayer          ! loop different than when w<0
-        do ij=1,ngrid
-
-         if(w(ij,l).gt.0.)then
-
-c         Regular scheme (transfered mass < 1 layer)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          if(w(ij,l).le.masse(ij,l))then
-            sigw=w(ij,l)/masse(ij,l)
-            wq(ij,l)=w(ij,l)*(q(ij,l)+0.5*(1.-sigw)*dzq(ij,l))
-            
-
-c         Extended scheme (transfered mass > 1 layer)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          else 
-            m=l
-            Mtot = masse(ij,m)
-            MQtot = masse(ij,m)*q(ij,m)
-            if(m.ge.nlayer)goto 88
-            do while(w(ij,l).gt.(Mtot+masse(ij,m+1)))
-                m=m+1
-                Mtot = Mtot + masse(ij,m)
-                MQtot = MQtot + masse(ij,m)*q(ij,m)
-                if(m.ge.nlayer)goto 88
-            end do
- 88         continue
-            if (m.lt.nlayer) then
-                sigw=(w(ij,l)-Mtot)/masse(ij,m+1)
-                wq(ij,l)=(MQtot + (w(ij,l)-Mtot)*
-     &          (q(ij,m+1)+0.5*(1.-sigw)*dzq(ij,m+1)) )
-            else
-                w(ij,l) = Mtot
-                wq(ij,l) = Mqtot 
-            end if
-          end if
-         end if
-        enddo
-       enddo
-
-c      2) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)     
-c      ===============================
-       goto 99 ! SKIPPING THIS PART FOR SEDIMENTATION 
-
-c      Surface flux up:
-       do ij=1,ngrid
-         if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid
-       end do
-
-       do l = 1,nlayer-1  ! loop different than when w>0
-        do ij=1,ngrid
-         if(w(ij,l+1).le.0)then
-
-c         Regular scheme (transfered mass < 1 layer)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          if(-w(ij,l+1).le.masse(ij,l))then
-            sigw=w(ij,l+1)/masse(ij,l)
-            wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
-c         Extended scheme (transfered mass > 1 layer)
-c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          else 
-             m = l-1
-             Mtot = masse(ij,m+1)
-             MQtot = masse(ij,m+1)*q(ij,m+1)
-             if (m.le.0)goto 77
-             do while(-w(ij,l+1).gt.(Mtot+masse(ij,m)))
-                m=m-1
-                Mtot = Mtot + masse(ij,m+1)
-                MQtot = MQtot + masse(ij,m+1)*q(ij,m+1)
-                if (m.le.0)goto 77
-             end do
- 77          continue
-
-             if (m.gt.0) then
-                sigw=(w(ij,l+1)+Mtot)/masse(ij,m)
-                wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*
-     &          (q(ij,m)-0.5*(1.+sigw)*dzq(ij,m))  )
-             else
-c               wq(ij,l+1)= (MQtot + (-w(ij,l+1)-Mtot)*qm(ij,1))
-                write(*,*) 'a rather weird situation in vlz_fi !'
-                stop
-             end if
-          endif
-         endif
-        enddo
-       enddo
- 99    continue
-
-      do l=1,nlayer
-         do ij=1,ngrid
-
-cccccccc lines below not used for sedimentation (No real flux)
-ccccc       newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l) 
-ccccc       q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
-ccccc&         /newmasse
-ccccc       masse(ij,l)=newmasse
-
-            q(ij,l)=q(ij,l) +  (wq(ij,l+1)-wq(ij,l))/masse(ij,l)
-
-         enddo
-      enddo
-
-
-      end
