Index: LMDZ5/trunk/libf/phylmd/cosp/cosp_constants.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/cosp_constants.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/cosp_constants.F90	(revision 2713)
@@ -35,6 +35,4 @@
 #include "cosp_defs.h"
 MODULE MOD_COSP_CONSTANTS
-
-    use netcdf, only: nf90_fill_real
     IMPLICIT NONE
 
@@ -54,9 +52,8 @@
     ! Missing value
     real,parameter :: R_UNDEF = -1.0E30
-!    real,parameter :: R_UNDEF = nf90_fill_real
 
     ! Number of possible output variables
-    integer,parameter :: N_OUT_LIST = 63
-    integer,parameter :: N3D = 8
+    integer,parameter :: N_OUT_LIST = 65
+    integer,parameter :: N3D = 10
     integer,parameter :: N2D = 14
     integer,parameter :: N1D = 40
@@ -108,6 +105,5 @@
                    -31.5,-28.5,-25.5,-22.5,-19.5,-16.5,-13.5,-10.5, -7.5, -4.5, &
                     -1.5,  1.5,  4.5,  7.5, 10.5, 13.5, 16.5, 19.5, 22.5, 25.5/)
-    real,parameter,dimension(2,LIDAR_NTEMP) :: LIDAR_PHASE_TEMP_BNDS=reshape(source=&
-                    (/-273.15,-90.,-90.,-87.,-87.,-84.,-84.,-81.,-81.,-78., &
+    real,parameter,dimension(2,LIDAR_NTEMP) :: LIDAR_PHASE_TEMP_BNDS=reshape(source=(/-273.15,-90.,-90.,-87.,-87.,-84.,-84.,-81.,-81.,-78., &
                    -78.,-75.,-75.,-72.,-72.,-69.,-69.,-66.,-66.,-63., &
                    -63.,-60.,-60.,-57.,-57.,-54.,-54.,-51.,-51.,-48., &
Index: LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90	(revision 2713)
@@ -2,5 +2,5 @@
 !   Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences
 ! All rights reserved.
-! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
+! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_modis_simulator.F90 $
 ! 
@@ -65,4 +65,6 @@
      !
      real, dimension(:, :, :), pointer :: Optical_Thickness_vs_Cloud_Top_Pressure
+     real, dimension(:, :, :), pointer :: Optical_Thickness_vs_ReffICE
+     real, dimension(:, :, :), pointer :: Optical_Thickness_vs_ReffLIQ
   end type COSP_MODIS 
   
@@ -115,5 +117,9 @@
         
     real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numModisPressureBins) :: & 
-       jointHistogram
+         jointHistogram
+    real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numMODISReffIceBins) :: & 
+         jointHistogram2
+    real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numMODISReffLiqBins) :: & 
+         jointHistogram3
     
     integer, dimension(count(gridBox%sunlit(:) >  0)) :: sunlit
@@ -214,16 +220,27 @@
                                 retrievedPhase(i, :), retrievedCloudTopPressure(i, :),      & 
                                 retrievedTau(i, :), retrievedSize(i, :))
