! (c) 2009-2010, Regents of the Unversity of Colorado
!   Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences
! All rights reserved.
! $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 $

! Redistribution and use in source and binary forms, with or without modification, are permitted 
! provided that the following conditions are met:

!     * Redistributions of source code must retain the above copyright notice, this list 
!       of conditions and the following disclaimer.
!     * Redistributions in binary form must reproduce the above copyright notice, this list
!       of conditions and the following disclaimer in the documentation and/or other materials 
!       provided with the distribution.
!     * Neither the name of the Met Office nor the names of its contributors may be used 
!       to endorse or promote products derived from this software without specific prior written 
!       permission.

! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! History:
!   May 2009 - Robert Pincus - Initial version
!   June 2009 - Steve Platnick and Robert Pincus - Simple radiative transfer for size retrievals
!   August 2009 - Robert Pincus - Consistency and bug fixes suggested by Rick Hemler (GFDL) 
!   November 2009 - Robert Pincus - Bux fixes and speed-ups after experience with Rick Hemler using AM2 (GFDL) 
!   January 2010 - Robert Pincus - Added high, middle, low cloud fractions 

! Notes on using the MODIS simulator: 
!  *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 microns, or 
!     optical thickness at 0.67 microns and ice- and liquid-water contents (in consistent units of 
!     your choosing)
!  *) Required input also includes the optical thickness and cloud top pressure 
!     derived from the ISCCP simulator run with parameter top_height = 1. 
!  *) Cloud particle sizes are specified as radii, measured in meters, though within the module we 
!     use units of microns. Where particle sizes are outside the bounds used in the MODIS retrieval
!     libraries (parameters re_water_min, re_ice_min, etc.) the simulator returns missing values (re_fill)

! When error conditions are encountered this code calls the function complain_and_die, supplied at the 
!   bottom of this module. Users probably want to replace this with something more graceful. 

module mod_modis_sim
  USE MOD_COSP_TYPES, only: R_UNDEF, ok_debug_cosp
  implicit none
  ! ------------------------------
  ! Algorithmic parameters

  real, parameter :: ice_density          = 0.93               ! liquid density is 1.  

  ! Retrieval parameters

  real, parameter :: min_OpticalThickness = 0.3,             & ! Minimum detectable optical thickness
                     CO2Slicing_PressureLimit = 700. * 100., & ! Cloud with higher pressures use thermal methods, units Pa
                     CO2Slicing_TauLimit = 1.,               & ! How deep into the cloud does CO2 slicing see? 
                     phase_TauLimit      = 1.,               & ! How deep into the cloud does the phase detection see?
                     size_TauLimit       = 2.,               & ! Depth of the re retreivals
                     phaseDiscrimination_Threshold = 0.7       ! What fraction of total extincton needs to be 
                                                               !  in a single category to make phase discrim. work? 
  real,    parameter :: re_fill= -999.
  integer, parameter :: phaseIsNone = 0, phaseIsLiquid = 1, phaseIsIce = 2, phaseIsUndetermined = 3
  
  logical, parameter :: useSimpleReScheme = .false. 

  ! These are the limits of the libraries for the MODIS collection 5 algorithms 
  !   They are also the limits used in the fits for g and w0

  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
! 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

  integer, private :: i 
  real, dimension(num_trial_res), parameter :: & 
        trial_re_w = re_water_min + (re_water_max - re_water_min)/(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /), &
        trial_re_i = re_ice_min   + (re_ice_max -   re_ice_min)  /(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /)
  
  ! Can't initialze these during compilation, but do in before looping columns in retrievals
  real, dimension(num_trial_res) ::  g_w, g_i, w0_w, w0_i
  ! ------------------------------
  ! Bin boundaries for the joint optical thickness/cloud top pressure histogram

  integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7

  real, private :: dummy_real 
  real, dimension(numTauHistogramBins + 1),      parameter :: &
    tauHistogramBoundaries = (/ min_OpticalThickness, 1.3, 3.6, 9.4, 23., 60., 10000. /) 
  real, dimension(numPressureHistogramBins + 1), parameter :: & ! Units Pa 
    pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /) 
  real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100.

  ! For output - nominal bin centers and  bin boundaries. On output pressure bins are highest to lowest. 

  integer, private :: k, l
  real, parameter, dimension(2, numTauHistogramBins) ::   &
    nominalTauHistogramBoundaries =                       &
        reshape(source = (/ tauHistogramBoundaries(1),    &
                            ((tauHistogramBoundaries(k), l = 1, 2), k = 2, numTauHistogramBins), &
                            100000. /),                    &
                shape = (/2,  numTauHistogramBins /) )
  real, parameter, dimension(numTauHistogramBins) ::                    &
    nominalTauHistogramCenters = (nominalTauHistogramBoundaries(1, :) + &
                                  nominalTauHistogramBoundaries(2, :) ) / 2.
  
  real, parameter, dimension(2, numPressureHistogramBins) :: &
    nominalPressureHistogramBoundaries =                     &
        reshape(source = (/ 100000.,                         &
                            ((pressureHistogramBoundaries(k), l = 1, 2), k = numPressureHistogramBins, 2, -1), &
                            0.  /), &
                shape = (/2,  numPressureHistogramBins /) )
  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 
  !     subroutine modis_L2_simulator_oneTau, or 
  !  2) Provide ice and liquid optical depths in each layer

  interface modis_L2_simulator
    module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus
  end interface 
