Index: trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90	(revision 3585)
@@ -5,13 +5,14 @@
 contains
 
-      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
-         aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col)
+      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq, &
+         dtau_aer,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col)
 
        use radinc_h, only : L_TAUMAX,naerkind
        use aerosol_mod, only: iaero_haze, i_haze
-       USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol
+       USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol,micro_indx
        use comcstfi_mod, only: g, pi, mugaz, avocado
        use geometry_mod, only: latitude
-       use callkeys_mod, only: kastprof
+       use callkeys_mod, only: kastprof, callmufi
+       use mp2m_diagnostics
        implicit none
 
@@ -28,6 +29,4 @@
 !     update J.-B. Madeleine (2008)
 !     dust removal, simplification by Robin Wordsworth (2009)
-!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
-!     Radiative Generic Condensable Species - Lucas Teinturier (2022)
 !
 !     Input
@@ -44,5 +43,5 @@
 !     Output
 !     ------
-!     aerosol            Aerosol optical depth in layer l, grid point ig
+!     dtau_aer           Aerosol optical depth in layer l, grid point ig
 !     tau_col            Total column optical depth at grid point ig
 !
@@ -54,9 +53,10 @@
       REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: zzlev(ngrid,nlayer)   ! Altitude at the layer boundaries.
       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
       REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
-      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
-      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance
+      REAL,INTENT(OUT) :: dtau_aer(ngrid,nlayer,naerkind) ! aerosol optical depth
+      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind)   ! aerosol effective radius
+      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind)  ! aerosol effective variance
       REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
       REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
@@ -69,4 +69,9 @@
       real dp_tropo(ngrid)
       real dp_layer(ngrid)
+
+      ! Microphysical tracers
+      real sig
+      real m0as(ngrid,nlayer)
+      real m0af(ngrid,nlayer)
 
       INTEGER l,ig,iq,iaer,ia
@@ -110,28 +115,47 @@
 
       if (iaero_haze.ne.0) then