-      end do
-      call modis_L3_simulator(retrievedPhase,              &
-                              retrievedCloudTopPressure,   &
-                              retrievedTau, retrievedSize, &
-                              cfTotal,         cfLiquid,         cfIce,         &
-                              cfHigh,          cfMid,            cfLow,         &
-                              meanTauTotal,    meanTauLiquid,    meanTauIce,    &
-                              meanLogTauTotal, meanLogTauLiquid, meanLogTauIce, &
-                                               meanSizeLiquid,   meanSizeIce,   &
-                              meanCloudTopPressure,                             &
-                                               meanLiquidWaterPath, meanIceWaterPath, &
-                              jointHistogram)
+     end do
+     
+      ! DJS2015: Call L3 modis simulator used by cospv2.0
+     ! call modis_L3_simulator(retrievedPhase,              &
+     !                         retrievedCloudTopPressure,   &
+     !                         retrievedTau, retrievedSize, &
+     !                         cfTotal,         cfLiquid,         cfIce,         &
+     !                         cfHigh,          cfMid,            cfLow,         &
+     !                         meanTauTotal,    meanTauLiquid,    meanTauIce,    &
+     !                         meanLogTauTotal, meanLogTauLiquid, meanLogTauIce, &
+     !                         meanSizeLiquid,   meanSizeIce,   &
+     !                         meanCloudTopPressure,                             &
+     !                         meanLiquidWaterPath, meanIceWaterPath, &
+     !                         jointHistogram)
+     call modis_column(nSunlit,nSubcols,retrievedPhase,retrievedCloudTopPressure,   &
+                        retrievedTau,retrievedSize,cfTotal,cfLiquid,cfIce,cfHigh,    &
+                        cfMid,cfLow,meanTauTotal,meanTauLiquid,meanTauIce,           &
+                        meanLogTauTotal,meanLogTauLiquid,meanLogTauIce,              &
+                        meanSizeLiquid,meanSizeIce,meanCloudTopPressure,             &
+                        meanLiquidWaterPath, meanIceWaterPath,                       &
+                        jointHistogram,jointHistogram2,jointHistogram3)
+      ! DJS2015: END
+      
       !
       ! Copy results into COSP structure
@@ -254,9 +271,10 @@
       
       modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(sunlit(:), 2:numModisTauBins+1, :) = jointHistogram(:, :, :)
+      modisSim%Optical_Thickness_vs_ReffICE(sunlit(:),2:numModisTauBins+1,:)              = jointHistogram2(:, :, :)
+      modisSim%Optical_Thickness_vs_ReffLIQ(sunlit(:),2:numModisTauBins+1,:)              = jointHistogram3(:, :, :)
       ! 
       ! Reorder pressure bins in joint histogram to go from surface to TOA 
       !
-      modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:,:,:) = &
-        modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, numModisPressureBins:1:-1)
+      modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:,:,:) = modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, numModisPressureBins:1:-1)
       if(nSunlit < nPoints) then 
         !
@@ -288,4 +306,6 @@
   
         modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(notSunlit(:), :, :) = R_UNDEF
+        modisSim%Optical_Thickness_vs_ReffICE(notSunlit(:), :, :) = R_UNDEF
+        modisSim%Optical_Thickness_vs_ReffLIQ(notSunlit(:), :, :) = R_UNDEF
       end if 
     else
@@ -318,4 +338,6 @@
   
       modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, :) = R_UNDEF
+      modisSim%Optical_Thickness_vs_ReffICE(:, :, :) = R_UNDEF
+      modisSim%Optical_Thickness_vs_ReffLIQ(:, :, :) = R_UNDEF
     end if 
 
@@ -363,4 +385,6 @@
       
     allocate(x%Optical_Thickness_vs_Cloud_Top_Pressure(nPoints, numModisTauBins+1, numModisPressureBins))
+    allocate(x%Optical_Thickness_vs_ReffICE(nPoints, numModisTauBins+1, numModisReffIceBins))
+    allocate(x%Optical_Thickness_vs_ReffLIQ(nPoints, numModisTauBins+1, numModisReffLiqBins))
     x%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, :) = R_UNDEF
   END SUBROUTINE CONSTRUCT_COSP_MODIS
@@ -400,4 +424,6 @@
     
     if(associated(x%Optical_Thickness_vs_Cloud_Top_Pressure)) deallocate(x%Optical_Thickness_vs_Cloud_Top_Pressure   ) 
+    if(associated(x%Optical_Thickness_vs_ReffIce)) deallocate(x%Optical_Thickness_vs_ReffIce) 
+    if(associated(x%Optical_Thickness_vs_ReffLiq)) deallocate(x%Optical_Thickness_vs_ReffLiq) 
   END SUBROUTINE FREE_COSP_MODIS
   ! -----------------------------------------------------