contains
  !------------------------------------------------------------------------------------------------
  ! MODIS simulator using specified liquid and ice optical thickness in each layer 

  !   Note: this simulator operates on all points; to match MODIS itself night-time 
  !     points should be excluded

  !   Note: the simulator requires as input the optical thickness and cloud top pressure 
  !     derived from the ISCCP simulator run with parameter top_height = 1. 
  !     If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing 
  !     and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that 
  !     alogrithm in this simulator we simply report the values from the ISCCP simulator. 

  subroutine modis_L2_simulator_twoTaus(                                       &
                                temp, pressureLayers, pressureLevels,          &
                                liquid_opticalThickness, ice_opticalThickness, &
                                waterSize, iceSize,                            & 
                                isccpTau, isccpCloudTopPressure,               &
                                retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)

    ! Grid-mean quantities at layer centers, starting at the model top
    !   dimension nLayers
    real, dimension(:),    intent(in ) :: temp,           & ! Temperature, K
                                          pressureLayers, & ! Pressure, Pa
                                          pressureLevels    ! Pressure at layer edges, Pa (dimension nLayers + 1) 
    ! Sub-column quantities
    !   dimension  nSubcols, nLayers
    real, dimension(:, :), intent(in ) :: liquid_opticalThickness, & ! Layer optical thickness @ 0.67 microns due to liquid
                                          ice_opticalThickness       ! ditto, due to ice
    real, dimension(:, :), intent(in ) :: waterSize,        & ! Cloud drop effective radius, microns
                                          iceSize             ! Cloud ice effective radius, microns
                                          
    ! Cloud properties retrieved from ISCCP using top_height = 1
    !    dimension nSubcols
    real, dimension(:),    intent(in ) :: isccpTau, &           ! Column-integrated optical thickness 
                                          isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa) 

    ! Properties retrieved by MODIS
    !   dimension nSubcols
    integer, dimension(:), intent(out) :: retrievedPhase               ! liquid/ice/other - integer, defined in module header
    real,    dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
                                          retrievedTau,              & ! unitless
                                          retrievedSize                ! microns 
    ! ---------------------------------------------------
    ! Local variables
    logical, dimension(size(retrievedTau))                     :: cloudMask
    real,    dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal
    real    :: integratedLiquidFraction
    integer :: i, nSubcols, nLevels

    ! ---------------------------------------------------
    nSubcols = size(liquid_opticalThickness, 1)
    nLevels  = size(liquid_opticalThickness, 2) 

    ! Initial error checks 

    if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), &
              size(isccpTau), size(isccpCloudTopPressure),              &
              size(retrievedPhase), size(retrievedCloudTopPressure),    &
              size(retrievedTau), size(retrievedSize) /) /= nSubcols )) &
       call complain_and_die("Differing number of subcolumns in one or more arrays") 
    
    if(any((/ size(ice_opticalThickness, 2), size(waterSize, 2), size(iceSize, 2),      &
              size(temp), size(pressureLayers), size(pressureLevels)-1 /) /= nLevels )) &
       call complain_and_die("Differing number of levels in one or more arrays") 
       
       
    if(any( (/ any(temp <= 0.), any(pressureLayers <= 0.),  &
               any(liquid_opticalThickness < 0.),           &
               any(ice_opticalThickness < 0.),              &
               any(waterSize < 0.), any(iceSize < 0.) /) )) &
       call complain_and_die("Input values out of bounds") 
             
    ! ---------------------------------------------------

    ! Compute the total optical thickness and the proportion due to liquid in each cell

    where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.) 
      tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :))
    elsewhere
      tauLiquidFraction(:, :) = 0. 
    end  where 
    tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) 

    ! Optical depth retrieval 
    !   This is simply a sum over the optical thickness in each layer 
    !   It should agree with the ISCCP values after min values have been excluded 

    retrievedTau(:) = sum(tauTotal(:, :), dim = 2)

    ! Cloud detection - does optical thickness exceed detection threshold? 

    cloudMask = retrievedTau(:) >= min_OpticalThickness

    ! Initialize initial estimates for size retrievals

    if(any(cloudMask) .and. .not. useSimpleReScheme) then 
      g_w(:)  = get_g_nir(  phaseIsLiquid, trial_re_w(:))
      w0_w(:) = get_ssa_nir(phaseIsLiquid, trial_re_w(:))
      g_i(:)  = get_g_nir(  phaseIsIce,    trial_re_i(:))
      w0_i(:) = get_ssa_nir(phaseIsIce,    trial_re_i(:))
    end if 
    
    do i = 1, nSubCols
      if(cloudMask(i)) then 

        ! Cloud top pressure determination 
        !   MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds
        !   lower than that. 
        !  For CO2 slicing we report the optical-depth weighted pressure, integrating to a specified 
        !    optical depth
        ! This assumes linear variation in p between levels. Linear in ln(p) is probably better, 
        !   though we'd need to deal with the lowest pressure gracefully. 

        retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), &
                                                          pressureLevels,           &
                                                          CO2Slicing_TauLimit)  

        ! Phase determination - determine fraction of total tau that's liquid 
        ! When ice and water contribute about equally to the extinction we can't tell 
        !   what the phase is 

        integratedLiquidFraction = weight_by_extinction(tauTotal(i, :),          &
                                                        tauLiquidFraction(i, :), &
                                                        phase_TauLimit)
        if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then 
          retrievedPhase(i) = phaseIsLiquid
        else if (integratedLiquidFraction <= 1.- phaseDiscrimination_Threshold) then 
          retrievedPhase(i) = phaseIsIce
        else 
          retrievedPhase(i) = phaseIsUndetermined
        end if 

        ! Size determination 

        if(useSimpleReScheme) then 
          !   This is the extinction-weighted size considering only the phase we've chosen 

          if(retrievedPhase(i) == phaseIsIce) then 
            retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :),  &
                                                    iceSize(i, :), &
                                                    (1. - integratedLiquidFraction) * size_TauLimit)
  
          else if(retrievedPhase(i) == phaseIsLiquid) then 
            retrievedSize(i) = weight_by_extinction(liquid_opticalThickness(i, :), &
                                                    waterSize(i, :), &
                                                    integratedLiquidFraction * size_TauLimit)
  
          else
            retrievedSize(i) = 0. 
          end if 
        else
          retrievedSize(i) = 1.0e-06*retrieve_re(retrievedPhase(i), retrievedTau(i), &
                         obs_Refl_nir = compute_nir_reflectance(liquid_opticalThickness(i, :), waterSize(i, :)*1.0e6, & 
                         ice_opticalThickness(i, :),      iceSize(i, :)*1.0e6))
        end if 
      else 

        ! Values when we don't think there's a cloud. 

        retrievedCloudTopPressure(i) = R_UNDEF 
        retrievedPhase(i) = phaseIsNone
        retrievedSize(i) = R_UNDEF 
        retrievedTau(i) = R_UNDEF 
      end if
    end do
    where((retrievedSize(:) < 0.).and.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill

    ! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS 
    !   mimics what MODIS does to first order. 
    !   Of course, ISCCP cloud top pressures are in mb. 

    where(cloudMask(:) .and. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) &
      retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100. 
    
  end subroutine modis_L2_simulator_twoTaus
  !------------------------------------------------------------------------------------------------

  ! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents; 
  !   we'll partition this into ice and liquid optical thickness and call the full MODIS simulator 

  subroutine modis_L2_simulator_oneTau(                                         &
                                temp, pressureLayers, pressureLevels,           &
                                opticalThickness, cloudWater, cloudIce,         &
                                waterSize, iceSize,                             & 
                                isccpTau, isccpCloudTopPressure,                &
                                retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
    ! Grid-mean quantities at layer centers, 
    !   dimension nLayers
    real, dimension(:),    intent(in ) :: temp,           & ! Temperature, K
                                          pressureLayers, & ! Pressure, Pa
                                          pressureLevels    ! Pressure at layer edges, Pa (dimension nLayers + 1) 
    ! Sub-column quantities
    !   dimension nLayers, nSubcols
    real, dimension(:, :), intent(in ) :: opticalThickness, & ! Layer optical thickness @ 0.67 microns
                                          cloudWater,       & ! Cloud water content, arbitrary units
                                          cloudIce            ! Cloud water content, same units as cloudWater
    real, dimension(:, :), intent(in ) :: waterSize,        & ! Cloud drop effective radius, microns
                                          iceSize             ! Cloud ice effective radius, microns

    ! Cloud properties retrieved from ISCCP using top_height = 1
    !    dimension nSubcols
    
    real, dimension(:),    intent(in ) :: isccpTau, &           ! Column-integrated optical thickness 
                                          isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa) 

    ! Properties retrieved by MODIS
    !   dimension nSubcols
    integer, dimension(:), intent(out) :: retrievedPhase               ! liquid/ice/other - integer
    real,    dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
                                          retrievedTau,              & ! unitless
                                          retrievedSize                ! microns (or whatever units 
                                                                       !   waterSize and iceSize are supplied in)
    ! ---------------------------------------------------
    ! Local variables
    real, dimension(size(opticalThickness, 1), size(opticalThickness, 2)) :: & 
           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. 
    elsewhere
      where (cloudWater(:, :) <= 0.) 
        tauLiquidFraction(:, :) = 0. 
      elsewhere 

        ! Geometic optics limit - tau as LWP/re  (proportional to LWC/re) 

! Modif AI 02 2018 
        tauLiquidFraction(:, :) = (cloudWater(:, :)/max(waterSize(:, :), seuil) ) / &
                                  (cloudWater(:, :)/max(waterSize(:, :), seuil) + & 
                                   cloudIce(:, :)/(ice_density * max(iceSize(:, :), seuil)) ) 
      end where
    end where
    liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :) 
    ice_opticalThickness   (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :)
    
    call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels,          &
                                    liquid_opticalThickness, ice_opticalThickness, &
                                    waterSize, iceSize,                            & 
                                    isccpTau, isccpCloudTopPressure,               &
                                    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.
    ! ########################################################################################
    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 
    ! when we evaluate the logarithm of a negative number, even though it's not included in 
    ! the sum. 
    ! ########################################################################################
    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.
    ! ########################################################################################
    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_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,       &
       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)

    ! Inputs; dimension nPoints, nSubcols

    integer, dimension(:, :),   intent(in)  :: phase
    real,    dimension(:, :),   intent(in)  :: cloud_top_pressure, optical_thickness, particle_size

    ! Outputs; dimension nPoints

    real,    dimension(:),      intent(out) :: &
       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
    ! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins 
    real,    dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure
    ! ---------------------------
    ! Local variables

    real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units  
    integer :: i, j
    integer :: nPoints, nSubcols 
    logical, dimension(size(phase, 1), size(phase, 2)) :: &
      cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask
    logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins     ) :: tauMask
    logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask
    ! ---------------------------
    
    nPoints  = size(phase, 1) 
    nSubcols = size(phase, 2) 

    ! Array conformance checks

    if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1),                                &
               size(Cloud_Fraction_Total_Mean),       size(Cloud_Fraction_Water_Mean),       size(Cloud_Fraction_Ice_Mean),    &
               size(Cloud_Fraction_High_Mean),        size(Cloud_Fraction_Mid_Mean),         size(Cloud_Fraction_Low_Mean),    &
               size(Optical_Thickness_Total_Mean),    size(Optical_Thickness_Water_Mean),    size(Optical_Thickness_Ice_Mean), &
               size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), &
               size(Optical_Thickness_Ice_MeanLog10),   size(Cloud_Particle_Size_Water_Mean),    &
               size(Cloud_Particle_Size_Ice_Mean),      size(Cloud_Top_Pressure_Total_Mean),     &
               size(Liquid_Water_Path_Mean),          size(Ice_Water_Path_Mean) /) /= nPoints))  &
      call complain_and_die("Some L3 arrays have wrong number of grid points") 
    if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /)  /= nSubcols)) &
      call complain_and_die("Some L3 arrays have wrong number of subcolumns") 

    ! Include only those pixels with successful retrievals in the statistics 

    validRetrievalMask(:, :) = particle_size(:, :) > 0.
    cloudMask      = phase(:, :) /= phaseIsNone   .and. validRetrievalMask(:, :)
    waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :)
    iceCloudMask   = phase(:, :) == phaseIsIce    .and. validRetrievalMask(:, :)

    ! Use these as pixel counts at first 

    Cloud_Fraction_Total_Mean(:) = real(count(cloudMask,      dim = 2))
    Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2))
    Cloud_Fraction_Ice_Mean(:)   = real(count(iceCloudMask,   dim = 2))
    
    Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2)) 
    Cloud_Fraction_Low_Mean(:)  = real(count(cloudMask .and. cloud_top_pressure >  lowCloudPressureLimit,  dim = 2)) 
    Cloud_Fraction_Mid_Mean(:)  = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:)

    ! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0. 

    where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1. 
    where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1.
    where (Cloud_Fraction_Ice_Mean   == 0) Cloud_Fraction_Ice_Mean   = -1.
    
    Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask,      dim = 2) / Cloud_Fraction_Total_Mean(:) 
    Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
    Optical_Thickness_Ice_Mean   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / Cloud_Fraction_Ice_Mean(:)
   
    ! 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 = sum(log10(abs(optical_thickness)), mask = cloudMask,      dim = 2) / &
                                        Cloud_Fraction_Total_Mean(:)
    Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / &
                                        Cloud_Fraction_Water_Mean(:)
    Optical_Thickness_Ice_MeanLog10   = sum(log10(abs(optical_thickness)), mask = iceCloudMask,   dim = 2) / &
                                        Cloud_Fraction_Ice_Mean(:)
   
    Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
    Cloud_Particle_Size_Ice_Mean   = sum(particle_size, mask = iceCloudMask,   dim = 2) / Cloud_Fraction_Ice_Mean(:)
    
    Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2))
    
    Liquid_Water_Path_Mean = LWP_conversion &
                             * sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) &
                             / Cloud_Fraction_Water_Mean(:)
    Ice_Water_Path_Mean    = LWP_conversion * ice_density &
                             * sum(particle_size * optical_thickness, mask = iceCloudMask,   dim = 2) &
                             / Cloud_Fraction_Ice_Mean(:)

    ! 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_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols)
    Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols)
    Cloud_Fraction_Ice_Mean(:)   = max(0., Cloud_Fraction_Ice_Mean(:)  /nSubcols)
    
    Cloud_Fraction_High_Mean(:)  = Cloud_Fraction_High_Mean(:) /nSubcols
    Cloud_Fraction_Mid_Mean(:)   = Cloud_Fraction_Mid_Mean(:)  /nSubcols
    Cloud_Fraction_Low_Mean(:)   = Cloud_Fraction_Low_Mean(:)  /nSubcols
    
    ! ----
    ! Joint histogram 

    do i = 1, numTauHistogramBins 
      where(cloudMask(:, :)) 
        tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. &
                           optical_thickness(:, :) <  tauHistogramBoundaries(i+1)
      elsewhere
        tauMask(:, :, i) = .false.
      end where
    end do 

    do i = 1, numPressureHistogramBins 
      where(cloudMask(:, :)) 
        pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. &
                                cloud_top_pressure(:, :) <  pressureHistogramBoundaries(i+1)
      elsewhere
        pressureMask(:, :, i) = .false.
      end where
    end do 
    
    do i = 1, numPressureHistogramBins
      do j = 1, numTauHistogramBins
        Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = & 
          real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
      end do 
    end do 
    
  end subroutine modis_L3_simulator
  !------------------------------------------------------------------------------------------------
  function cloud_top_pressure(tauIncrement, pressure, tauLimit) 
    real, dimension(:), intent(in) :: tauIncrement, pressure
    real,               intent(in) :: tauLimit
    real                           :: cloud_top_pressure

    ! Find the extinction-weighted pressure. Assume that pressure varies linearly between 
    !   layers and use the trapezoidal rule.

    real :: deltaX, totalTau, totalProduct
    integer :: i 
    
    totalTau = 0.; totalProduct = 0. 
    do i = 2, size(tauIncrement)
      if(totalTau + tauIncrement(i) > tauLimit) then 
        deltaX = tauLimit - totalTau
        totalTau = totalTau + deltaX

        ! Result for trapezoidal rule when you take less than a full step
        !   tauIncrement is a layer-integrated value

        totalProduct = totalProduct           &
                     + pressure(i-1) * deltaX &
                     + (pressure(i) - pressure(i-1)) * deltaX**2/(2. * tauIncrement(i)) 
      else
        totalTau =     totalTau     + tauIncrement(i) 
        totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2.
      end if 
      if(totalTau >= tauLimit) exit
    end do 
    cloud_top_pressure = totalProduct/totalTau
  end function cloud_top_pressure
  !------------------------------------------------------------------------------------------------
  function weight_by_extinction(tauIncrement, f, tauLimit) 
    real, dimension(:), intent(in) :: tauIncrement, f
    real,               intent(in) :: tauLimit
    real                           :: weight_by_extinction

    ! Find the extinction-weighted value of f(tau), assuming constant f within each layer

    real    :: deltaX, totalTau, totalProduct
    integer :: i 
    
    totalTau = 0.; totalProduct = 0. 
    do i = 1, size(tauIncrement)
      if(totalTau + tauIncrement(i) > tauLimit) then 
        deltaX       = tauLimit - totalTau
        totalTau     = totalTau     + deltaX
        totalProduct = totalProduct + deltaX * f(i) 
      else
        totalTau     = totalTau     + tauIncrement(i) 
        totalProduct = totalProduct + tauIncrement(i) * f(i) 
      end if 
      if(totalTau >= tauLimit) exit
    end do 
    weight_by_extinction = totalProduct/totalTau
  end function weight_by_extinction
  !------------------------------------------------------------------------------------------------
  pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size) 
    real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
    real                           :: compute_nir_reflectance
    
    real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
                                        tau, g, w0
    !----------------------------------------
    water_g(:)  = get_g_nir(  phaseIsLiquid, water_size) 
    water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size) 
    ice_g(:)    = get_g_nir(  phaseIsIce,    ice_size) 
    ice_w0(:)   = get_ssa_nir(phaseIsIce,    ice_size) 

    ! Combine ice and water optical properties

    g(:) = 0; w0(:) = 0. 
    tau(:) = ice_tau(:) + water_tau(:) 
    where (tau(:) > 0) 
      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)
  end function compute_nir_reflectance
  !------------------------------------------------------------------------------------------------
  ! Retreivals
  !------------------------------------------------------------------------------------------------
  elemental function retrieve_re (phase, tau, obs_Refl_nir)
      integer, intent(in) :: phase
      real,    intent(in) :: tau, obs_Refl_nir
      real                :: retrieve_re

      ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in 
      !   MODIS band 7 (near IR)
      ! Uses 
      !  fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables 
      !  two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0
      !  adding-doubling for total reflectance 



      ! Local variables

      real, parameter :: min_distance_to_boundary = 0.01
      real    :: re_min, re_max, delta_re
      integer :: i 
      
      real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir
      ! --------------------------
    
    if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then 
      if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then
        re_min = re_water_min
        re_max = re_water_max
        trial_re(:) = trial_re_w
        g(:)   = g_w(:) 
        w0(:)  = w0_w(:)
      else
        re_min = re_ice_min
        re_max = re_ice_max
        trial_re(:) = trial_re_i
        g(:)   = g_i(:) 
        w0(:)  = w0_i(:)
      end if

      ! 1st attempt at index: w/coarse re resolution

      predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
      retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) 

      ! If first retrieval works, can try 2nd iteration using greater re resolution 