-        iaer=iaero_haze
-!       1. Initialization
-         aerosol(1:ngrid,1:nlayer,iaer)=0.0
-!       2. Opacity calculation
-         DO ig=1, ngrid
-               DO l=1,nlayer-1 ! to stop the rad tran bug
-                  ! if fractal, radius doit etre equivalent sphere radius
-                  aerosol0 =                         &
-                       (  0.75 * QREFvis3d(ig,l,iaer) /        &
-                       ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
-                       ( pq(ig,l,i_haze) + 1.E-10 ) *         &
-                       ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                  aerosol0           = max(aerosol0,1.e-10)
-                  aerosol0           = min(aerosol0,L_TAUMAX)
-                  aerosol(ig,l,iaer) = aerosol0
+         if (callmufi) then
+            ! Convert intensive microphysical tracers to extensive [m-2]
+            m0as(:,:) = pq(:,:,micro_indx(1)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
+            m0af(:,:) = pq(:,:,micro_indx(3)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g
+            
+            ! Spherical aerosols
+            sig = 0.2
+            dtau_aer(:,:,1) = m0as(:,:) * QREFvis3d(:,:,1) * pi * mp2m_rc_sph(:,:)**2 * exp(2*sig**2)
+            ! Fractal aerosols
+            sig = 0.35
+            dtau_aer(:,:,2) = m0af(:,:) * QREFvis3d(:,:,2) * pi * mp2m_rc_fra(:,:)**2 * exp(2*sig**2)
+            
+            ! write(*,*) 'dtau_as  :', MINVAL(dtau_aer(:,:,1)), '-', MAXVAL(dtau_aer(:,:,1))
+            ! write(*,*) 'dtau_af  :', MINVAL(dtau_aer(:,:,2)), '-', MAXVAL(dtau_aer(:,:,2))
+         
+         else
+            do iaer = 1, naerkind
+               ! 1. Initialization
+               dtau_aer(1:ngrid,1:nlayer,iaer)=0.0
+               ! 2. Opacity calculation
+               DO ig = 1, ngrid
+                  DO l = 1, nlayer-1 ! to stop the rad tran bug
+                     ! If fractal, radius doit etre equivalent sphere radius
+                     ! Eq. 2.37 - Madeleine's PhD (2011).
+                     aerosol0 =                         &
+                        (  0.75 * QREFvis3d(ig,l,iaer) /        &
+                        ( rho_q(i_haze) * reffrad(ig,l,iaer) )  ) *   &
+                        ( pq(ig,l,i_haze) + 1.E-10 ) *         &
+                        ( pplev(ig,l) - pplev(ig,l+1) ) / g
+                     aerosol0            = max(aerosol0,1.e-10)
+                     aerosol0            = min(aerosol0,L_TAUMAX)
+                     dtau_aer(ig,l,iaer) = aerosol0
+                  ENDDO
                ENDDO
-         ENDDO
-         !QREF est le meme dans toute la colonne par def si size uniforme
-         !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)
-         !print*, 'TB17: rho_q=',rho_q(i_haze)
-         !print*, 'TB17: reffrad=',reffrad(1,:,1)
-         !print*, 'TB17: pq=',pq(1,:,i_haze)
-         !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
-   end if ! if haze aerosols
+               !QREF est le meme dans toute la colonne par def si size uniforme
+               !print*, 'TB17: QREFvis3d=',QREFvis3d(1,:,1)
+               !print*, 'TB17: rho_q=',rho_q(i_haze)
+               !print*, 'TB17: reffrad=',reffrad(1,:,1)
+               !print*, 'TB17: pq=',pq(1,:,i_haze)
+               !print*, 'TB17: deltap=',pplev(1,1) - pplev(1,nlayer)
+            enddo ! end iaer = 1, naerkind
+         endif ! end callmufi
+      endif ! if haze aerosols
 
 ! --------------------------------------------------------------------------
@@ -142,5 +166,5 @@
       do l=1,nlayer
          do ig=1,ngrid
-            tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
+            tau_col(ig) = tau_col(ig) + dtau_aer(ig,l,iaer)
          end do
       end do
@@ -150,6 +174,6 @@
       do l=1,nlayer
          do iaer = 1, naerkind
-            if(aerosol(ig,l,iaer).gt.1.e3)then
-               print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
+            if(dtau_aer(ig,l,iaer).gt.1.e3)then
+               print*,'WARNING: dtau_aer=',dtau_aer(ig,l,iaer)
                print*,'at ig=',ig,',  l=',l,', iaer=',iaer
                print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
@@ -164,5 +188,5 @@
          print*,'WARNING: tau_col=',tau_col(ig)
          print*,'at ig=',ig
-         print*,'aerosol=',aerosol(ig,:,:)
+         print*,'dtau_aer=',dtau_aer(ig,:,:)
          print*,'QREFvis3d=',QREFvis3d(ig,:,:)
          print*,'reffrad=',reffrad(ig,:,:)
Index: trunk/LMDZ.PLUTO/libf/phypluto/aeroptproperties.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/aeroptproperties.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/aeroptproperties.F90	(revision 3585)
@@ -279,9 +279,8 @@
               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)
+              !omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
+              !omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
             ENDDO
           ENDDO
-
 
           if (firstcall) then
@@ -297,9 +296,11 @@
       IF (firstcall) THEN
 
-!       1.1 Pi!
+        ! 1.1 Pi
+        !-------
         pi = 2. * asin(1.e0)
 
-!       1.2 Effective radius
-        refftab(1)    = refftabmin
+        ! 1.2 Effective radius
+        !---------------------
+        refftab(1)           = refftabmin
         refftab(refftabsize) = refftabmax
 
@@ -311,5 +312,6 @@
         enddo
 
-!       1.3 Effective variance
+        ! 1.3 Effective variance
+        !-----------------------
         if(nuefftabsize.eq.1)then ! addded by RDW
            print*,'Warning: no variance range in aeroptproperties'
@@ -344,9 +346,11 @@
       ENDIF ! of IF (firstcall)
 
-!       1.4 Radius middle point and range for Gauss integration
+        ! 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:
+        ! 1.5 Interpolating data at the Gauss quadrature points:
+        !-------------------------------------------------------
         DO gausind=1,ngau
           drad=radiusr*radgaus(gausind)
@@ -829,6 +833,4 @@
 !==================================================================
 
-
-
       ENDDO ! idomain
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90	(revision 3585)
@@ -7,5 +7,5 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
-          zzlay,tsurf,fract,dist_star,aerosol,muvar,           &
+          zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar,    &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
@@ -36,5 +36,6 @@
                               optichaze,haze_radproffix,&
                               methane,carbox,cooling,nlte,strobel,&
-                              ch4fix,vmrch4_proffix,vmrch4fix
+                              ch4fix,vmrch4_proffix,vmrch4fix,&
+                              callmufi
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -44,4 +45,5 @@
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
       use radcommon_h, only: wavev,wavei
+      use mp2m_diagnostics
       implicit none
 
@@ -81,4 +83,5 @@
       REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
       REAL,INTENT(IN) :: zzlay(ngrid,nlayer)       ! Mid-layer altitude
+      REAL,INTENT(IN) :: zzlev(ngrid,nlayer)       ! Altitude at the layer boundaries.
       REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
       REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
@@ -89,18 +92,18 @@
 
       ! OUTPUT
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
-      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.
-      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)             ! Incident LW flux to surf (W/m2).
-      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)             ! Incident SW flux to surf (W/m2)
-      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)          ! Absorbed SW flux by the surface (W/m2). By MT2015.
-      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)              ! Outgoing LW flux to space (W/m2).
-      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)              ! SW flux absorbed by the planet (W/m2).
-      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)              ! Incident top of atmosphere SW flux (W/m2).
-      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)        ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1).
-      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1).
-      REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV)        ! Surface SW radiation 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
+      REAL,INTENT(OUT) :: dtau_aer(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
+      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.
+      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)              ! Incident LW flux to surf (W/m2).
+      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)              ! Incident SW flux to surf (W/m2)
+      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)           ! Absorbed SW flux by the surface (W/m2). By MT2015.
+      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)               ! Outgoing LW flux to space (W/m2).
+      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)               ! SW flux absorbed by the planet (W/m2).
+      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)               ! Incident top of atmosphere SW flux (W/m2).
+      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)         ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1).
+      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)         ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1).
+      REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV)         ! Surface SW radiation 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
       REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
       REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
@@ -186,7 +189,6 @@
 !$OMP THREADPRIVATE(QREFvis3d,QREFir3d)
 
-
       ! Miscellaneous :
-      real*8  temp,temp1,temp2,pweight
+      real*8  temp,temp1,temp2,pweight,sig
       character(len=10) :: tmp1
       character(len=10) :: tmp2
@@ -331,7 +333,4 @@
          if (is_master) call system('rm -f surf_vals_long.out')
 
-         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-
-
 !--------------------------------------------------
 !             Set up correlated k
@@ -454,24 +453,47 @@
 !     Effective radius and variance of the aerosols
 !--------------------------------------------------
+      ! Radiative Hazes
       if (optichaze) then
-         do iaer=1,naerkind
-            if ((iaer.eq.iaero_haze)) then
-               call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer), &
-                            nueffrad(1,1,iaer))
+         if (callmufi) then
+            ! Spherical aerosols
+            sig = 0.2
+            where (mp2m_rc_sph(:,:) > 1e-10)
+               reffrad(:,:,1) = mp2m_rc_sph(:,:) * exp(5.*sig**2 / 2.)
+            elsewhere
+               reffrad(:,:,1) = 0d0     
+            endwhere
+            if (exp(sig**2) - 1 > 0.1) then
+               nueffrad(:,:,1) = exp(sig**2) - 1
+            else
+               nueffrad(:,:,1) = 0.1    
             endif