@@ -447,5 +473,10 @@
     
     copy%Optical_Thickness_vs_Cloud_Top_Pressure(copy_start:copy_end, :, :) = &
-                          orig%Optical_Thickness_vs_Cloud_Top_Pressure(orig_start:orig_end, :, :)
+         orig%Optical_Thickness_vs_Cloud_Top_Pressure(orig_start:orig_end, :, :)
+    copy%Optical_Thickness_vs_ReffIce(copy_start:copy_end, :, :) = &
+         orig%Optical_Thickness_vs_ReffIce(orig_start:orig_end, :, :)
+    copy%Optical_Thickness_vs_ReffLiq(copy_start:copy_end, :, :) = &
+         orig%Optical_Thickness_vs_ReffLiq(orig_start:orig_end, :, :)
+
   END SUBROUTINE COSP_MODIS_CPSECTION
   ! -----------------------------------------------------
Index: LMDZ5/trunk/libf/phylmd/cosp/cosp_output_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/cosp_output_mod.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/cosp_output_mod.F90	(revision 2713)
@@ -9,4 +9,6 @@
   USE MOD_COSP_TYPES
   use MOD_COSP_Modis_Simulator, only : cosp_modis
+  use mod_modis_sim, only : numMODISReffIceBins, reffICE_binCenters, &
+                            numMODISReffLiqBins, reffLIQ_binCenters
 
 ! cosp_output_mod
@@ -17,7 +19,9 @@
 !$OMP THREADPRIVATE(cosp_outfilekeys, cosp_nidfiles)
       INTEGER, DIMENSION(3), SAVE  :: nhoricosp,nvert,nvertmcosp,nvertcol,nvertbze, &
-                                      nvertsratio,nvertisccp,nvertp,nverttemp,nvertmisr
+                                      nvertsratio,nvertisccp,nvertp,nverttemp,nvertmisr, &
+                                      nvertReffIce,nvertReffLiq
       REAL, DIMENSION(3), SAVE                :: zoutm_cosp
 !$OMP THREADPRIVATE(nhoricosp, nvert,nvertmcosp,nvertcol,nvertsratio,nvertbze,nvertisccp,nvertp,zoutm_cosp,nverttemp,nvertmisr)
+!$OMP THREADPRIVATE(nvertReffIce,nvertReffLiq)
       REAL, SAVE                   :: zdtimemoy_cosp
 !$OMP THREADPRIVATE(zdtimemoy_cosp) 
@@ -176,4 +180,8 @@
   TYPE(ctrl_outcosp), SAVE :: o_clmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
          "clmodis", "MODIS Cloud Area Fraction", "%", (/ ('', i=1, 3) /))
+  TYPE(ctrl_outcosp), SAVE :: o_crimodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
+         "crimodis", "Optical_Thickness_vs_ReffIce from Modis", "%", (/ ('',i=1, 3) /))          
+  TYPE(ctrl_outcosp), SAVE :: o_crlmodis = ctrl_outcosp((/ .TRUE., .TRUE.,.TRUE. /), &
+         "crlmodis", "Optical_Thickness_vs_ReffLiq from Modis", "%", (/ ('',i=1, 3) /))         
 
 ! Rttovs simulator
@@ -325,4 +333,7 @@
 !   CALL wxios_add_vaxis("dbze", DBZE_BINS, dbze_ax) 
 !   CALL wxios_add_vaxis("scatratio", SR_BINS, sratio_ax)
+   CALL wxios_add_vaxis("ReffIce", numMODISReffIceBins, reffICE_binCenters)
+   CALL wxios_add_vaxis("ReffLiq", numMODISReffLiqBins, reffLIQ_binCenters)
+
 #endif
    
@@ -366,4 +377,10 @@
 
       CALL histvert(cosp_nidfiles(iff),"cth16","altitude","m",MISR_N_CTH,MISR_CTH,nvertmisr(iff))
