Index: trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90	(revision 3940)
+++ trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90	(revision 3947)
@@ -144,5 +144,5 @@
   DTAUKV(:,:,:)    = 0.0
   dtauv_aer(:,:,:) = 0.0
-  wbarv_aer(:,:,:) = 0.0
+  wbarv_aer(:,:,:) = 1.0
 
   if(callmufi) then
@@ -353,4 +353,5 @@
       ! Aerosol absorption correction --> impact on single scattering albedo.
       if (callmufi) then
+         ! Spherical aerosols
          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) / &
@@ -363,4 +364,5 @@
             btemp(L,NW) = TAUAEROLK(K,NW,1) + TAUAEROLK(K+1,NW,1)
          endif
+         ! Fractal aerosols
          if ((TAEROS(K,NW,2).gt.0.d0) .and. TAEROS(K+1,NW,2).gt.0.d0) then
             btemp(L,NW) = btemp(L,NW) + &
@@ -390,4 +392,5 @@
      ! Aerosol absorption correction --> impact on single scattering albedo.
      if (callmufi) then
+      ! Spherical aerosols
       if (TAEROS(K,NW,1).gt.0.d0) then
          btemp(L,NW) = TAUAEROLK(K,NW,1) / &
@@ -397,4 +400,5 @@
          btemp(L,NW) = TAUAEROLK(K,NW,1)
       endif
+      ! Fractal aerosols
       if (TAEROS(K,NW,2).gt.0.d0) then
          btemp(L,NW) = btemp(L,NW) + &
@@ -442,18 +446,36 @@
       DTAUV_AER(L,nw,iaer) = TAEROS(K,nw,iaer) + TAEROS(K+1,nw,iaer)
 
-      wbarv_prime = (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)) / &
-                    (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer) + Fabs_aer(iaer)*(1.-QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)))
+      IF (QXVAER(K,nw,iaer) > 0.0D0) THEN
+         wbarv_prime = (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)) / &
+                       (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer) + Fabs_aer(iaer)*(1.-QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)))
+      ELSE
+         wbarv_prime = 1.0
+      ENDIF
       WBARV_AER(L,nw,iaer) = wbarv_prime * TAEROS(K,nw,iaer)
-      wbarv_prime = (QSVAER(K+1,nw,iaer)/QXVAER(K+1,nw,iaer)) / &
-                    (QSVAER(K+1,nw,iaer)/QXVAER(K+1,nw,iaer) + Fabs_aer(iaer)*(1.-QSVAER(K+1,nw,iaer)/QXVAER(K+1,nw,iaer)))
+      IF (QXVAER(K+1,nw,iaer) > 0.0D0) THEN
+         wbarv_prime = (QSVAER(K+1,nw,iaer)/QXVAER(K+1,nw,iaer)) / &
+                       (QSVAER(K+1,nw,iaer)/QXVAER(K+1,nw,iaer) + Fabs_aer(iaer)*(1.-QSVAER(K+1,nw,iaer)/QXVAER(K+1,nw,iaer)))
+      ELSE
+         wbarv_prime = 1.0
+      ENDIF
       WBARV_AER(L,nw,iaer) = WBARV_AER(L,nw,iaer) + (wbarv_prime * TAEROS(K+1,nw,iaer))
-      WBARV_AER(L,nw,iaer) = WBARV_AER(L,nw,iaer) / DTAUV_AER(L,nw,iaer)
-      END DO ! L vertical loop
-      ! Last level
-      L              = L_NLAYRAD
-      K              = 2*L+1
-      DTAUV_AER(L,nw,iaer) = TAEROS(K,nw,iaer)
-      WBARV_AER(L,nw,iaer) = (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)) / &
-                             (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer) + Fabs_aer(iaer)*(1.-QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)))
+      IF (DTAUV_AER(L,nw,iaer) > 0.0D0) THEN
+         WBARV_AER(L,nw,iaer) = WBARV_AER(L,nw,iaer) / DTAUV_AER(L,nw,iaer)
+      ELSE
+         WBARV_AER(L,nw,iaer) = 1.0
+      ENDIF
+     END DO ! L vertical loop
+
+     ! Last level
+     !-----------
+     L = L_NLAYRAD
+     K = 2*L+1
+     DTAUV_AER(L,nw,iaer) = TAEROS(K,nw,iaer)
+     IF (QXVAER(K,nw,iaer) > 0.0D0) THEN
+        WBARV_AER(L,nw,iaer) = (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)) / &
+                               (QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer) + Fabs_aer(iaer)*(1.-QSVAER(K,nw,iaer)/QXVAER(K,nw,iaer)))
+     ELSE
+       WBARV_AER(L,nw,iaer) = 1.0
+     ENDIF
    END DO ! nw spectral loop
   END DO ! iaer Gauss loop