Index: trunk/LMDZ.PLUTO/libf/phypluto/optci.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/optci.F90	(revision 3880)
+++ trunk/LMDZ.PLUTO/libf/phypluto/optci.F90	(revision 3881)
@@ -15,5 +15,5 @@
                      igas_CH4, igas_N2
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody
+  use callkeys_mod, only: kastprof,continuum,graybody,callmufi
   use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
   use tpindex_mod, only: tpindex
@@ -94,4 +94,7 @@
   !real*8 rho !! see test below
 
+  ! Variables for aerosol absorption
+  real*8 Fabs_aer(NAERKIND)
+
   integer igas, jgas
 
@@ -140,4 +143,11 @@
   lkcoef(:,:)   = 0.0
 
+  if(callmufi) then
+   Fabs_aer(1) = 1.2
+   Fabs_aer(2) = 1.3
+  else
+   Fabs_aer(:) = 1.0
+  endif
+
   do K=2,L_LEVELS
      DPR(k) = PLEV(K)-PLEV(K-1)
@@ -176,7 +186,12 @@
         do K=2,L_LEVELS
            TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
-        end do                    ! levels
-     END DO
-  end do
+           ! Aerosol absorption correction --> impact on opacity.
+           if (callmufi .and. (TAEROS(K,NW,IAER).gt.0.d0)) then
+            TAEROS(K,NW,IAER) = TAEROS(K,NW,IAER) * &
+                                ((QSIAER(K,NW,IAER)/QXIAER(K,NW,IAER)) + Fabs_aer(IAER)*(1.-(QSIAER(K,NW,IAER)/QXIAER(K,NW,IAER))))
+           endif
+        end do ! L_LEVELS
+     END DO ! L_NSPECTI
+  end do ! naerkind
 
   do NW=1,L_NSPECTI
@@ -322,17 +337,61 @@
     ENDDO
   end do
-  
-  DO NW=1,L_NSPECTI
-     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))
+
+  DO NW = 1, L_NSPECTI
+     DO L = 1, L_NLAYRAD-1
+        K = 2*L + 1
+        ! Aerosol absorption correction --> impact on single scattering albedo.
+        if (callmufi) then
+           if ((TAEROS(K,NW,1).gt.0.d0) .and. TAEROS(K+1,NW,1).gt.0.d0) then
+              btemp(L,NW) = TAUAEROLK(K,NW,1) / &
+                          (QSIAER(K,NW,1)/QXIAER(K,NW,1) + &
+                             Fabs_aer(1)*(1.-QSIAER(K,NW,1)/QXIAER(K,NW,1))) + &
+                          TAUAEROLK(K+1,NW,1) / &
+                          (QSIAER(K+1,NW,1)/QXIAER(K+1,NW,1) + &
+                             Fabs_aer(1)*(1.-QSIAER(K+1,NW,1)/QXIAER(K+1,NW,1)))
+           else
+              btemp(L,NW) = TAUAEROLK(K,NW,1) + TAUAEROLK(K+1,NW,1)
+           endif
+           if ((TAEROS(K,NW,2).gt.0.d0) .and. TAEROS(K+1,NW,2).gt.0.d0) then
+              btemp(L,NW) = btemp(L,NW) + &
+                          (TAUAEROLK(K,NW,2) / &
+                          (QSIAER(K,NW,2)/QXIAER(K,NW,2) + &
+                             Fabs_aer(2)*(1.-QSIAER(K,NW,2)/QXIAER(K,NW,2)))) + &
+                          (TAUAEROLK(K+1,NW,2) / &
+                          (QSIAER(K+1,NW,2)/QXIAER(K+1,NW,2) + &
+                             Fabs_aer(2)*(1.-QSIAER(K+1,NW,2)/QXIAER(K+1,NW,2))))
+           else
+              btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2) + TAUAEROLK(K+1,NW,2)
+           endif
+        else
+           btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
+        endif ! callmufi
      END DO ! L vertical loop
      
      ! Last level
-     L           = L_NLAYRAD
-     K           = 2*L+1    
-     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
-     
-  END DO                    ! NW spectral loop
+     !-----------
+     L = L_NLAYRAD
+     K = 2*L+1
+     ! Aerosol absorption correction --> impact on single scattering albedo.
+     if (callmufi) then
+      if (TAEROS(K,NW,1).gt.0.d0) then
+         btemp(L,NW) = TAUAEROLK(K,NW,1) / &
+                       (QSIAER(K,NW,1)/QXIAER(K,NW,1) + &
+                        Fabs_aer(1)*(1.-QSIAER(K,NW,1)/QXIAER(K,NW,1)))
+      else
+         btemp(L,NW) = TAUAEROLK(K,NW,1)
+      endif
+      if (TAEROS(K,NW,2).gt.0.d0) then
+         btemp(L,NW) = btemp(L,NW) + &
+                       (TAUAEROLK(K,NW,2) / &
+                       (QSIAER(K,NW,2)/QXIAER(K,NW,2) + &
+                        Fabs_aer(2)*(1.-QSIAER(K,NW,2)/QXIAER(K,NW,2))))
+      else
+         btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2)
+      endif
+     else
+      btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
+     endif ! callmufi
+  END DO ! NW spectral loop
   
 
@@ -366,5 +425,5 @@
      
      ! Last level