-         end do !iaer=1,naerkind.
-         if (haze_radproffix) then
+            ! Fractal aerosols
+            sig = 0.35
+            where (mp2m_rc_fra(:,:) > 1e-10)
+               reffrad(:,:,2) = mp2m_rc_fra(:,:) * exp(5.*sig**2 / 2.)
+               elsewhere
+               reffrad(:,:,2) = 0d0
+               endwhere
+            if (exp(sig**2) - 1 > 0.1) then
+               nueffrad(:,:,2) = exp(sig**2) - 1
+            else
+               nueffrad(:,:,2) = 0.1    
+            endif
+           
+         else
+            do iaer=1,naerkind
+               if ((iaer.eq.iaero_haze)) then
+               call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer),nueffrad(1,1,iaer))
+               endif
+            end do ! iaer = 1, naerkind.
             if (haze_radproffix) then
-                call haze_reffrad_fix(ngrid,nlayer,zzlay,&
-                    reffrad,nueffrad)
-            endif
-
-            print*, 'haze_radproffix=T : fixed profile for haze rad'
-         else
-            print*,'reffrad haze:',reffrad(1,1,iaero_haze)
-            print*,'nueff haze',nueffrad(1,1,iaero_haze)
-         endif
-      endif
-
+               call haze_reffrad_fix(ngrid,nlayer,zzlay,reffrad,nueffrad)
+               if (is_master) print*, 'haze_radproffix=T : fixed profile for haze rad'
+            else
+               if (is_master) print*,'reffrad haze:',reffrad(1,1,iaero_haze)
+               if (is_master) print*,'nueff haze',nueffrad(1,1,iaero_haze)
+            endif ! end haze_radproffix
+         endif ! end callmufi
+      endif ! end radiative haze
 
       ! How much light do we get ?
@@ -489,5 +511,5 @@
 
          ! Get aerosol optical depths.
-         call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
+         call aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq,dtau_aer,      &
             reffrad,nueffrad,QREFvis3d,QREFir3d,                             &
             tau_col)
@@ -665,7 +687,7 @@
                pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
                        (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
-               ! As 'aerosol' is at reference (visible) wavelenght we scale it as
+               ! As 'dtau_aer' is at reference (visible) wavelenght we scale it as
                ! it will be multplied by qxi/v in optci/v
-               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
+               temp=dtau_aer(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)
@@ -986,5 +1008,5 @@
             print*,'fluxtop_dn=',fluxtop_dn(ig)
             print*,'acosz=',acosz
-            print*,'aerosol=',aerosol(ig,:,:)
+            print*,'dtau_aer=',dtau_aer(ig,:,:)
             print*,'temp=   ',pt(ig,:)
             print*,'pplay=  ',pplay(ig,:)
Index: trunk/LMDZ.PLUTO/libf/phypluto/callcorrk_pluto_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callcorrk_pluto_mod.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callcorrk_pluto_mod.F90	(revision 3585)
@@ -7,5 +7,5 @@
     subroutine callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,   &
           albedo,emis,mu0,pplev,pplay,pt,   &
-          zzlay,tsurf,fract,dist_star,aerosol,    &
+          zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,    &
           dtlw,dtsw,fluxsurf_lw,    &
           fluxsurf_sw,fluxtop_lw,fluxtop_sw,fluxtop_dn, &
@@ -26,5 +26,5 @@
       use callkeys_mod, only: optichaze,ch4fix,cooling,methane,nlte,&
                               strobel,vmrch4_proffix,specOLR,vmrch4fix,&
-                              haze_radproffix
+                              haze_radproffix,callmufi
       use optcv_pluto_mod, only: optcv_pluto
       use optci_pluto_mod, only: optci_pluto
@@ -32,5 +32,5 @@
       use sfluxv_pluto_mod, only: sfluxv_pluto
       use mod_phys_lmdz_para, only : is_master
-
+      use mp2m_diagnostics
 
       implicit none
@@ -61,11 +61,12 @@
       INTEGER icount
       INTEGER ngrid,nlayer
-      REAL aerosol(ngrid,nlayer,naerkind) ! aerosol opacity tau
-      REAL albedo(ngrid)                    ! SW albedo
-      REAL emis(ngrid)                      ! LW emissivity
-      REAL pplay(ngrid,nlayer)            ! pres. level in GCM mid of layer
-      REAL zzlay(ngrid,nlayer)            ! altitude at the middle of the layers
-      REAL pplev(ngrid,nlayer+1)          ! pres. level at GCM layer boundaries
-
+      REAL dtau_aer(ngrid,nlayer,naerkind) ! aerosol opacity tau
+      REAL albedo(ngrid)                   ! SW albedo
+      REAL emis(ngrid)                     ! LW emissivity
+      REAL pplay(ngrid,nlayer)             ! pres. level in GCM mid of layer
+      REAL zzlay(ngrid,nlayer)             ! Altitude at the middle of the layers
+      REAL zzlev(ngrid,nlayer)             ! Altitude at the layer boundaries.
+      REAL pplev(ngrid,nlayer+1)           ! pres. level at GCM layer boundaries
+ 
       REAL pt(ngrid,nlayer)               ! air temperature (K)
       REAL tsurf(ngrid)                     ! surface temperature (K)
@@ -135,4 +136,6 @@
       INTEGER ig,l,k,nw,iaer,irad
 
+      real*8 sig
+
       real fluxtoplanet
       real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
@@ -234,32 +237,42 @@
          call setspv            ! basic visible properties
 
+         !--------------------------------------------------
+         !     Effective radius and variance of the aerosols
+         !--------------------------------------------------
          ! Radiative Hazes
          if (optichaze) then
-
-           ! AF24: TODO check duplicate suaer_corrk called from physiq_mod
-          !  print*,'optichaze: starting suaer_corrk'
-          !  call suaer_corrk       ! set up aerosol optical properties
-          !  print*,'ending suaer_corrk'
-
-           !--------------------------------------------------
-           !     Effective radius and variance of the aerosols
-           !--------------------------------------------------
-           do iaer=1,naerkind
+           if (callmufi) then
+            ! Spherical aerosols
+            sig = 0.2
+            WHERE(mp2m_rc_sph(:,:) > 1e-10)
+               reffrad(:,:,1) = mp2m_rc_sph(:,:) * exp(5.*sig**2 / 2.)
+            ELSEWHERE
+               reffrad(:,:,1) = 0d0         
+            ENDWHERE
+            nueffrad(:,:,1) = exp(sig**2) - 1
+            ! Fractal aerosols
+            sig = 0.35
+            WHERE(mp2m_rc_fra(:,:) > 1e-10)
+               reffrad(:,:,2) = mp2m_rc_fra(:,:) * exp(5.*sig**2 / 2.)
+            ELSEWHERE
+               reffrad(:,:,2) = 0d0
+            ENDWHERE
+            nueffrad(:,:,2) = exp(sig**2) - 1
+           
+           else
+            do iaer=1,naerkind
               if ((iaer.eq.iaero_haze)) then
-               call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer),    &
-                  nueffrad(1,1,iaer))
-                ! write(*,*) "Not supported yet"
-                ! STOP
+              call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer),nueffrad(1,1,iaer))
               endif
