Index: /LMDZ6/trunk/libf/phylmd/cosp/cosp_output_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/cosp_output_mod.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/cosp_output_mod.F90	(revision 3241)
@@ -12,4 +12,5 @@
                             numMODISReffLiqBins, reffLIQ_binCenters
 
+     IMPLICIT NONE
 ! cosp_output_mod
       INTEGER, PRIVATE             :: i
Index: /LMDZ6/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 3241)
@@ -6,4 +6,6 @@
    USE cosp_output_mod
   
+   IMPLICIT NONE
+
    INTEGER, SAVE  :: itau_iocosp
 !$OMP THREADPRIVATE(itau_iocosp)
@@ -29,5 +31,5 @@
     USE xios, only: xios_update_calendar, xios_field_is_active
 #endif
-
+  IMPLICIT NONE  
 !!! Variables d'entree
   integer               :: itap, Nlevlmdz, Ncolumns, Npoints
@@ -45,5 +47,5 @@
 
 !!! Variables locales
-  integer               :: icl
+  integer               :: icl,k,ip
   logical               :: ok_sync
   integer               :: itau_wcosp
@@ -78,7 +80,7 @@
 ! On n'ecrit pas quand itap=1 (cosp)
 
-   if (prt_level >= 10) then
-         WRITE(lunout,*)'DO iinit=1, iinitend ',iinitend
-   endif
+!   if (prt_level >= 10) then
+!         WRITE(lunout,*)'DO iinit=1, iinitend ',iinitend
+!   endif
 
 !!#ifdef CPP_XIOS
@@ -97,58 +99,4 @@
 !!!! Sorties Calipso
  if (cfg%Llidar_sim) then
-! Pb des valeurs indefinies, on les met a 0
-! A refaire proprement
-!  do k = 1,Nlevout
-!     do ip = 1,Npoints
-!     if(stlidar%lidarcld(ip,k).eq.R_UNDEF)then
-!      stlidar%lidarcld(ip,k)=missing_val
-!     endif
-!     if(stlidar%proftemp(ip,k).eq.R_UNDEF)then !TIBO
-!      stlidar%proftemp(ip,k)=missing_val       !TIBO
-!     endif                                     !TIBO
-!     enddo
-
-!     do ii= 1,Ncolumns                               !TIBO
-!      do ip = 1,Npoints                              !TIBO
-!!       if(stlidar%profSR(ip,ii,k).eq.R_UNDEF)then    !TIBO
-!!        stlidar%profSR(ip,ii,k)=missing_val          !TIBO
-!       if(stlidar%profSR(ip,k,ii).eq.R_UNDEF)then    !TIBO2
-!        stlidar%profSR(ip,k,ii)=missing_val          !TIBO2
-!       endif                                         !TIBO
-!      enddo                                          !TIBO
-!     enddo                                           !TIBO
-!  enddo
-
-!  do ip = 1,Npoints
-!   do k = 1,Nlevlmdz
-!     if(sglidar%beta_mol(ip,k).eq.R_UNDEF)then
-!      sglidar%beta_mol(ip,k)=missing_val
-!     endif
-
-!     do ii= 1,Ncolumns
-!       if(sglidar%beta_tot(ip,ii,k).eq.R_UNDEF)then
-!        sglidar%beta_tot(ip,ii,k)=missing_val
-!       endif
-!     enddo
-
-!    enddo    !k = 1,Nlevlmdz
-!   enddo     !ip = 1,Npoints
-
-!   do k = 1,LIDAR_NCAT
-!    do ip = 1,Npoints
-!     if(stlidar%cldlayer(ip,k).eq.R_UNDEF)then
-!       stlidar%cldlayer(ip,k)=missing_val
-!     endif
-!    enddo
-!   enddo
-
-!   do k = 1,LIDAR_NTYPE                       !OPAQ
-!    do ip = 1,Npoints                         !OPAQ
-!     if(stlidar%cldtype(ip,k).eq.R_UNDEF)then !OPAQ
-!       stlidar%cldtype(ip,k)=missing_val      !OPAQ
-!     endif                                    !OPAQ
-!    enddo                                     !OPAQ
-!   enddo                                      !OPAQ
- 
 !!! AI 02 2018 
 ! Traitement missing_val
@@ -278,45 +226,4 @@
 !!! Sorties Isccp
  if (cfg%Lisccp_sim) then