-     
+     !-----------
      L              = L_NLAYRAD
      K              = 2*L+1
Index: trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90	(revision 3880)
+++ trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90	(revision 3881)
@@ -14,5 +14,5 @@
                      igas_CH4, igas_N2
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis
+  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,callmufi
   use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
   use tpindex_mod, only: tpindex
@@ -100,4 +100,6 @@
   real*8 dz(L_LEVELS)
 
+  ! Variables for aerosol absorption
+  real*8 Fabs_aer(NAERKIND)
 
   integer igas, jgas
@@ -137,4 +139,11 @@
   lkcoef(:,:)   = 0.0
   DTAUKV(:,:,:) = 0.0
+
+  if(callmufi) then
+   Fabs_aer(1) = 0.7
+   Fabs_aer(2) = 2.5
+  else
+   Fabs_aer(:) = 1.0
+  endif
 
   do K=2,L_LEVELS
@@ -160,5 +169,5 @@
 
   ! Spectral dependance of aerosol absorption
-            !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
+       !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
 	    !   but visible does not handle very well diffusion in first layer.
 	    !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
@@ -170,8 +179,13 @@
         do K=5,L_LEVELS
            TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
-        end do                    ! levels
-     end do
-  end do
-
+           ! Aerosol absorption correction --> impact on opacity.
+           if (callmufi .and. (TAEROS(K,NW,IAER).gt.0.d0)) then
+            TAEROS(K,NW,IAER) = TAEROS(K,NW,IAER) * &
+                                ((QSVAER(K,NW,IAER)/QXVAER(K,NW,IAER)) + Fabs_aer(IAER)*(1.-(QSVAER(K,NW,IAER)/QXVAER(K,NW,IAER))))
+           endif
+        end do ! L_LEVELS
+     end do ! L_NSPECTV
+  end do ! naerkind
+  
   ! Rayleigh scattering
   do NW=1,L_NSPECTV
@@ -328,23 +342,69 @@
   DO NW=1,L_NSPECTV
      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))
-        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
-	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)
-	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
+        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))
+      ! Aerosol absorption correction --> impact on single scattering albedo.
+      if (callmufi) then
+         if ((TAEROS(K,NW,1).gt.0.d0) .and. TAEROS(K+1,NW,1).gt.0.d0) then
+            btemp(L,NW) = TAUAEROLK(K,NW,1) / &
+                          (QSVAER(K,NW,1)/QXVAER(K,NW,1) + &
+                           Fabs_aer(1)*(1.-QSVAER(K,NW,1)/QXVAER(K,NW,1))) + &
+                          TAUAEROLK(K+1,NW,1) / &
+                          (QSVAER(K+1,NW,1)/QXVAER(K+1,NW,1) + &
+                           Fabs_aer(1)*(1.-QSVAER(K+1,NW,1)/QXVAER(K+1,NW,1)))
+         else
+            btemp(L,NW) = TAUAEROLK(K,NW,1) + TAUAEROLK(K+1,NW,1)
+         endif
+         if ((TAEROS(K,NW,2).gt.0.d0) .and. TAEROS(K+1,NW,2).gt.0.d0) then
+            btemp(L,NW) = btemp(L,NW) + &
+                          (TAUAEROLK(K,NW,2) / &
+                          (QSVAER(K,NW,2)/QXVAER(K,NW,2) + &
+                           Fabs_aer(2)*(1.-QSVAER(K,NW,2)/QXVAER(K,NW,2)))) + &
+                          (TAUAEROLK(K+1,NW,2) / &
+                          (QSVAER(K+1,NW,2)/QXVAER(K+1,NW,2) + &
+                           Fabs_aer(2)*(1.-QSVAER(K+1,NW,2)/QXVAER(K+1,NW,2))))
+         else
+            btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2) + TAUAEROLK(K+1,NW,2)
+         endif
+      else
+         btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
+      endif ! callmufi
+      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)
+	   COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
      END DO ! L vertical loop
 
      ! Last level
-     L           = L_NLAYRAD
-     K           = 2*L+1
+     !-----------
+     L = L_NLAYRAD
+     K = 2*L+1
      atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
-     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
+     ! Aerosol absorption correction --> impact on single scattering albedo.
+     if (callmufi) then
+      if (TAEROS(K,NW,1).gt.0.d0) then
+         btemp(L,NW) = TAUAEROLK(K,NW,1) / &
+                       (QSVAER(K,NW,1)/QXVAER(K,NW,1) + &
+                        Fabs_aer(1)*(1.-QSVAER(K,NW,1)/QXVAER(K,NW,1)))
+      else
+         btemp(L,NW) = TAUAEROLK(K,NW,1)
+      endif
+      if (TAEROS(K,NW,2).gt.0.d0) then
+         btemp(L,NW) = btemp(L,NW) + &
+                       (TAUAEROLK(K,NW,2) / &
+                       (QSVAER(K,NW,2)/QXVAER(K,NW,2) + &
+                        Fabs_aer(2)*(1.-QSVAER(K,NW,2)/QXVAER(K,NW,2))))
+      else
+         btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2)
+      endif
+     else
+      btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
+     endif ! callmufi
      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)
      COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
-
-
-  END DO                    ! NW spectral loop
+  END DO ! NW spectral loop
 
   DO NG=1,L_NGAUSS