-           end do !iaer=1,naerkind.
-           if (haze_radproffix) then
-              call haze_reffrad_fix(ngrid,nlayer,zzlay, &
-                  reffrad,nueffrad)
+            end do ! iaer = 1, naerkind.
+            if (haze_radproffix) then
+              call haze_reffrad_fix(ngrid,nlayer,zzlay,reffrad,nueffrad)
               if (is_master) print*, 'haze_radproffix=T : fixed profile for haze rad'
-           else
+            else
               if (is_master) print*,'reffrad haze:',reffrad(1,1,iaero_haze)
               if (is_master) print*,'nueff haze',nueffrad(1,1,iaero_haze)
-           endif
-         endif ! radiative haze
+            endif ! end haze_radproffix
+           endif ! end callmufi
+         endif ! end radiative haze
 
          Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
@@ -305,5 +318,5 @@
 
         ! Get aerosol optical depths.
-        call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
+        call aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq,dtau_aer,      &
             reffrad,nueffrad,QREFvis3d,QREFir3d,                             &
             tau_col)
@@ -471,5 +484,5 @@
                    (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
                if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then
-                 temp=aerosol(ig,L_NLAYRAD-k,iaer)/ &
+                 temp=dtau_aer(ig,L_NLAYRAD-k,iaer)/ &
                    QREFvis3d(ig,L_NLAYRAD-k,iaer)
                else
Index: trunk/LMDZ.PLUTO/libf/phypluto/callsedim_pluto.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callsedim_pluto.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callsedim_pluto.F90	(revision 3585)
@@ -3,5 +3,4 @@
                      pq, pdqfi, pdqsed,pdqs_sed,nq,pphi)
 
-      use radinc_h, only : naerkind
       use tracer_h, only: igcm_ch4_ice,igcm_co_ice,radius,rho_q
       use comcstfi_mod, only:  g