+  
+      CALL histvert(cosp_nidfiles(iff),"ReffIce","Effective_particle_size_Ice","microns",numMODISReffIceBins, reffICE_binCenters, &
+                    nvertReffIce(iff))                                         
+     
+      CALL histvert(cosp_nidfiles(iff),"ReffLiq","Effective_particle_size_Liq","microns",numMODISReffLiqBins, reffLIQ_binCenters, &                                  
+                    nvertReffLiq(iff))
 
 !      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))
Index: LMDZ5/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 2713)
@@ -372,4 +372,17 @@
  endif
 
+    where(modis%Optical_Thickness_vs_ReffIce == R_UNDEF) &
+          modis%Optical_Thickness_vs_ReffIce = Cosp_fill_value
+
+    where(modis%Optical_Thickness_vs_ReffLiq == R_UNDEF) &
+          modis%Optical_Thickness_vs_ReffLiq = Cosp_fill_value
+
+    do icl=1,7
+   CALL histwrite3d_cosp(o_crimodis, &
+     modis%Optical_Thickness_vs_ReffIce(:,icl,:),nvertReffIce,icl)
+   CALL histwrite3d_cosp(o_crlmodis, &
+     modis%Optical_Thickness_vs_ReffLiq(:,icl,:),nvertReffLiq,icl)
+    enddo
+
  IF(.NOT.cosp_varsdefined) THEN
 !$OMP MASTER
@@ -521,5 +534,11 @@
       ELSE IF (nvertsave.eq.nvertmisr(iff)) THEN
           klevs=MISR_N_CTH
-           nam_axvert="cth16"
+          nam_axvert="cth16"
+      ELSE IF (nvertsave.eq.nvertReffIce(iff)) THEN
+          klevs= numMODISReffIceBins
+          nam_axvert="ReffIce"
+      ELSE IF (nvertsave.eq.nvertReffLiq(iff)) THEN
+          klevs= numMODISReffLiqBins
+          nam_axvert="ReffLiq"
       ELSE
            klevs=Nlevout
Index: LMDZ5/trunk/libf/phylmd/cosp/cosp_types.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/cosp_types.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/cosp_types.F90	(revision 2713)
@@ -51,5 +51,5 @@
                 Lcltmodis,Lclwmodis,Lclimodis,Lclhmodis,Lclmmodis,Lcllmodis,Ltautmodis,Ltauwmodis,Ltauimodis,Ltautlogmodis, &
                 Ltauwlogmodis,Ltauilogmodis,Lreffclwmodis,Lreffclimodis,Lpctmodis,Llwpmodis, &
-                Liwpmodis,Lclmodis
+                Liwpmodis,Lclmodis,Lcrimodis,Lcrlmodis
 
      character(len=32) :: out_list(N_OUT_LIST)
Index: LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90	(revision 2713)
@@ -2,5 +2,5 @@
 !   Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences
 ! All rights reserved.
-! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
+! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
 ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $
 ! 
@@ -79,6 +79,7 @@
   real,    parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90.
   integer, parameter :: num_trial_res = 15             ! increase to make the linear pseudo-retrieval of size more accurate
-  logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? 
-  
+! DJS2015: Remove unused parameter
+!  logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? 
+! DJS2015 END  
   !
   ! Precompute near-IR optical params vs size for retrieval scheme
@@ -125,4 +126,28 @@
     nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + &
                                        nominalPressureHistogramBoundaries(2, :) ) / 2.
