Index: LMDZ6/trunk/libf/phylmd/cosp/mod_modis_sim.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cosp/mod_modis_sim.F90	(revision 3233)
+++ LMDZ6/trunk/libf/phylmd/cosp/mod_modis_sim.F90	(revision 3241)
@@ -52,5 +52,5 @@
 !
 module mod_modis_sim
-  USE MOD_COSP_TYPES, only: R_UNDEF
+  USE MOD_COSP_TYPES, only: R_UNDEF, ok_debug_cosp
   implicit none
   ! ------------------------------
@@ -384,6 +384,13 @@
            liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction
     
+    real :: seuil
     ! ---------------------------------------------------
     
+    if (ok_debug_cosp) then
+       seuil=1.e-9
+    else
+       seuil=0.0
+    endif
+
     where(cloudIce(:, :) <= 0.) 
       tauLiquidFraction(:, :) = 1. 
@@ -395,6 +402,8 @@
         ! Geometic optics limit - tau as LWP/re  (proportional to LWC/re) 
         !
-        tauLiquidFraction(:, :) = (cloudWater(:, :)/waterSize(:, :)) / &
-                                  (cloudWater(:, :)/waterSize(:, :) + cloudIce(:, :)/(ice_density * iceSize(:, :)) ) 
+! Modif AI 02 2018 
+        tauLiquidFraction(:, :) = (cloudWater(:, :)/max(waterSize(:, :), seuil) ) / &
+                                  (cloudWater(:, :)/max(waterSize(:, :), seuil) + & 
+                                   cloudIce(:, :)/(ice_density * max(iceSize(:, :), seuil)) ) 
       end where
     end where
@@ -502,11 +511,23 @@
     ! Compute mean optical thickness.
     ! ########################################################################################
-    Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask,      dim = 2) / &
-                                              Cloud_Fraction_Total_Mean(1:nPoints) 
-    Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
-                                              Cloud_Fraction_Water_Mean(1:nPoints)
-    Optical_Thickness_Ice_Mean(1:nPoints)   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / &
-                                              Cloud_Fraction_Ice_Mean(1:nPoints)
-       
+    where (Cloud_Fraction_Total_Mean == 0)
+       Optical_Thickness_Total_Mean      = R_UNDEF
+    elsewhere
+       Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask,      dim = 2) / &
+                                                 Cloud_Fraction_Total_Mean(1:nPoints) 
+    endwhere
+    where (Cloud_Fraction_Water_Mean == 0)
+       Optical_Thickness_Water_Mean      = R_UNDEF
+    elsewhere
+       Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
+                                                 Cloud_Fraction_Water_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Ice_Mean == 0)
+       Optical_Thickness_Ice_Mean        = R_UNDEF
+    elsewhere
+       Optical_Thickness_Ice_Mean(1:nPoints)   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / &
+                                                 Cloud_Fraction_Ice_Mean(1:nPoints)
+    endwhere
+   
     ! ########################################################################################
     ! We take the absolute value of optical thickness here to satisfy compilers that complains 
@@ -514,21 +535,52 @@
     ! the sum. 
     ! ########################################################################################
-    Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
-         dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
-    Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
-         dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
-    Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
-         dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
-    Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
-         Cloud_Fraction_Water_Mean(1:nPoints)
-    Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask,   dim = 2) / &
-         Cloud_Fraction_Ice_Mean(1:nPoints)
-    Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
-         max(1, count(cloudMask, dim = 2))
-    Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
-         mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
-    Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,&
-         mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
-
+    where (Cloud_Fraction_Total_Mean == 0)
+       Optical_Thickness_Total_MeanLog10 = R_UNDEF
+    elsewhere
+       Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
+           dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Water_Mean == 0)
+       Optical_Thickness_Water_MeanLog10 = R_UNDEF
+    elsewhere
+       Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
+            dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Ice_Mean == 0)
+       Optical_Thickness_Ice_MeanLog10   = R_UNDEF
+    elsewhere 
+       Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
+            dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Water_Mean == 0)
+       Cloud_Particle_Size_Water_Mean    = R_UNDEF
+    elsewhere
+       Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
+            Cloud_Fraction_Water_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Ice_Mean == 0)
+       Cloud_Particle_Size_Ice_Mean      = R_UNDEF  
+    elsewhere
+       Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask,   dim = 2) / &
+            Cloud_Fraction_Ice_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Total_Mean == 0)
+       Cloud_Top_Pressure_Total_Mean     = R_UNDEF
+    elsewhere
+       Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
+            max(1, count(cloudMask, dim = 2))
+    endwhere
+    where (Cloud_Fraction_Water_Mean == 0)
+       Liquid_Water_Path_Mean            = R_UNDEF
+    elsewhere
+       Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
+            mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
+    endwhere
+    where (Cloud_Fraction_Ice_Mean == 0)
+       Ice_Water_Path_Mean               = R_UNDEF
+    elsewhere
+       Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,&
+             mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
+    endwhere
     ! ########################################################################################
     ! Normalize pixel counts to fraction.
@@ -544,21 +596,4 @@
     ! Set clear-scenes to undefined
     ! ########################################################################################
-    where (Cloud_Fraction_Total_Mean == 0)
-       Optical_Thickness_Total_Mean      = R_UNDEF
-       Optical_Thickness_Total_MeanLog10 = R_UNDEF
-       Cloud_Top_Pressure_Total_Mean     = R_UNDEF
-    endwhere
-    where (Cloud_Fraction_Water_Mean == 0)
-       Optical_Thickness_Water_Mean      = R_UNDEF
-       Optical_Thickness_Water_MeanLog10 = R_UNDEF
-       Cloud_Particle_Size_Water_Mean    = R_UNDEF
-       Liquid_Water_Path_Mean            = R_UNDEF
-    endwhere
-    where (Cloud_Fraction_Ice_Mean == 0)
-       Optical_Thickness_Ice_Mean        = R_UNDEF
-       Optical_Thickness_Ice_MeanLog10   = R_UNDEF
-       Cloud_Particle_Size_Ice_Mean      = R_UNDEF
-       Ice_Water_Path_Mean               = R_UNDEF
-    endwhere
     where (Cloud_Fraction_High_Mean == 0)  Cloud_Fraction_High_Mean = R_UNDEF
     where (Cloud_Fraction_Mid_Mean == 0)   Cloud_Fraction_Mid_Mean = R_UNDEF