Index: trunk/LMDZ.PLUTO/libf/phypluto/ch4cloud.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/ch4cloud.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/ch4cloud.F90	(revision 3585)
@@ -1,3 +1,3 @@
-       SUBROUTINE ch4cloud(ngrid,nlay,naersize, ptimestep,  &
+       SUBROUTINE ch4cloud(ngrid,nlay,ptimestep,  &
                      pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, &
                      pq,pdq,pdqcloud,pdqscloud,pdtcloud,  &
@@ -47,5 +47,4 @@
       real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
       real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
-      integer naersize   ! nombre de traceurs radiativement actifs (=naerkind)
       integer nq         ! nombre de traceurs
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/cocloud.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/cocloud.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/cocloud.F90	(revision 3585)
@@ -1,3 +1,3 @@
-       SUBROUTINE cocloud(ngrid,nlay,naersize, ptimestep,   &
+       SUBROUTINE cocloud(ngrid,nlay,ptimestep,   &
                      pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,   &
                      pq,pdq,pdqcloud,pdqscloud,pdtcloud,  &
@@ -44,5 +44,4 @@
       real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
       real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
-      integer naersize   ! nombre de traceurs radiativement actifs (=naerkind)
       integer nq         ! nombre de traceurs
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90	(revision 3585)
@@ -6,5 +6,4 @@
       pdqc,pdicen2)
 
-  use radinc_h, only : naerkind
   use comgeomfi_h
   use comcstfi_mod, only: g, r, cpp, pi
Index: trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90	(revision 3585)
@@ -19,4 +19,6 @@
       character(LEN=100),save :: config_mufi ='datagcm/microphysics/config.cfg'
 !$OMP THREADPRIVATE(config_mufi)
+      character(len=300),save :: aersprop_file
+      character(len=300),save :: aerfprop_file
 
       ! surfdir stores planetary topography, albedo, etc. (surface.nc files)
Index: trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/kcm1d.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/kcm1d.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/kcm1d.F90	(revision 3585)
@@ -78,5 +78,5 @@
 
   real dTstrat
-  real,allocatable :: aerosol(:,:) ! aerosol tau (kg/kg)
+  real,allocatable :: dtau_aer(:,:) ! aerosol tau (kg/kg)
   real OLR_nu(1,L_NSPECTI)
   real OSR_nu(1,L_NSPECTV)
@@ -320,9 +320,9 @@
   !write(*,*) 1,llm,nq,0,86400.0,1,1.0,latitude,longitude,cell_area,rad,g,r,cpp
 
-  ! initialise naerkind (from callphys.def) and allocate aerosol(:,:)
+  ! initialise naerkind (from callphys.def) and allocate dtau_aer(:,:)
   naerkind=0 !default
   call getin("naerkind",naerkind)
-  allocate(aerosol(llm,naerkind))
-  aerosol(:,:)=0
+  allocate(dtau_aer(llm,naerkind))
+  dtau_aer(:,:)=0
 
   do iq=1,nq
@@ -374,5 +374,5 @@
           albedo_wv,albedo_equivalent,                    &
           emis,mu0,plev,play,temp,                        &
-          tsurf,fract,dist_star,aerosol,muvar,            &
+          tsurf,fract,dist_star,dtau_aer,muvar,            &
           dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,              & 
           fluxsurfabs_sw,fluxtop_lw,                      &
@@ -385,5 +385,5 @@
      !     albedo_wv,'D',albedo_equivalent,'E',                    &
      !     emis,'F',mu0,'G',plev,'H',play,'I',temp,'J',                        &
-     !     tsurf,'K',fract,'L',dist_star,'M',aerosol,'N',muvar,'O',            &
+     !     tsurf,'K',fract,'L',dist_star,'M',dtau_aer,'N',muvar,'O',            &
      !     dtlw,'P',dtsw,'Q',fluxsurf_lw,'R',fluxsurf_sw,'S',              & 
      !     fluxsurfabs_sw,'T',fluxtop_lw,'U',                      &
@@ -426,5 +426,5 @@
   call callcorrk(1,nlayer,q,nq,qsurf,                          &
        albedo_wv,albedo_equivalent,emis,mu0,plev,play,temp,    &
-       tsurf,fract,dist_star,aerosol,muvar,                    &
+       tsurf,fract,dist_star,dtau_aer,muvar,                    &
        dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,       &
        fluxtop_lw, fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,GSR_nu, &
@@ -437,5 +437,5 @@
      !     albedo_wv,'D',albedo_equivalent,'E',                    &
      !     emis,'F',mu0,'G',plev,'H',play,'I',temp,'J',                        &
-     !     tsurf,'K',fract,'L',dist_star,'M',aerosol,'N',muvar,'O',            &
+     !     tsurf,'K',fract,'L',dist_star,'M',dtau_aer,'N',muvar,'O',            &
      !     dtlw,'P',dtsw,'Q',fluxsurf_lw,'R',fluxsurf_sw,'S',              & 
      !     fluxsurfabs_sw,'T',fluxtop_lw,'U',                      &
Index: trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/rcm1d.F	(revision 3585)
@@ -6,5 +6,5 @@
       use mod_grid_phy_lmdz, only : regular_lonlat
       use infotrac, only: nqtot, tname
-      use tracer_h, only: noms, is_condensable
+      use tracer_h, only: noms
       use surfdat_h, only: albedodat, phisfi, dryness,
      &                     zmea, zstd, zsig, zgam, zthe,
Index: trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/hazecloud.F90	(revision 3585)
@@ -5,5 +5,4 @@
 !          zdqphot_ch4,zdqconv_prec,declin,zdqhaze_col)
 
-      use radinc_h, only : naerkind
       use comgeomfi_h
       use comcstfi_mod, only: pi, g
Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90	(revision 3585)
@@ -12,6 +12,7 @@
   use radinc_h, only: ini_radinc_h, naerkind
   use radcommon_h, only: ini_radcommon_h
-  use radii_mod, only: radfixed, Nmix_n2
-  use datafile_mod, only: datadir,config_mufi,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file
+  use radii_mod, only: radfixed
+  use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file, &
+                          config_mufi, aersprop_file, aerfprop_file
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
   use comgeomfi_h, only: totarea, totarea_planet
@@ -687,4 +688,14 @@
      call getin_p("config_mufi",config_mufi)
      if (is_master) write(*,*)" config_mufi = ",config_mufi
+
+     if (is_master) write(*,*) "Spherical aerosol optical properties datafile"
+     aersprop_file="optprop_rannou_r2-200nm_nu003.dat"  ! default file
+     call getin_p("aersprop_file",aersprop_file)
+     if (is_master) write(*,*) trim(rname)//" aersprop_file = ",trim(aersprop_file)
+
+     if (is_master) write(*,*) "Fractal aerosol optical properties datafile"
+     aerfprop_file="optprop_rannou_fractal_r010nm_N1_1e4_d2.dat"  ! default file
+     call getin_p("aerfprop_file",aerfprop_file)
+     if (is_master) write(*,*) trim(rname)//" aerfprop_file = ",trim(aerfprop_file)
  
      if (is_master) write(*,*) "Use haze production from CH4 photolysis or production rate?"
@@ -1285,4 +1296,19 @@
        call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
      endif
+     if ((callmufi).and.(haze)) then
+       call abort_physic(rname, 'if haze are on, microphysics should be desactivated!', 1)
+     endif
+     if ((haze).and.(naerkind.gt.1)) then
+      call abort_physic(rname, 'if haze are on, naerkind must be equal to 1!', 1)
+     endif
+     if ((callmufi).and..not.(naerkind.gt.1)) then
+      call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1)
+     endif
+     if (.not.(callmufi.or.haze).and.(optichaze)) then
+      call abort_physic(rname, 'if microphysics and haze are off, optichaze must be desactivated!', 1)
+     endif
+     if ((callmufi.and.call_haze_prod_pCH4).and..not.(methane)) then
+      call abort_physic(rname, 'if haze production from CH4 photolysis is on, methane must be activated!', 1)
+     endif
      if (haze_proffix.and.sedimentation) then
          call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
Index: trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90	(revision 3585)
@@ -123,12 +123,4 @@
          ALLOCATE(is_recomb_qotf(nqtot))
        ENDIF
-       IF (.NOT. allocated(is_condensable)) allocate(is_condensable(nq)) !LT
-       IF (.NOT. allocated(is_rgcs)) allocate(is_rgcs(nq)) !LT
-       IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq))
-       IF (.NOT. allocated(constants_delta_gasH)) allocate(constants_delta_gasH(nq))
-       IF (.NOT. allocated(constants_Tref)) allocate(constants_Tref(nq))
-       IF (.NOT. allocated(constants_Pref)) allocate(constants_Pref(nq))
-       IF (.NOT. allocated(constants_RLVTT_generic)) allocate(constants_RLVTT_generic(nq))
-       IF (.NOT. allocated(constants_RCPV_generic)) allocate(constants_RCPV_generic(nq))
 
        !! initialization