+  ! DJS2015 START: Add bin descriptions for joint-histograms of partice-sizes and optical depth. This is
+  !                identical to what is done in COSPv.2.0.0 for histogram bin initialization. 
+  integer :: j
+  integer,parameter :: &
+       numMODISReffLiqBins = 6, & ! Number of bins for tau/ReffLiq joint-histogram
+       numMODISReffIceBins = 6    ! Number of bins for tau/ReffICE joint-histogram
+  real,parameter,dimension(numMODISReffLiqBins+1) :: &
+       reffLIQ_binBounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/)
+  real,parameter,dimension(numMODISReffIceBins+1) :: &
+       reffICE_binBounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/)
+  real,parameter,dimension(2,numMODISReffIceBins) :: &
+       reffICE_binEdges = reshape(source=(/reffICE_binBounds(1),((reffICE_binBounds(k),  &
+                                  l=1,2),k=2,numMODISReffIceBins),reffICE_binBounds(numMODISReffIceBins+1)/),  &
+                                  shape = (/2,numMODISReffIceBins/)) 
+  real,parameter,dimension(2,numMODISReffLiqBins) :: &
+       reffLIQ_binEdges = reshape(source=(/reffLIQ_binBounds(1),((reffLIQ_binBounds(k),  &
+                                  l=1,2),k=2,numMODISReffLiqBins),reffLIQ_binBounds(numMODISReffIceBins+1)/),  &
+                                  shape = (/2,numMODISReffLiqBins/))             
+  real,parameter,dimension(numMODISReffIceBins) :: &
+       reffICE_binCenters = (reffICE_binEdges(1,:)+reffICE_binEdges(2,:))/2.
+  real,parameter,dimension(numMODISReffLiqBins) :: &
+       reffLIQ_binCenters = (reffLIQ_binEdges(1,:)+reffLIQ_binEdges(2,:))/2.
+  ! DJS2015 END
+
   ! ------------------------------
   ! There are two ways to call the MODIS simulator: 
@@ -384,4 +409,230 @@
                                 
   end subroutine modis_L2_simulator_oneTau