! 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 
    
  end function retrieve_re
  ! --------------------------------------------
  pure function interpolate_to_min(x, y, yobs)
    real, dimension(:), intent(in) :: x, y 
    real,               intent(in) :: yobs
    real                           :: interpolate_to_min

    ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
    !   y must be monotonic in x

    real, dimension(size(x)) :: diff
    integer                  :: nPoints, minDiffLoc, lowerBound, upperBound
    ! ---------------------------------
    nPoints = size(y)
    diff(:) = y(:) - yobs
    minDiffLoc = minloc(abs(diff), dim = 1) 
    
    if(minDiffLoc == 1) then 
      lowerBound = minDiffLoc
      upperBound = minDiffLoc + 1
    else if(minDiffLoc == nPoints) then
      lowerBound = minDiffLoc - 1
      upperBound = minDiffLoc
    else
      if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
        lowerBound = minDiffLoc-1
        upperBound = minDiffLoc
      else 
        lowerBound = minDiffLoc
        upperBound = minDiffLoc + 1
      end if 
    end if 
    
    if(diff(lowerBound) * diff(upperBound) < 0) then     

      ! Interpolate the root position linearly if we bracket the root

      interpolate_to_min = x(upperBound) - & 
                           diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
    else 
      interpolate_to_min = re_fill
    end if 
    

  end function interpolate_to_min
  ! --------------------------------------------
  ! Optical properties
  ! --------------------------------------------
  elemental function get_g_nir (phase, re)

    ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function 
    !   of size for ice and water
    ! Fits from Steve Platnick

    integer, intent(in) :: phase
    real,    intent(in) :: re
    real :: get_g_nir 

    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)
      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 
    
  end function get_g_nir

  ! --------------------------------------------
    elemental function get_ssa_nir (phase, re)
        integer, intent(in) :: phase
        real,    intent(in) :: re
        real                :: get_ssa_nir

        ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function 
        !   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 /)
        
        ! 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)
          if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
        else
          get_ssa_nir = fit_to_cubic(re, ice_coefficients)
          if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
          if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
        end if 

    end function get_ssa_nir
   ! --------------------------------------------
  pure function fit_to_cubic(x, coefficients) 
    real,               intent(in) :: x
    real, dimension(:), intent(in) :: coefficients
    real                           :: fit_to_cubic
    
    
    fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
 end function fit_to_cubic
   ! --------------------------------------------
  pure function fit_to_quadratic(x, coefficients) 
    real,               intent(in) :: x
    real, dimension(:), intent(in) :: coefficients
    real                           :: fit_to_quadratic
    
    
    fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
 end function fit_to_quadratic
  ! --------------------------------------------
  ! Radiative transfer
  ! --------------------------------------------
  pure function compute_toa_reflectace(tau, g, w0)
    real, dimension(:), intent(in) :: tau, g, w0
    real                           :: compute_toa_reflectace
    
    logical, dimension(size(tau))         :: cloudMask
    integer, dimension(count(tau(:) > 0)) :: cloudIndicies
    real,    dimension(count(tau(:) > 0)) :: Refl,     Trans
    real                                  :: Refl_tot, Trans_tot
    integer                               :: i
    ! ---------------------------------------

    ! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation

    cloudMask = tau(:) > 0. 
    cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask) 
    do i = 1, size(cloudIndicies)
      call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
    end do 
                    
    call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot)  
    
    compute_toa_reflectace = Refl_tot
    
  end function compute_toa_reflectace
  ! --------------------------------------------
  pure subroutine two_stream(tauint, gint, w0int, ref, tra) 
    real, intent(in)  :: tauint, gint, w0int
    real, intent(out) :: ref, tra

    ! Compute reflectance in a single layer using the two stream approximation 
    !   The code itself is from Lazaros Oreopoulos via Steve Platnick 

    ! ------------------------
    ! Local variables 
    !   for delta Eddington code
    !   xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
    integer, parameter :: beam = 2
    real,    parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
    real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
            rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th

    ! Compute reflectance and transmittance in a single layer using the two stream approximation 
    !   The code itself is from Lazaros Oreopoulos via Steve Platnick 

    f   = gint**2
    tau = (1 - w0int * f) * tauint
    w0  = (1 - f) * w0int / (1 - w0int * f)
    g   = (gint - f) / (1 - f)

    ! delta-Eddington (Joseph et al. 1976)
    gamma1 =  (7 - w0* (4 + 3 * g)) / 4.0
    gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
    gamma3 =  (2 - 3*g*xmu) / 4.0
    gamma4 =   1 - gamma3

    if (w0int > minConservativeW0) then
      ! Conservative scattering
      if (beam == 1) then
          rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu)))
          ref = rh / (1 + gamma1 * tau)
          tra = 1 - ref       
      else if(beam == 2) then
          ref = gamma1*tau/(1 + gamma1*tau)
          tra = 1 - ref
      endif
    else
      ! Non-conservative scattering
      a1 = gamma1 * gamma4 + gamma2 * gamma3
      a2 = gamma1 * gamma3 + gamma2 * gamma4

      rk = sqrt(gamma1**2 - gamma2**2)
      
      r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
      r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
      r3 = 2 * rk *(gamma3 - a2 * xmu)
      r4 = (1 - (rk * xmu)**2) * (rk + gamma1)
      r5 = (1 - (rk * xmu)**2) * (rk - gamma1)
      
      t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
      t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
      t3 = 2 * rk * (gamma4 + a1 * xmu)
      t4 = r4
      t5 = r5

      beta = -r5 / r4         
      
      e1 = min(rk * tau, 500.) 
      e2 = min(tau / xmu, 500.) 
      
      if (beam == 1) then
         den = r4 * exp(e1) + r5 * exp(-e1)
         ref  = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den
         den = t4 * exp(e1) + t5 * exp(-e1)
         th  = exp(-e2)
         tra = th-th*w0*(t1*exp(e1)-t2*exp(-e1)-t3*exp(e2))/den
      elseif (beam == 2) then
         ef1 = exp(-e1)
         ef2 = exp(-2*e1)
         ref = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2))
         tra = (2*rk*ef1)/((rk+gamma1)*(1-beta*ef2))
      endif
    end if
  end subroutine two_stream
  ! --------------------------------------------------
  elemental function two_stream_reflectance(tauint, gint, w0int) 
    real, intent(in) :: tauint, gint, w0int
    real             :: two_stream_reflectance

    ! Compute reflectance in a single layer using the two stream approximation 
    !   The code itself is from Lazaros Oreopoulos via Steve Platnick 

    ! ------------------------
    ! Local variables 
    !   for delta Eddington code
    !   xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
    integer, parameter :: beam = 2
    real,    parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
    real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
            rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
    ! ------------------------


    f   = gint**2
    tau = (1 - w0int * f) * tauint
    w0  = (1 - f) * w0int / (1 - w0int * f)
    g   = (gint - f) / (1 - f)

    ! delta-Eddington (Joseph et al. 1976)
    gamma1 =  (7 - w0* (4 + 3 * g)) / 4.0
    gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
    gamma3 =  (2 - 3*g*xmu) / 4.0
    gamma4 =   1 - gamma3

    if (w0int > minConservativeW0) then
      ! Conservative scattering
      if (beam == 1) then
          rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu)))
          two_stream_reflectance = rh / (1 + gamma1 * tau)
      elseif (beam == 2) then
          two_stream_reflectance = gamma1*tau/(1 + gamma1*tau)
      endif
        
    else    !

        ! Non-conservative scattering
         a1 = gamma1 * gamma4 + gamma2 * gamma3
         a2 = gamma1 * gamma3 + gamma2 * gamma4

         rk = sqrt(gamma1**2 - gamma2**2)
         
         r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
         r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
         r3 = 2 * rk *(gamma3 - a2 * xmu)
         r4 = (1 - (rk * xmu)**2) * (rk + gamma1)
         r5 = (1 - (rk * xmu)**2) * (rk - gamma1)
         
         t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
         t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
         t3 = 2 * rk * (gamma4 + a1 * xmu)
         t4 = r4
         t5 = r5

         beta = -r5 / r4         
         
         e1 = min(rk * tau, 500.) 
         e2 = min(tau / xmu, 500.) 
         
         if (beam == 1) then
           den = r4 * exp(e1) + r5 * exp(-e1)
           two_stream_reflectance  = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den
         elseif (beam == 2) then
           ef1 = exp(-e1)
           ef2 = exp(-2*e1)
           two_stream_reflectance = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2))
         endif
           
      end if
  end function two_stream_reflectance 
  ! --------------------------------------------
    pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot)      
      real,    dimension(:), intent(in)  :: Refl,     Tran
      real,                  intent(out) :: Refl_tot, Tran_tot

      ! Use adding/doubling formulas to compute total reflectance and transmittance from layer values

      integer :: i
      real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative
      
      Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)    
      
      do i=2, size(Refl)
          ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
          Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i))
          Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1 - Refl_cumulative(i-1) * Refl(i))
      end do
      
      Refl_tot = Refl_cumulative(size(Refl))
      Tran_tot = Tran_cumulative(size(Refl))

    end subroutine adding_doubling
  ! --------------------------------------------------
  subroutine complain_and_die(message) 
    character(len = *), intent(in) :: message
    
    write(6, *) "Failure in MODIS simulator" 
    write(6, *)  trim(message) 
    stop
  end subroutine complain_and_die
  !------------------------------------------------------------------------------------------------
end module mod_modis_sim