@@ -148,15 +140,4 @@
        radius(:)=0.
        qext(:)=0.
-
-       ! For condensable tracers, by Lucas Teinturier and Noé Clément (2022)
-
-       is_condensable(:)= 0
-       is_rgcs(:) = 0
-       constants_mass(:)=0
-       constants_delta_gasH(:)=0
-       constants_Tref(:)=0
-       constants_Pref(:)=0
-       constants_RLVTT_generic(:)=0
-       constants_RCPV_generic(:)=0
 
        rho_q(:) = 0. !need to be init here if we want to read it from modern traceur with get_tracdat
@@ -403,38 +384,7 @@
       if (is_master) close(407)
 
-      ! Get specific data of condensable tracers
-      do iq=1,nq
-        if((is_condensable(iq)==1)) then
-                write(*,*) "There is a specie which is condensable, for generic condensation : ", noms(iq)
-                write(*,*) 'looking specie parameters for : ',noms(iq)(1:len(trim(noms(iq)))-4)
-                ! call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4))
-                ! constants_mass(iq)=m
-                ! constants_delta_gasH(iq)=delta_gasH
-                ! constants_Tref(iq)=Tref
-                ! constants_Pref(iq)=Pref
-                ! constants_RLVTT_generic(iq)=RLVTT_generic
-                ! constants_RCPV_generic(iq)=RCPV_generic
-        else
-                write(*,*) "This tracer is not condensable, for generic condensation :  : ", noms(iq)
-                write(*,*) "We keep condensable constants at zero"
-        endif !(is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))
-      enddo ! iq=1,nq
-
       ! Calculate number of species in the chemistry
       nesp = sum(is_chim)
       write(*,*) 'Number of species in the chemistry nesp = ',nesp
-
-      ! Calculate number of generic tracers
-      ngt = sum(is_condensable)
-      write(*,*) 'Number of generic tracer is  ngt = ',ngt
-
-      ! Calculate number of radiative generic condensable species
-      n_rgcs = sum(is_rgcs)
-      write(*,*)'Number of Radiative Generic Condensable Species is n_rgcs = ',n_rgcs
-      if (n_rgcs> ngt/2) then
-        write(*,*) 'You have more Radiative Generic Condensable Species than Generic Condensable Species'
-        write(*,*)'This is not possible: check your Modern traceur.def'
-        call abort_physic("initracer, issue with # of RGCS and GCS")
-      endif
 
       ! Calculate number of microphysical tracer
@@ -461,19 +411,25 @@
       write(*,*) "lw_n2 = ", lw_n2
 
-      if (haze) then
+      if (callmufi) then
+        if (optichaze) then
+          iaero_haze = 2
+          write(*,*) 'Microphysical moment model'
+          write(*,*) '--- number of haze aerosol = ', iaero_haze
+        endif ! end optichaze
+      
+      elseif (haze) then
         ! the sedimentation radius remains radius(igcm_haze)
         if (fractal) then
-           nmono=nb_monomer
+          nmono=nb_monomer
         else
-           nmono=1
+          nmono=1
         endif ! end fractal
-
-        ia=0
+        
+        ia = 0
         if (optichaze) then
-          ia=ia+1
-          iaero_haze=ia
+          ia = ia + 1
+          iaero_haze = ia
           write(*,*) '--- number of haze aerosol = ', iaero_haze
-
-          block=0  ! Only one type of haze is active : the first one set in traceur.def
+          block = 0 ! Only one type of haze is active : the first one set in traceur.def
           do iq=1,nq
             tracername=noms(iq)
@@ -484,8 +440,8 @@
               write(*,*) "i_haze=",i_haze
               write(*,*) "Careful: if you set many haze traceurs in &
-                traceur.def,only ",tracername," will be radiatively active &
-                (first one in traceur.def)"
+              traceur.def,only ",tracername," will be radiatively active &
+              (first one in traceur.def)"
             endif
-          enddo
+            enddo
         endif ! end optichaze
       endif ! end callmufi or haze
@@ -602,14 +558,4 @@
               is_recomb_qotf(iq)
           end if
-          !option is_condensable (LT)
-          if (index(tracline,'is_condensable=') /=0) then
-            read(tracline(index(tracline,'is_condensable=') &
-              +len('is_condensable='):),*) is_condensable(iq)
-            write(*,*) ' Parameter value (traceur.def) :'// &
-              ' is_condensable=', is_condensable(iq)
-          else
-              write(*,*) ' Parameter value (default)     :'// &
-              ' is_condensable=', is_condensable(iq)
-          endif
           !option radius
           if (index(tracline,'radius=') .ne. 0) then
@@ -632,14 +578,4 @@
               ' rho=', rho_q(iq)
           endif
-          !option is_rgcs
-          if (index(tracline,'is_rgcs') .ne. 0) then
-            read(tracline(index(tracline,'is_rgcs=') &
-              +len('is_rgcs='):),*) is_rgcs(iq)
-            write(*,*)'Parameter value (traceur.def) :'// &
-              'is_rgcs=',is_rgcs(iq)
-          else
-            write(*,*)'Parameter value (default) : '// &
-              'is_rgcs = ',is_rgcs(iq)
-          endif
       end subroutine get_tracdat
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3585)
@@ -21,5 +21,4 @@
       use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
       use suaer_corrk_mod, only: suaer_corrk