+
+  ! ########################################################################################
+  subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size,      &
+       Cloud_Fraction_Total_Mean,         Cloud_Fraction_Water_Mean,         Cloud_Fraction_Ice_Mean,        &
+       Cloud_Fraction_High_Mean,          Cloud_Fraction_Mid_Mean,           Cloud_Fraction_Low_Mean,        &
+       Optical_Thickness_Total_Mean,      Optical_Thickness_Water_Mean,      Optical_Thickness_Ice_Mean,     &
+       Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
+       Cloud_Particle_Size_Water_Mean,    Cloud_Particle_Size_Ice_Mean,      Cloud_Top_Pressure_Total_Mean,  &
+       Liquid_Water_Path_Mean,            Ice_Water_Path_Mean,                                               &    
+       Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
+    
+    ! INPUTS
+    integer,intent(in) :: &
+         nPoints,                           & ! Number of horizontal gridpoints
+         nSubCols                             ! Number of subcolumns
+    integer,intent(in), dimension(:,:) ::  &
+!ds    integer,intent(in), dimension(nPoints, nSubCols) ::  &
+         phase                             
+    real,intent(in),dimension(:,:) ::  &
+!ds    real,intent(in),dimension(nPoints, nSubCols) ::  &
+         cloud_top_pressure,                &
+         optical_thickness,                 &
+         particle_size
+ 
+    ! OUTPUTS 
+    real,intent(inout),dimension(:)  ::   & !
+!ds    real,intent(inout),dimension(nPoints)  ::   & !
+         Cloud_Fraction_Total_Mean,         & !
+         Cloud_Fraction_Water_Mean,         & !
+         Cloud_Fraction_Ice_Mean,           & !
+         Cloud_Fraction_High_Mean,          & !
+         Cloud_Fraction_Mid_Mean,           & !
+         Cloud_Fraction_Low_Mean,           & !
+         Optical_Thickness_Total_Mean,      & !
+         Optical_Thickness_Water_Mean,      & !
+         Optical_Thickness_Ice_Mean,        & !
+         Optical_Thickness_Total_MeanLog10, & !
+         Optical_Thickness_Water_MeanLog10, & !
+         Optical_Thickness_Ice_MeanLog10,   & !
+         Cloud_Particle_Size_Water_Mean,    & !
+         Cloud_Particle_Size_Ice_Mean,      & !
+         Cloud_Top_Pressure_Total_Mean,     & !
+         Liquid_Water_Path_Mean,            & !
+         Ice_Water_Path_Mean                  !
+    real,intent(inout),dimension(:,:,:) :: &
+!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: &
+         Optical_Thickness_vs_Cloud_Top_Pressure
+    real,intent(inout),dimension(:,:,:) :: &    
+!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: &    
+         Optical_Thickness_vs_ReffIce
+    real,intent(inout),dimension(:,:,:) :: &    
+!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: &    
+         Optical_Thickness_vs_ReffLiq         
+
+    ! LOCAL VARIABLES
+    real, parameter :: &
+         LWP_conversion = 2./3. * 1000. ! MKS units  
+    integer :: i, j
+    logical, dimension(nPoints,nSubCols) :: &
+         cloudMask,      &
+         waterCloudMask, &
+         iceCloudMask,   &
+         validRetrievalMask
+    real,dimension(nPoints,nSubCols) :: &
+         tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
+    
+    ! ########################################################################################
+    ! Include only those pixels with successful retrievals in the statistics 
+    ! ########################################################################################
+    validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
+    cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and.       &
+         validRetrievalMask(1:nPoints,1:nSubCols)
+    waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. &
+         validRetrievalMask(1:nPoints,1:nSubCols)
+    iceCloudMask(1:nPoints,1:nSubCols)   = phase(1:nPoints,1:nSubCols) == phaseIsIce .and.    &
+         validRetrievalMask(1:nPoints,1:nSubCols)
+
+    ! ########################################################################################
+    ! Use these as pixel counts at first 
+    ! ########################################################################################
+    Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask,      dim = 2))
+    Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
+    Cloud_Fraction_Ice_Mean(1:nPoints)   = real(count(iceCloudMask,   dim = 2))
+    Cloud_Fraction_High_Mean(1:nPoints)  = real(count(cloudMask .and. cloud_top_pressure <=          &
+                                           highCloudPressureLimit, dim = 2)) 
+    Cloud_Fraction_Low_Mean(1:nPoints)   = real(count(cloudMask .and. cloud_top_pressure >           &
+                                           lowCloudPressureLimit,  dim = 2)) 
+    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
+                                           - Cloud_Fraction_Low_Mean(1:nPoints)
+
+    ! ########################################################################################
+    ! 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)
+       
+    ! ########################################################################################
+    ! We take the absolute value of optical thickness here to satisfy compilers that complains 
+    ! when we evaluate the logarithm of a negative number, even though it's not included in 
+    ! 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)
+
+    ! ########################################################################################
+    ! Normalize pixel counts to fraction. The first three cloud fractions have been set to -1 
+    ! in cloud-free areas, so set those places to 0.
+    ! ########################################################################################
+    Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
+    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)  /nSubcols
+    Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)  /nSubcols
+
+    ! ########################################################################################
+    ! 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
+    where (Cloud_Fraction_Low_Mean == 0)   Cloud_Fraction_Low_Mean = R_UNDEF
+
+    ! ########################################################################################
+    ! Joint histogram  
+    ! ########################################################################################
+
+    ! Loop over all points
+    tauWRK(1:nPoints,1:nSubCols)     = optical_thickness(1:nPoints,1:nSubCols)
+    ctpWRK(1:nPoints,1:nSubCols)     = cloud_top_pressure(1:nPoints,1:nSubCols)
+    reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
+    reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
+    do j=1,nPoints
+
+       ! Fill clear and optically thin subcolumns with fill
+       where(.not. cloudMask(j,1:nSubCols)) 
+          tauWRK(j,1:nSubCols) = -999.
+          ctpWRK(j,1:nSubCols) = -999.
+       endwhere
+       ! Joint histogram of tau/CTP
+       call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
+                   tauHistogramBoundaries,numTauHistogramBins,&
+                   pressureHistogramBoundaries,numPressureHistogramBins,&
+                   Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins))
+       ! Joint histogram of tau/ReffICE
+       call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols,               &
+                   tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds,         &
+                   numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins))
+       ! Joint histogram of tau/ReffLIQ
+       call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols,               &
+                   tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds,         &
+                   numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins))                   
+
+    enddo   
+    Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = &
+         Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols
+    Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = &
+         Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols
+    Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = &
+         Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols 
+
+  end subroutine modis_column
+  ! ######################################################################################
+  ! SUBROUTINE hist2D
+  ! ######################################################################################
+  subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
+    implicit none
+    
+    ! INPUTS
+    integer, intent(in) :: &
+         npts,  & ! Number of data points to be sorted
+         nbin1, & ! Number of bins in histogram direction 1 
+         nbin2    ! Number of bins in histogram direction 2
+    real,intent(in),dimension(npts) :: &
+         var1,  & ! Variable 1 to be sorted into bins
+         var2     ! variable 2 to be sorted into bins
+    real,intent(in),dimension(nbin1+1) :: &
+         bin1     ! Histogram bin 1 boundaries
+    real,intent(in),dimension(nbin2+1) :: &
+         bin2     ! Histogram bin 2 boundaries
+    ! OUTPUTS
+    real,intent(out),dimension(nbin1,nbin2) :: &
+         jointHist
+    
+    ! LOCAL VARIABLES
+    integer :: ij,ik
+    
+    do ij=2,nbin1+1
+       do ik=2,nbin2+1
+          jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
+               var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))        
+       enddo
+    enddo
+  end subroutine hist2D
+  
   !------------------------------------------------------------------------------------------------
   subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size,            &