-! Traitement des valeurs indefinies
-!   do ip = 1,Npoints
-!    if(isccp%totalcldarea(ip).eq.R_UNDEF)then
-!      isccp%totalcldarea(ip)=missing_val
-!    endif
-!    if(isccp%meanptop(ip).eq.R_UNDEF)then
-!      isccp%meanptop(ip)=missing_val
-!    endif
-!    if(isccp%meantaucld(ip).eq.R_UNDEF)then
-!      isccp%meantaucld(ip)=missing_val
-!    endif
-!    if(isccp%meanalbedocld(ip).eq.R_UNDEF)then
-!     isccp%meanalbedocld(ip)=missing_val
-!    endif
-!    if(isccp%meantb(ip).eq.R_UNDEF)then
-!     isccp%meantb(ip)=missing_val
-!    endif
-!    if(isccp%meantbclr(ip).eq.R_UNDEF)then
-!       isccp%meantbclr(ip)=missing_val
-!    endif
-
-!    do k=1,7
-!     do ii=1,7
-!     if(isccp%fq_isccp(ip,ii,k).eq.R_UNDEF)then
-!      isccp%fq_isccp(ip,ii,k)=missing_val
-!     endif
-!     enddo
-!    enddo
-
-!    do ii=1,Ncolumns
-!     if(isccp%boxtau(ip,ii).eq.R_UNDEF)then
-!       isccp%boxtau(ip,ii)=missing_val
-!     endif
-!    enddo
-
-!    do ii=1,Ncolumns
-!     if(isccp%boxptop(ip,ii).eq.R_UNDEF)then
-!      isccp%boxptop(ip,ii)=missing_val
-!     endif
-!    enddo
-!   enddo
   where(isccp%totalcldarea == R_UNDEF) isccp%totalcldarea = missing_val
   where(isccp%meanptop == R_UNDEF) isccp%meanptop = missing_val
@@ -349,13 +256,4 @@
 !!! MISR simulator
  if (cfg%Lmisr_sim) then
-!   do ip=1,Npoints
-!     do ii=1,7
-!       do k=1,MISR_N_CTH
-!        if(misr%fq_MISR(ip,ii,k).eq.R_UNDEF)then
-!            misr%fq_MISR(ip,ii,k)=missing_val
-!        endif
-!       enddo
-!     enddo
-!   enddo
    where(misr%fq_MISR == R_UNDEF) misr%fq_MISR = missing_val
 
@@ -400,49 +298,4 @@
         modis%Ice_Water_Path_Mean = missing_val
 
-!  de ip=1,Npoints
-!    if(modis%Cloud_Fraction_Low_Mean(ip).eq.R_UNDEF)then
-!      modis%Cloud_Fraction_Low_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Fraction_High_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Fraction_High_Mean(ip)=missing_val 
-!    endif
-!    if(modis%Cloud_Fraction_Mid_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Fraction_Mid_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Fraction_Total_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Fraction_Total_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Fraction_Water_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Fraction_Water_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Fraction_Ice_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Fraction_Ice_Mean(ip)=missing_val
-!    endif
-!    if(modis%Optical_Thickness_Total_Mean(ip).eq.R_UNDEF)then
-!       modis%Optical_Thickness_Total_Mean(ip)=missing_val
-!    endif
-!    if(modis%Optical_Thickness_Water_Mean(ip).eq.R_UNDEF)then
-!       modis%Optical_Thickness_Water_Mean(ip)=missing_val
-!    endif
-!    if(modis%Optical_Thickness_Ice_Mean(ip).eq.R_UNDEF)then
-!       modis%Optical_Thickness_Ice_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Particle_Size_Water_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Particle_Size_Water_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Particle_Size_Ice_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Particle_Size_Ice_Mean(ip)=missing_val
-!    endif
-!    if(modis%Cloud_Top_Pressure_Total_Mean(ip).eq.R_UNDEF)then
-!       modis%Cloud_Top_Pressure_Total_Mean(ip)=missing_val
-!    endif
-!    if(modis%Liquid_Water_Path_Mean(ip).eq.R_UNDEF)then
-!       modis%Liquid_Water_Path_Mean(ip)=missing_val
-!    endif
-!    if(modis%Ice_Water_Path_Mean(ip).eq.R_UNDEF)then
-!       modis%Ice_Water_Path_Mean(ip)=missing_val
-!    endif
-!  enddo
-
     where(modis%Optical_Thickness_Total_LogMean == R_UNDEF) &
           modis%Optical_Thickness_Total_LogMean = missing_val
@@ -472,13 +325,4 @@
    CALL histwrite2d_cosp(o_iwpmodis,modis%Ice_Water_Path_Mean)
 
-!   do ip=1,Npoints
-!     do ii=1,7
-!       do k=1,7
-!       if(modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k).eq.R_UNDEF)then
-!          modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k)=missing_val
-!        endif
-!       enddo
-!     enddo
-!    enddo 
     where(modis%Optical_Thickness_vs_Cloud_Top_Pressure == R_UNDEF) &
           modis%Optical_Thickness_vs_Cloud_Top_Pressure = missing_val