-      use radii_mod, only: su_aer_radii,haze_reffrad_fix
       use aerosol_mod, only: i_haze, haze_prof
       use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
@@ -35,5 +34,5 @@
                           igcm_co_gas,igcm_co_ice,igcm_prec_haze,lw_n2,lw_ch4,lw_co,&
                           alpha_lift, alpha_devil, qextrhor, &
-                          nesp, is_chim, is_condensable
+                          nesp, is_chim
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyetat0_mod, only: phyetat0,tab_cntrl_mod
@@ -286,6 +285,6 @@
 !     for the "naerkind" optically active aerosols:
 
-      real,save,allocatable :: aerosol(:,:,:) ! Aerosols
-!$OMP THREADPRIVATE(aerosol)
+      real,save,allocatable :: dtau_aer(:,:,:) ! Aerosols
+!$OMP THREADPRIVATE(dtau_aer)
       real zh(ngrid,nlayer)               ! Potential temperature (K).
       real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
@@ -464,7 +463,4 @@
       real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
 
-      real,save,allocatable :: reffcol(:,:)
-!$OMP THREADPRIVATE(reffcol)
-
 !  Non-oro GW tendencies
       REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
@@ -535,6 +531,5 @@
          ! allocate related local arrays
          ! (need be allocated instead of automatic because of "naerkind")
-         allocate(aerosol(ngrid,nlayer,naerkind))
-         allocate(reffcol(ngrid,naerkind))
+         allocate(dtau_aer(ngrid,nlayer,naerkind))
 
 #ifdef CPP_XIOS
@@ -1022,6 +1017,6 @@
                if (oldplutocorrk) then
                   call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,          &
-                               albedo(:,1),emis,mu0,pplev,pplay,pt,                   &
-                               zzlay,tsurf,fract,dist_star,aerosol,              &
+                               albedo(:,1),emis,mu0,pplev,pplay,pt,              &
+                               zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,       &
                                zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,   &
                                fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
@@ -1034,5 +1029,5 @@
                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,  &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
-                              zzlay,tsurf,fract,dist_star,aerosol,muvar,                &
+                              zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar,   &
                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
                               fluxsurfabs_sw,fluxtop_lw,                          &
@@ -1436,5 +1431,5 @@
          zdqch4cloud(:,:,:)=0.
          if ((methane).and.(metcloud).and.(.not.fast)) THEN
