patch.v1.4.0_to_v1.4.2.tar

patch.cosp_modis_simulator.F90.txt0000644000446700022600000001650313211604734017365 0ustar  dswalesclimate--- a/cosp_modis_simulator.F90
+++ b/cosp_modis_simulator.F90
@@ -64,6 +64,8 @@ MODULE MOD_COSP_Modis_Simulator
      ! Also need the ISCCP-type optical thickness/cloud top pressure histogram
      !
      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 
   
 contains
@@ -114,7 +116,11 @@ contains
         meanLiquidWaterPath, meanIceWaterPath
         
     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
     integer, dimension(count(gridBox%sunlit(:) <= 0)) :: notSunlit
@@ -213,18 +219,29 @@ contains
                                 isccpTau(i, :), isccpCloudTopPressure(i, :),                &
                                 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
       !
@@ -253,11 +270,12 @@ contains
       modisSim%Ice_Water_Path_Mean   (sunlit(:)) = meanIceWaterPath
       
       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 
         !
         ! Where it's night and we haven't done the retrievals the values are undefined
@@ -287,6 +305,8 @@ contains
         modisSim%Ice_Water_Path_Mean   (notSunlit(:)) = R_UNDEF
   
         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
       !
@@ -317,6 +337,8 @@ contains
       modisSim%Ice_Water_Path_Mean   (:) = R_UNDEF
   
       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 
 
   end subroutine COSP_Modis_Simulator
@@ -362,7 +384,11 @@ contains
     allocate(x%Ice_Water_Path_Mean(x%nPoints)) 
       
     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
+    x%Optical_Thickness_vs_ReffICE(:, :, :)            = R_UNDEF
+    x%Optical_Thickness_vs_ReffLIQ(:, :, :)            = R_UNDEF
   END SUBROUTINE CONSTRUCT_COSP_MODIS
 
   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -399,6 +425,8 @@ contains
     if(associated(x%Ice_Water_Path_Mean   )) deallocate(x%Ice_Water_Path_Mean   ) 
     
     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
   ! -----------------------------------------------------
 
@@ -446,7 +474,12 @@ contains
     copy%Ice_Water_Path_Mean   (copy_start:copy_end) = orig%Ice_Water_Path_Mean  (orig_start:orig_end)
     
     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
   ! -----------------------------------------------------
 
patch.cosp_simulator.F90.txt0000644000446700022600000000176313211604772016176 0ustar  dswalesclimate--- a/cosp_simulator.F90
+++ b/cosp_simulator.F90
@@ -204,6 +204,14 @@ SUBROUTINE COSP_SIMULATOR(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,m
     where(modis%Optical_Thickness_vs_Cloud_Top_Pressure /= R_UNDEF) modis%Optical_Thickness_vs_Cloud_Top_Pressure = &
                                                       modis%Optical_Thickness_vs_Cloud_Top_Pressure*100.0
   endif
+  if (cfg%Lcrimodis) then
+     where(modis%Optical_Thickness_vs_ReffICE /= R_UNDEF) modis%Optical_Thickness_vs_ReffICE = &
+                                                      modis%Optical_Thickness_vs_ReffICE*100.0
+  endif
+  if (cfg%Lcrlmodis) then
+     where(modis%Optical_Thickness_vs_ReffLIQ /= R_UNDEF) modis%Optical_Thickness_vs_ReffLIQ = &
+                                                      modis%Optical_Thickness_vs_ReffLIQ*100.0
+  endif
   if (cfg%Lcltmodis) then
      where(modis%Cloud_Fraction_Total_Mean /= R_UNDEF) modis%Cloud_Fraction_Total_Mean = modis%Cloud_Fraction_Total_Mean*100.0
   endif
patch.cosp_types.F90.txt0000644000446700022600000000102613211604720015304 0ustar  dswalesclimate--- a/cosp_types.F90
+++ b/cosp_types.F90
@@ -50,7 +50,7 @@ MODULE MOD_COSP_TYPES
                 Lfracout,LlidarBetaMol532,Ltbrttov, &
                 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)
   END TYPE COSP_CONFIG
patch.modis_simulator.F90.txt0000644000446700022600000005342313211605011016330 0ustar  dswalesclimate--- a/MODIS_simulator/modis_simulator.F90
+++ b/MODIS_simulator/modis_simulator.F90
@@ -78,8 +78,9 @@ module mod_modis_sim
   !
   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
   !
@@ -124,6 +125,30 @@ module mod_modis_sim
   real, parameter, dimension(numPressureHistogramBins) ::                         &
     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: 
   !  1) Provide total optical thickness and liquid/ice water content and we'll partition tau in 
@@ -383,6 +408,234 @@ contains
                                     retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
                                 
   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.
+    ! ########################################################################################
+    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
+    Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
+    Cloud_Fraction_Ice_Mean(1:nPoints)   = Cloud_Fraction_Ice_Mean(1:nPoints)   /nSubcols
+    Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_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,            &
        Cloud_Fraction_Total_Mean,       Cloud_Fraction_Water_Mean,       Cloud_Fraction_Ice_Mean,       &
@@ -609,10 +862,10 @@ contains
     g(:) = 0; w0(:) = 0. 
     tau(:) = ice_tau(:) + water_tau(:) 
     where (tau(:) > 0) 
-      g(:)  = (water_tau(:) * water_g(:)                + ice_tau(:) * ice_g(:)            ) / & 
-              tau(:) 
-      w0(:) = (water_tau(:) * water_g(:) * water_w0(:)  + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
-              (g(:) * tau(:))
+      w0(:) = (water_tau(:) * water_w0(:)  + ice_tau(:) * ice_w0(:)) / &
+              tau(:)
+      g(:) = (water_tau(:) * water_g(:) * water_w0(:)  + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
+             (w0(:) * tau(:))
     end where
     
     compute_nir_reflectance = compute_toa_reflectace(tau, g, w0)
@@ -665,17 +918,19 @@ contains
       !
       ! 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
     end if 
@@ -738,22 +993,22 @@ contains
     integer, intent(in) :: phase
     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)
     end if 
@@ -770,11 +1025,10 @@ contains
         !   of size for ice and water
         ! 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)
           if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
@@ -1017,4 +1271,4 @@ contains
     stop
   end subroutine complain_and_die
   !------------------------------------------------------------------------------------------------
-end module mod_modis_sim
\ No newline at end of file
+end module mod_modis_sim
README0000644000446700022600000000211413211605337011604 0ustar  dswalesclimate####################################################################################	 
COSPv1.4.2 patch for v1.4.0.
Purpose: This patch updates the MODIS simulator to create additional diagnostics for
         joint-PDFs of optical-depth and effective radius, for both ice and liquid
	 particles. Additionally, an error in the MODIS simulator has been corrected.
	 This error impacts only mixed-phase clouds visible at the top of the
	 atmosphere.
Dustin Swales (dusitn.swales@noaa.gov)
####################################################################################	 
	 
There are 4 COSP source files impacted by this update:
1) MODIS_simulator/modis_simulator.F90
2) cosp_types.F90
3) cosp_modis_simulator.F90
4) cosp_simulator.F90

The following commands (patch FileToPatch Patch)  will apply the patches:
patch MODIS_simulator/modis_simulator.F90 patch.modis_simulator.F90.txt
patch cosp_types.F90                      patch.cosp_types.F90.txt
patch cosp_modis_simulator.F90            patch.cosp_modis_simulator.F90.txt
patch cosp_simulator.F90                  patch.cosp_simulator.F90.txt