Index: /LMDZ6/trunk/libf/phylmd/cosp/lidar_simulator.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/lidar_simulator.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/lidar_simulator.F90	(revision 3241)
@@ -121,4 +121,5 @@
 !---------------------------------------------------------------------------------
 
+      USE MOD_COSP_CONSTANTS, only : ok_debug_cosp
       IMPLICIT NONE
       REAL :: SRsat
@@ -201,4 +202,6 @@
       REAL pnorm_perp_liq(npoints,nlev) ! perpendicular lidar backscattered signal power for liq
 
+      REAL :: seuil 
+
 ! Output variable
       REAL pnorm_perp_tot (npoints,nlev) ! perpendicular lidar backscattered signal power
@@ -235,4 +238,9 @@
          Gamma = 0.009
 
+  if (ok_debug_cosp) then
+     seuil=1.e-15
+  else
+     seuil=0.0
+  endif 
 !------------------------------------------------------------
 !---- 1. Preliminary definitions and calculations :
@@ -494,6 +502,6 @@
 
       DO k= nlev, 1, -1
-	      pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
-	      pnorm_perp_liq(:,k) = 1000*Beta * pnorm_liq(:,k)**2 + Gamma * pnorm_liq(:,k) ! Liquid particles
+         pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
+         pnorm_perp_liq(:,k) = 1000*Beta * pnorm_liq(:,k)**2 + Gamma * pnorm_liq(:,k) ! Liquid particles
       ENDDO
 
@@ -523,6 +531,6 @@
           tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1) 
         WHERE (tautot_lay_liq(:).GT.0.)
-         beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ EXP(-2.0*tautot_liq(:,k+1)) * (2.*tautot_lay_liq(:)) &
-            & / (1.-exp(-2.0*tautot_lay_liq(:)))
+         beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ max(seuil,EXP(-2.0*tautot_liq(:,k+1))) &
+            & * (2.*tautot_lay_liq(:)) / (1.-exp(-2.0*tautot_lay_liq(:)))
 
         ELSEWHERE
Index: /LMDZ6/trunk/libf/phylmd/cosp/m_mrgrnk.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/m_mrgrnk.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/m_mrgrnk.F90	(revision 3241)
@@ -1,3 +1,4 @@
 Module m_mrgrnk
+IMPLICIT NONE
 Integer, Parameter :: kdp = selected_real_kind(15)
 public :: mrgrnk
Index: /LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_constants.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_constants.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_constants.F90	(revision 3241)
@@ -38,4 +38,7 @@
 
     character(len=32) :: COSP_VERSION='COSP v1.4'
+
+    ! Mode debug ou prod (AI 0302018)
+    logical :: ok_debug_cosp = .false.
 
     ! Indices to address arrays of LS and CONV hydrometeors
Index: /LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_utils.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_utils.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_utils.F90	(revision 3241)
@@ -61,5 +61,12 @@
     integer :: i,j,k
     real :: sigma,one_over_xip1,xi,rho0,rho,lambda_x,gamma_4_3_2,delta
-    
+    real :: seuil 
+
+    if (ok_debug_cosp) then
+       seuil=1.e-15
+    else
+       seuil=0.0
+    endif
+
     mxratio = 0.0
 
@@ -80,5 +87,6 @@
                         mxratio(i,j,k)=mxratio(i,j,k)/rho
                         ! Compute effective radius
-                        if ((reff(i,j,k) <= 0.0).and.(flux(i,k) /= 0.0)) then
+!                        if ((reff(i,j,k) <= 0.0).and.(flux(i,k) /= 0.0)) then
+                        if ((reff(i,j,k) <= 0.0).and.(flux(i,k) > seuil)) then
                            lambda_x = (a_x*c_x*((rho0/rho)**g_x)*n_ax*gamma1/flux(i,k))**(1./delta)
                            reff(i,j,k) = gamma_4_3_2/lambda_x
Index: /LMDZ6/trunk/libf/phylmd/cosp/mod_modis_sim.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/mod_modis_sim.F90	(revision 3240)
+++ /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
Index: /LMDZ6/trunk/libf/phylmd/cosp/radar_simulator_types.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/radar_simulator_types.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/radar_simulator_types.F90	(revision 3241)
@@ -4,5 +4,5 @@
 ! Part of QuickBeam v1.03 by John Haynes
 ! Updated by Roj Marchand June 2010 
-
+  implicit none
   integer, parameter ::       &
   maxhclass = 20         ,& ! max number of hydrometeor classes
Index: /LMDZ6/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90	(revision 3240)
+++ /LMDZ6/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90	(revision 3241)
@@ -10,5 +10,5 @@
     USE xios, ONLY: xios_field_is_active
 #endif
-
+  implicit none
   character(len=*),intent(in) :: cosp_nl
   type(cosp_config),intent(out) :: cfg