-            call ch4cloud(ngrid,nlayer,naerkind,ptimestep,  &
+            call ch4cloud(ngrid,nlayer,ptimestep,  &
                       pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,  &
                       pq,pdq,zdqch4cloud,zdqsch4cloud,zdtch4cloud,   &
@@ -1469,5 +1464,5 @@
          zdqcocloud(:,:,:)=0.
          IF ((carbox).and.(monoxcloud).and.(.not.fast)) THEN
-            call cocloud(ngrid,nlayer,naerkind,ptimestep,   &
+            call cocloud(ngrid,nlayer,ptimestep,   &
                       pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,  &
                       pq,pdq,zdqcocloud,zdqscocloud,zdtcocloud,   &
@@ -1845,5 +1840,5 @@
                print*,'fluxtop_dn=',fluxtop_dn(ig)
                print*,'tau_col=',tau_col(ig)
-               print*,'aerosol=',aerosol(ig,:,:)
+               print*,'dtau_aer=',dtau_aer(ig,:,:)
                print*,'temp=   ',pt(ig,:)
                print*,'pplay=  ',pplay(ig,:)
Index: trunk/LMDZ.PLUTO/libf/phypluto/prodhaze.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/prodhaze.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/prodhaze.F90	(revision 3585)
@@ -3,5 +3,4 @@
           flym_bot,flym_sol_bot,flym_ipm_bot,flym_sol,flym_ipm)
 
-      use radinc_h, only : naerkind
       use comgeomfi_h
       use comcstfi_mod, only: pi, g
Index: trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/radii_mod.F90	(revision 3585)
@@ -4,8 +4,4 @@
 !  module to centralize the radii calculations for aerosols
 !==================================================================
-
-!     N2 cloud properties (initialized in inifis)
-      real,save :: Nmix_n2 ! Number mixing ratio of N2 ice particles
-!$OMP THREADPRIVATE(Nmix_n2)
 
       ! flag to specify if we assume a constant fixed radius for particles
@@ -30,5 +26,5 @@
       use datafile_mod, only: hazerad_file
       use aerosol_mod, only: iaero_haze, i_haze
-      use tracer_h, only: radius, nqtot, is_rgcs, nmono
+      use tracer_h, only: radius, nqtot, nmono
       Implicit none
 
@@ -67,5 +63,5 @@
       USE tracer_h, only:rho_n2,nmono
       use comcstfi_mod, only: pi
-      use aerosol_mod, only: iaero_haze, i_haze
+      use aerosol_mod, only: iaero_haze
       use datafile_mod
       Implicit none
Index: trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90	(revision 3585)
@@ -6,5 +6,4 @@
       use comgeomfi_h
       use geometry_mod, only: latitude, longitude, cell_area
-      use radinc_h, only : naerkind
       use tracer_h, only: igcm_n2,igcm_ch4_ice,igcm_co_ice
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90	(revision 3585)
@@ -2,5 +2,4 @@
 
       use geometry_mod, only: latitude, longitude, cell_area
-      use radinc_h, only : naerkind
       use surfdat_h, only: phisfi
       use tracer_h, only: igcm_n2
Index: trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90	(revision 3585)
@@ -10,5 +10,6 @@
       use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
       use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
-      use datafile_mod, only: datadir, aerdir, hazeprop_file
+      use datafile_mod, only: datadir, aerdir, &
+                              hazeprop_file, aersprop_file, aerfprop_file
 
       ! outputs
@@ -17,5 +18,5 @@
       use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
       use aerosol_mod, only: iaero_haze
-      use callkeys_mod, only: tplanet
+      use callkeys_mod, only: tplanet, callmufi
       use tracer_h, only: noms
 
@@ -129,28 +130,50 @@
       allocate(file_id(naerkind,2))
 
-      do iaer=1,naerkind
-         !     naerkind=1: haze
-         if (iaer.eq.1) then
-
-         ! Only one table of optical properties :
-         if (is_master) write(*,*)'Suaer haze optical properties, using: ', &
-                           TRIM(hazeprop_file)
-
-         ! visible
-         file_id(iaer,1)=TRIM(hazeprop_file)
-         ! infrared
-         file_id(iaer,2)=file_id(iaer,1)
-
-         lamrefvis(iaer)=0.607E-6   ! reference wavelength for opacity vis pivot wavelength Cheng et al 2017
-         lamrefir(iaer)=2.E-6   !  reference wavelength for opacity IR (in the LEISA range)
-
+
+      if (callmufi) then
+         if (is_master) then
+            write(*,*)'Suaer spherical aerosols optical properties, using: ', &
+                       TRIM(aersprop_file)
+            write(*,*)'Suaer fractal aerosols optical properties, using: ', &
+                       TRIM(aerfprop_file) 
          endif
-      enddo
-
-      ! IF (naerkind .gt. 1) THEN ! AF24: ?
-      !    print*,'naerkind = ',naerkind
-      !    print*,'but we only have data for 1 type, exiting.'
-      !    call abort
-      ! ENDIF
+         ! Visible
+         file_id(1,1)=TRIM(aersprop_file)
+         file_id(2,1)=TRIM(aerfprop_file)
+         ! Infrared
+         file_id(1,2)=file_id(1,1)
+         file_id(2,2)=file_id(2,1)
+
+         do iaer=1,naerkind
+            lamrefvis(iaer) = 0.607E-6 ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
+            lamrefir(iaer)  = 2.E-6    ! Reference wavelength for opacity IR (in the LEISA range)
+         enddo ! end iaer=1,naerkind
+         
+      else 
+         do iaer=1,naerkind
+            if (iaer.eq.1) then
+               ! Only one table of optical properties :
+               if (is_master) write(*,*)'Suaer haze optical properties, using: ', &
+                                       TRIM(hazeprop_file)
+               ! Visible
+               file_id(iaer,1)=TRIM(hazeprop_file)
+               ! Infrared
+               file_id(iaer,2)=file_id(iaer,1)
+
+               lamrefvis(iaer) = 0.607E-6 ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
+               lamrefir(iaer)  = 2.E-6    ! Reference wavelength for opacity IR (in the LEISA range)
+            endif ! end iaer.eq.1
+
+            if (naerkind.gt.1) then
+               ! Visible
+               file_id(iaer,1)=file_id(1,1)
+               ! Infrared
+               file_id(iaer,2)=file_id(1,1)
+
+               lamrefvis(iaer) = lamrefvis(1) ! Reference wavelength for opacity vis pivot wavelength Cheng et al 2017
+               lamrefir(iaer)  = lamrefir(1)  ! Reference wavelength for opacity IR (in the LEISA range)
+            endif ! end naerkind.gt.1
+         enddo ! end iaer=1,naerkind
+      endif ! end callmufi
 
 !------------------------------------------------------------------
Index: trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90	(revision 3583)
+++ trunk/LMDZ.PLUTO/libf/phypluto/tracer_h.F90	(revision 3585)
@@ -13,7 +13,5 @@
        integer, save :: nqtot  ! total number of tracers
        integer, save :: nesp   ! number of species in the chemistry
-       integer, save :: ngt    ! number of generic tracers
-       integer, save :: n_rgcs ! number of Radiative Generic Condensable Species
-!$OMP THREADPRIVATE(nqtot,nesp,ngt,n_rgcs)
+!$OMP THREADPRIVATE(nqtot,nesp)
 
        logical :: moderntracdef=.false. ! Standard or modern traceur.def
@@ -54,17 +52,5 @@
        integer, save, allocatable :: is_recomb_qotf(:) ! 1 if tracer recombination is done on-the-fly, else 0 (if 1, must have is_recomb_qset=0)
 !$OMP THREADPRIVATE(is_recomb,is_recomb_qset,is_recomb_qotf)
-       integer, save, allocatable :: is_condensable(:)      ! 1 if tracer is generic, else 0 (added LT)
-       integer,save,allocatable :: is_rgcs(:)               ! 1 if tracer is a radiative generic condensable specie, else 0 (added LT 2022)
-!$OMP THREADPRIVATE(is_condensable,is_rgcs)   !also added by LT
        ! Lists of constants for condensable tracers
-       real, save, allocatable :: constants_mass(:)                 ! molecular mass of the specie (g/mol)
-       real, save, allocatable :: constants_delta_gasH(:)           ! Enthalpy of vaporization (J/mol)
-       real, save, allocatable :: constants_Tref(:)                 ! Ref temperature for Clausis-Clapeyron (K)
-       real, save, allocatable :: constants_Pref(:)                 ! Reference pressure for Clausius Clapeyron (Pa)
-       real, save, allocatable :: constants_RLVTT_generic(:)        ! Latent heat of vaporization (J/kg)
-       real, save, allocatable :: constants_RCPV_generic(:)         ! specific heat capacity of the tracer vapor at Tref
-!$OMP THREADPRIVATE(constants_mass,constants_delta_gasH,constants_Tref)
-!$OMP THREADPRIVATE(constants_Pref)
-!$OMP THREADPRIVATE(constants_RLVTT_generic,constants_RCPV_generic)
 
 ! tracer indexes: these are initialized in initracer and should be 0 if the