@@ -666,15 +917,17 @@
       ! If first retrieval works, can try 2nd iteration using greater re resolution 
       !
-      if(use_two_re_iterations .and. retrieve_re > 0.) then
-        re_min = retrieve_re - delta_re
-        re_max = retrieve_re + delta_re
-        delta_re = (re_max - re_min)/real(num_trial_res-1)
-  
-        trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) 
-        g(:)  = get_g_nir(  phase, trial_re(:))
-        w0(:) = get_ssa_nir(phase, trial_re(:))
-        predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
-        retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) 
-      end if
+! DJS2015: Remove unused piece of code      
+!      if(use_two_re_iterations .and. retrieve_re > 0.) then
+!        re_min = retrieve_re - delta_re
+!        re_max = retrieve_re + delta_re
+!        delta_re = (re_max - re_min)/real(num_trial_res-1)
+!  
+!        trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) 
+!        g(:)  = get_g_nir(  phase, trial_re(:))
+!        w0(:) = get_ssa_nir(phase, trial_re(:))
+!        predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
+!        retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) 
+!      end if
+! DJS2015 END
     else 
       retrieve_re = re_fill
@@ -739,20 +992,20 @@
     real,    intent(in) :: re
     real :: get_g_nir 
-    
-    real, dimension(3), parameter :: ice_coefficients   = (/ 0.7432,  4.5563e-3, -2.8697e-5 /), & 
-                               small_water_coefficients = (/ 0.8027, -1.0496e-2,  1.7071e-3 /), & 
-                                 big_water_coefficients = (/ 0.7931,  5.3087e-3, -7.4995e-5 /) 
-    
-    ! approx. fits from MODIS Collection 5 LUT scattering calculations
-    if(phase == phaseIsLiquid) then
-      if(re < 8.) then 
-        get_g_nir = fit_to_quadratic(re, small_water_coefficients)
-        if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
-      else
-        get_g_nir = fit_to_quadratic(re,   big_water_coefficients)
-        if(re > re_water_max) get_g_nir = fit_to_quadratic(re_water_max, big_water_coefficients)
-      end if 
+
+    real, dimension(3), parameter :: ice_coefficients         = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
+                                     small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
+    real, dimension(4), parameter :: big_water_coefficients   = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
+
+    ! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
+    if(phase == phaseIsLiquid) then 
+       if(re < 7.) then
+          get_g_nir = fit_to_quadratic(re, small_water_coefficients)
+          if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
+       else
+          get_g_nir = fit_to_cubic(re, big_water_coefficients)
+          if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
+       end if
     else
-      get_g_nir = fit_to_quadratic(re, ice_coefficients)
+       get_g_nir = fit_to_quadratic(re, ice_coefficients)
       if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
       if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
@@ -771,9 +1024,8 @@
         ! Fits from Steve Platnick
         !
+        real, dimension(4), parameter :: ice_coefficients   = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
+        real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
         
-        real, dimension(4), parameter :: ice_coefficients   = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /)
-        real, dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /) 
-        
-        ! approx. fits from MODIS Collection 5 LUT scattering calculations
+        ! approx. fits from MODIS Collection 6 LUT scattering calculations
         if(phase == phaseIsLiquid) then
           get_ssa_nir = fit_to_quadratic(re, water_coefficients)
Index: LMDZ5/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90	(revision 2652)
+++ LMDZ5/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90	(revision 2713)
@@ -41,5 +41,5 @@
              Lcltmodis,Lclwmodis,Lclimodis,Lclhmodis,Lclmmodis,Lcllmodis,Ltautmodis,Ltauwmodis,Ltauimodis,Ltautlogmodis, &
              Ltauwlogmodis,Ltauilogmodis,Lreffclwmodis,Lreffclimodis,Lpctmodis,Llwpmodis, &
-             Liwpmodis,Lclmodis
+             Liwpmodis,Lclmodis,Lcrimodis,Lcrlmodis
 
   namelist/COSP_OUTPUT/Lradar_sim,Llidar_sim,Lisccp_sim,Lmodis_sim,Lmisr_sim,Lrttov_sim, &
@@ -57,5 +57,5 @@
              Lcltmodis,Lclwmodis,Lclimodis,Lclhmodis,Lclmmodis,Lcllmodis,Ltautmodis,Ltauwmodis,Ltauimodis,Ltautlogmodis, &
              Ltauwlogmodis,Ltauilogmodis,Lreffclwmodis,Lreffclimodis,Lpctmodis,Llwpmodis, &
-             Liwpmodis,Lclmodis
+             Liwpmodis,Lclmodis,Lcrimodis,Lcrlmodis
    
   do i=1,N_OUT_LIST
@@ -137,4 +137,7 @@
   CALL bcast(Lclmodis)
   CALL bcast(Ltbrttov)
+  CALL bcast(Lcrimodis)
+  CALL bcast(Lcrlmodis)
+
 !$OMP BARRIER
 
@@ -223,4 +226,6 @@
     Liwpmodis=.false.
     Lclmodis=.false.
+    Lcrimodis=.false.
+    Lcrlmodis=.false.
   endif
   if (Lmodis_sim) Lisccp_sim = .true.
@@ -381,4 +386,8 @@
   i = i+1
   if (Lclmodis)         cfg%out_list(i) = 'clmodis'
+  i = i+1
+  if (Lcrimodis)         cfg%out_list(i) = 'crimodis'
+  i = i+1
+  if (Lcrlmodis)         cfg%out_list(i) = 'crlmodis'
     
   if (i /= N_OUT_LIST) then
@@ -459,5 +468,7 @@
   cfg%Liwpmodis=Liwpmodis
   cfg%Lclmodis=Lclmodis
-  
+  cfg%Lcrimodis=Lcrimodis
+  cfg%Lcrlmodis=Lcrlmodis
+
  END SUBROUTINE READ_COSP_OUTPUT_NL
 
