source: LMDZ6/trunk/libf/phylmd/cospv2/modis_simulator.F90 @ 4096

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Copyright (c) 2015, Regents of the University of Colorado
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are 
! permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of 
!    conditions and the following disclaimer.
!
! 2. 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.
!
! 3. Neither the name of the copyright holder 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 HOLDER 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
! May 2015:      Dustin Swales - Modified for COSPv2.0
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! 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_CONFIG, only: R_UNDEF,modis_histTau,modis_histPres,numMODISTauBins,       &
                             numMODISPresBins,numMODISReffIceBins,numMODISReffLiqBins,   &
                             modis_histReffIce,modis_histReffLiq
  USE COSP_KINDS,      ONLY: wp
  use MOD_COSP_STATS,  ONLY: hist2D

  implicit none
  ! ##########################################################################
  ! Retrieval parameters
   integer, parameter :: &
        num_trial_res = 15              ! Increase to make the linear pseudo-retrieval of size more accurate
   
   real(wp) :: &
       min_OpticalThickness,          & ! Minimum detectable optical thickness
       CO2Slicing_PressureLimit,      & ! Cloud with higher pressures use thermal methods, units Pa
       CO2Slicing_TauLimit,           & ! How deep into the cloud does CO2 slicing see? 
       phase_TauLimit,                & ! How deep into the cloud does the phase detection see?
       size_TauLimit,                 & ! Depth of the re retreivals
       phaseDiscrimination_Threshold, & ! What fraction of total extincton needs to be in a single
                                        ! category to make phase discrim. work? 
       re_fill,                       & !
       re_water_min,                  & ! Minimum effective radius (liquid)
       re_water_max,                  & ! Maximum effective radius (liquid)
       re_ice_min,                    & ! Minimum effective radius (ice)
       re_ice_max,                    & ! Minimum effective radius (ice)
       highCloudPressureLimit,        & ! High cloud pressure limit (Pa)
       lowCloudPressureLimit            ! Low cloud pressure limit (Pa)
  integer :: &
       phaseIsNone,                   & !
       phaseIsLiquid,                 & !
       phaseIsIce,                    & !
       phaseIsUndetermined              !
 
  real(wp),dimension(num_trial_res) :: &
       trial_re_w, & ! Near-IR optical params vs size for retrieval scheme (liquid)
       trial_re_i    ! Near-IR optical params vs size for retrieval scheme (ice)
  real(wp),dimension(num_trial_res) :: &
       g_w,        & ! Assymettry parameter for size retrieval (liquid)
       g_i,        & ! Assymettry parameter for size retrieval (ice)
       w0_w,       & ! Single-scattering albedo for size retrieval (liquid)
       w0_i          ! Single-scattering albedo for size retrieval (ice)
  ! Algorithmic parameters
  real(wp),parameter :: &
     ice_density = 0.93_wp ! Liquid density is 1. 
      
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_subcolumn(nSubCols, nLevels, pressureLevels, optical_thickness,       & 
                         tauLiquidFraction, g, w0,isccpCloudTopPressure,                 &
                         retrievedPhase, retrievedCloudTopPressure,                      &
                         retrievedTau,   retrievedSize)

    ! INPUTS
    integer,intent(in) :: &
         nSubCols,                  & ! Number of subcolumns
         nLevels                      ! Number of levels         
    real(wp),dimension(nLevels+1),intent(in) :: &
         pressureLevels               ! Gridmean pressure at layer edges (Pa)                  
    real(wp),dimension(nSubCols,nLevels),intent(in) :: &
         optical_thickness,         & ! Subcolumn optical thickness @ 0.67 microns.
         tauLiquidFraction,         & ! Liquid water fraction
         g,                         & ! Subcolumn assymetry parameter  
         w0                           ! Subcolumn single-scattering albedo 
    real(wp),dimension(nSubCols),intent(in) :: &
         isccpCloudTopPressure        ! ISCCP retrieved cloud top pressure (Pa)

    ! OUTPUTS
    integer, dimension(nSubCols), intent(inout) :: &
         retrievedPhase               ! MODIS retrieved phase (liquid/ice/other)              
    real(wp),dimension(nSubCols), intent(inout) :: &
         retrievedCloudTopPressure, & ! MODIS retrieved CTP (Pa)
         retrievedTau,              & ! MODIS retrieved optical depth (unitless)              
         retrievedSize                ! MODIS retrieved particle size (microns)              

    ! LOCAL VARIABLES
    logical, dimension(nSubCols)      :: &
         cloudMask
    real(wp)                          :: &
         integratedLiquidFraction,       &
         obs_Refl_nir
    real(wp),dimension(num_trial_res) :: &
         predicted_Refl_nir
    integer                           :: &
         i

    ! ########################################################################################
    !                           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(1:nSubCols) = sum(optical_thickness(1:nSubCols,1:nLevels), dim = 2)

    ! ########################################################################################
    !                                 Cloud detection
    ! does optical thickness exceed detection threshold? 
    ! ########################################################################################
    cloudMask = retrievedTau(1:nSubCols) >= min_OpticalThickness
    
    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(nLevels,(/ 0._wp, optical_thickness(i,1:nLevels) /), &
                                                            pressureLevels(1:nLevels),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(nLevels,optical_thickness(i,1:nLevels),       &
                                                          tauLiquidFraction(i, 1:nLevels), &
                                                          phase_TauLimit)
          if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then 
             retrievedPhase(i) = phaseIsLiquid
          else if (integratedLiquidFraction <= 1._wp- phaseDiscrimination_Threshold) then 
             retrievedPhase(i) = phaseIsIce
          else 
             retrievedPhase(i) = phaseIsUndetermined
          end if
        
          ! ##################################################################################
          !                                 Size determination 
          ! ##################################################################################
          
          ! Compute observed reflectance
          obs_Refl_nir = compute_toa_reflectace(nLevels,optical_thickness(i,1:nLevels), g(i,1:nLevels), w0(i,1:nLevels))

          ! Compute predicted reflectance
          if(any(retrievedPhase(i) == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then 
             if (retrievedPhase(i) == phaseIsLiquid .OR. retrievedPhase(i) == phaseIsUndetermined) then
                predicted_Refl_nir(1:num_trial_res) = two_stream_reflectance(retrievedTau(i), &
                     g_w(1:num_trial_res), w0_w(1:num_trial_res))
                retrievedSize(i) = 1.0e-06_wp*interpolate_to_min(trial_re_w(1:num_trial_res), &
                     predicted_Refl_nir(1:num_trial_res), obs_Refl_nir)
             else
                predicted_Refl_nir(1:num_trial_res) = two_stream_reflectance(retrievedTau(i), &
                     g_i(1:num_trial_res), w0_i(1:num_trial_res))
                retrievedSize(i) = 1.0e-06_wp*interpolate_to_min(trial_re_i(1:num_trial_res), &
                     predicted_Refl_nir(1:num_trial_res), obs_Refl_nir)
             endif
          else 
             retrievedSize(i) = re_fill
          endif
       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(1:nSubCols) < 0.).and.(retrievedSize(1:nSubCols) /= R_UNDEF)) &
         retrievedSize(1:nSubCols) = 1.0e-06_wp*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(1:nSubCols) .and. retrievedCloudTopPressure(1:nSubCols) > CO2Slicing_PressureLimit) &
         retrievedCloudTopPressure(1:nSubCols) = isccpCloudTopPressure! * 100._wp
    
  end subroutine modis_subcolumn

  ! ########################################################################################
  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(nPoints, nSubCols) ::  &
         phase                             
    real(wp),intent(in),dimension(nPoints, nSubCols) ::  &
         cloud_top_pressure,                &
         optical_thickness,                 &
         particle_size
 
    ! OUTPUTS 
    real(wp),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(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISPresBins) :: &
         Optical_Thickness_vs_Cloud_Top_Pressure
    real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISReffIceBins) :: &    
         Optical_Thickness_vs_ReffIce
    real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISReffLiqBins) :: &    
         Optical_Thickness_vs_ReffLiq         

    ! LOCAL VARIABLES
    real(wp), parameter :: &
         LWP_conversion = 2._wp/3._wp * 1000._wp ! MKS units  
    integer :: j
    logical, dimension(nPoints,nSubCols) :: &
         cloudMask,      &
         waterCloudMask, &
         iceCloudMask,   &
         validRetrievalMask
    real(wp),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 column amounts.
    ! ########################################################################################
    where(Cloud_Fraction_Total_Mean(1:nPoints) > 0)
       Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask,      dim = 2) / &
            Cloud_Fraction_Total_Mean(1:nPoints)
       Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
            dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
    elsewhere
       Optical_Thickness_Total_Mean      = R_UNDEF
       Optical_Thickness_Total_MeanLog10 = R_UNDEF
    endwhere
    where(Cloud_Fraction_Water_Mean(1:nPoints) > 0)
       Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
            Cloud_Fraction_Water_Mean(1:nPoints)
       Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
            mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
       Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
            dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
       Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
            Cloud_Fraction_Water_Mean(1:nPoints)
    elsewhere
       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(1:nPoints) > 0)
       Optical_Thickness_Ice_Mean(1:nPoints)   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / &
            Cloud_Fraction_Ice_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) 
       Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
            dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
       Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask,   dim = 2) / &
            Cloud_Fraction_Ice_Mean(1:nPoints)    
    elsewhere
       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
    Cloud_Top_Pressure_Total_Mean  = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
                                     max(1, count(cloudMask, dim = 2))

    ! ########################################################################################
    ! 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
    
    ! ########################################################################################
    ! Joint histograms
    ! ########################################################################################
    ! 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._wp
          ctpWRK(j,1:nSubCols) = -999._wp
       endwhere
       ! Joint histogram of tau/CTP
       call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
                   modis_histTau,numMODISTauBins,&
                   modis_histPres,numMODISPresBins,&
                   Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numMODISTauBins,1:numMODISPresBins))
       ! Joint histogram of tau/ReffICE
       call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols,               &
                   modis_histTau,numMODISTauBins,modis_histReffIce,         &
                   numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numMODISTauBins,1:numMODISReffIceBins))
       ! Joint histogram of tau/ReffLIQ
       call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols,               &
                   modis_histTau,numMODISTauBins,modis_histReffLiq,         &
                   numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numMODISTauBins,1:numMODISReffLiqBins))                   

    enddo   
    Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numMODISTauBins,1:numMODISPresBins) = &
         Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numMODISTauBins,1:numMODISPresBins)/nSubCols
    Optical_Thickness_vs_ReffIce(1:nPoints,1:numMODISTauBins,1:numMODISReffIceBins) = &
         Optical_Thickness_vs_ReffIce(1:nPoints,1:numMODISTauBins,1:numMODISReffIceBins)/nSubCols
    Optical_Thickness_vs_ReffLiq(1:nPoints,1:numMODISTauBins,1:numMODISReffLiqBins) = &
         Optical_Thickness_vs_ReffLiq(1:nPoints,1:numMODISTauBins,1:numMODISReffLiqBins)/nSubCols 
                 

    ! Unit conversion
    where(Optical_Thickness_vs_Cloud_Top_Pressure /= R_UNDEF) &
      Optical_Thickness_vs_Cloud_Top_Pressure = Optical_Thickness_vs_Cloud_Top_Pressure*100._wp
    where(Optical_Thickness_vs_ReffIce /= R_UNDEF) Optical_Thickness_vs_ReffIce = Optical_Thickness_vs_ReffIce*100._wp
    where(Optical_Thickness_vs_ReffLiq /= R_UNDEF) Optical_Thickness_vs_ReffLiq = Optical_Thickness_vs_ReffLiq*100._wp
    where(Cloud_Fraction_Total_Mean /= R_UNDEF) Cloud_Fraction_Total_Mean = Cloud_Fraction_Total_Mean*100._wp
    where(Cloud_Fraction_Water_Mean /= R_UNDEF) Cloud_Fraction_Water_Mean = Cloud_Fraction_Water_Mean*100._wp
    where(Cloud_Fraction_Ice_Mean /= R_UNDEF)   Cloud_Fraction_Ice_Mean = Cloud_Fraction_Ice_Mean*100._wp
    where(Cloud_Fraction_High_Mean /= R_UNDEF)  Cloud_Fraction_High_Mean = Cloud_Fraction_High_Mean*100._wp
    where(Cloud_Fraction_Mid_Mean /= R_UNDEF)   Cloud_Fraction_Mid_Mean = Cloud_Fraction_Mid_Mean*100._wp
    where(Cloud_Fraction_Low_Mean /= R_UNDEF)   Cloud_Fraction_Low_Mean = Cloud_Fraction_Low_Mean*100._wp

  end subroutine modis_column

  ! ########################################################################################
  function cloud_top_pressure(nLevels,tauIncrement, pressure, tauLimit) 
    ! INPUTS
    integer, intent(in)                    :: nLevels
    real(wp),intent(in),dimension(nLevels) :: tauIncrement, pressure
    real(wp),intent(in)                    :: tauLimit
    ! OUTPUTS
    real(wp)                               :: cloud_top_pressure
    ! LOCAL VARIABLES
    real(wp)                               :: deltaX, totalTau, totalProduct
    integer                                :: i 
    
    ! Find the extinction-weighted pressure. Assume that pressure varies linearly between 
    !   layers and use the trapezoidal rule.
    totalTau = 0._wp; totalProduct = 0._wp
    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._wp * tauIncrement(i)) 
      else
        totalTau =     totalTau     + tauIncrement(i) 
        totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2._wp
      end if 
      if(totalTau >= tauLimit) exit
    end do 

    if (totalTau > 0._wp) then
       cloud_top_pressure = totalProduct/totalTau
    else
       cloud_top_pressure = 0._wp
    endif
    
  end function cloud_top_pressure

  ! ########################################################################################
  function weight_by_extinction(nLevels,tauIncrement, f, tauLimit) 
    ! INPUTS
    integer, intent(in)                    :: nLevels
    real(wp),intent(in),dimension(nLevels) :: tauIncrement, f
    real(wp),intent(in)                    :: tauLimit
    ! OUTPUTS
    real(wp)                               :: weight_by_extinction
    ! LOCAL VARIABLES
    real(wp)                               :: deltaX, totalTau, totalProduct
    integer                                :: i 
    
    ! Find the extinction-weighted value of f(tau), assuming constant f within each layer
    totalTau = 0._wp; totalProduct = 0._wp
    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 

    if (totalTau > 0._wp) then
       weight_by_extinction = totalProduct/totalTau
    else
       weight_by_extinction = 0._wp
    endif
    
  end function weight_by_extinction

  ! ########################################################################################
  pure function interpolate_to_min(x, y, yobs)
    ! INPUTS
    real(wp),intent(in),dimension(num_trial_res) :: x, y 
    real(wp),intent(in)                          :: yobs
    ! OUTPUTS
    real(wp)                                     :: interpolate_to_min
    ! LOCAL VARIABLES
    real(wp), dimension(num_trial_res)           :: diff
    integer                                      :: nPoints, minDiffLoc, lowerBound, upperBound
    
    ! 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
 
    nPoints = size(y)
    diff(1:num_trial_res) = y(1:num_trial_res) - 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_old (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

    ! INPUTS
    integer, intent(in) :: phase
    real(wp),intent(in) :: re
    ! OUTPUTS
    real(wp)            :: get_g_nir_old 
    ! LOCAL VARIABLES(parameters)
    real(wp), 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_old = fit_to_quadratic(re, small_water_coefficients)
        if(re < re_water_min) get_g_nir_old = fit_to_quadratic(re_water_min, small_water_coefficients)
      else
        get_g_nir_old = fit_to_quadratic(re,   big_water_coefficients)
        if(re > re_water_max) get_g_nir_old = fit_to_quadratic(re_water_max, big_water_coefficients)
      end if 
    else
      get_g_nir_old = fit_to_quadratic(re, ice_coefficients)
      if(re < re_ice_min) get_g_nir_old = fit_to_quadratic(re_ice_min, ice_coefficients)
      if(re > re_ice_max) get_g_nir_old = fit_to_quadratic(re_ice_max, ice_coefficients)
    end if 
    
  end function get_g_nir_old

  ! ########################################################################################
  elemental function get_ssa_nir_old (phase, re)
    ! 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
    
    ! INPUTS
    integer, intent(in) :: phase
    real(wp),intent(in) :: re
    ! OUTPUTS
    real(wp)            :: get_ssa_nir_old
    ! LOCAL VARIABLES (parameters)
    real(wp), dimension(4), parameter :: ice_coefficients   = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /)
    real(wp), dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /) 
    
    ! approx. fits from MODIS Collection 5 LUT scattering calculations
    if(phase == phaseIsLiquid) then
       get_ssa_nir_old = fit_to_quadratic(re, water_coefficients)
       if(re < re_water_min) get_ssa_nir_old = fit_to_quadratic(re_water_min, water_coefficients)
       if(re > re_water_max) get_ssa_nir_old = fit_to_quadratic(re_water_max, water_coefficients)
    else
       get_ssa_nir_old = fit_to_cubic(re, ice_coefficients)
       if(re < re_ice_min) get_ssa_nir_old = fit_to_cubic(re_ice_min, ice_coefficients)
       if(re > re_ice_max) get_ssa_nir_old = fit_to_cubic(re_ice_max, ice_coefficients)
    end if
    
  end function get_ssa_nir_old
  
  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(wp),    intent(in) :: re
    real(wp) :: get_g_nir 

    real(wp), 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(wp), 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(wp),    intent(in) :: re
        real(wp)                :: 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(wp), dimension(4), parameter :: ice_coefficients   = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
        real(wp), 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) 
    ! INPUTS
    real(wp),               intent(in) :: x
    real(wp), dimension(4), intent(in) :: coefficients
    ! OUTPUTS
    real(wp)                           :: 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) 
    ! INPUTS
    real(wp),               intent(in) :: x
    real(wp), dimension(3), intent(in) :: coefficients
    ! OUTPUTS
    real(wp)                           :: 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(nLevels,tau, g, w0)
    ! This wrapper reports reflectance only and strips out non-cloudy elements from the 
    ! calculation
    
    ! INPUTS
    integer,intent(in)                     :: nLevels
    real(wp),intent(in),dimension(nLevels) :: tau, g, w0
    ! OUTPUTS
    real(wp)                               :: compute_toa_reflectace
    ! LOCAL VARIABLES
    logical, dimension(nLevels)                   :: cloudMask
    integer, dimension(count(tau(1:nLevels) > 0)) :: cloudIndicies
    real(wp),dimension(count(tau(1:nLevels) > 0)) :: Refl,Trans
    real(wp)                                      :: Refl_tot, Trans_tot
    integer                                       :: i

    cloudMask(1:nLevels) = tau(1:nLevels) > 0. 
    cloudIndicies = pack((/ (i, i = 1, nLevels) /), 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(count(tau(1:nLevels) > 0),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) 
    ! Compute reflectance in a single layer using the two stream approximation 
    !   The code itself is from Lazaros Oreopoulos via Steve Platnick 
    ! INPUTS
    real(wp), intent(in)  :: tauint, gint, w0int
    ! OUTPUTS
    real(wp), intent(out) :: ref, tra
    ! 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(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
    real(wp) :: 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._wp - w0int * f) * tauint
    w0  = (1._wp - f) * w0int / (1._wp - w0int * f)
    g   = (gint - f) / (1._wp - f)

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

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

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

      beta = -r5 / r4         
  
      e1 = min(rk * tau, 500._wp) 
      e2 = min(tau / xmu, 500._wp) 
      
      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._wp-ef2))/((rk+gamma1)*(1._wp-beta*ef2))
         tra = (2._wp*rk*ef1)/((rk+gamma1)*(1._wp-beta*ef2))
      endif
    end if
  end subroutine two_stream

  ! ########################################################################################
  elemental function two_stream_reflectance(tauint, gint, w0int)
    ! Compute reflectance in a single layer using the two stream approximation 
    !   The code itself is from Lazaros Oreopoulos via Steve Platnick 
    
    ! INPUTS
    real(wp), intent(in) :: tauint, gint, w0int
    ! OUTPUTS
    real(wp)             :: two_stream_reflectance
    ! 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(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
    real(wp) :: 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._wp - w0int * f) * tauint
    w0  = (1._wp - f) * w0int / (1._wp - w0int * f)
    g   = (gint - f) / (1._wp - f)

    ! delta-Eddington (Joseph et al. 1976)
    gamma1 =  (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp
    gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp
    gamma3 =  (2._wp - 3._wp*g*xmu) / 4._wp
    gamma4 =   1._wp - 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._wp + gamma1 * tau)
      elseif (beam == 2) then
          two_stream_reflectance = gamma1*tau/(1._wp + 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._wp - rk * xmu) * (a2 + rk * gamma3)
         r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3)
         r3 = 2._wp * rk *(gamma3 - a2 * xmu)
         r4 = (1._wp - (rk * xmu)**2) * (rk + gamma1)
         r5 = (1._wp - (rk * xmu)**2) * (rk - gamma1)
         
         t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4)
         t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4)
         t3 = 2._wp * rk * (gamma4 + a1 * xmu)
         t4 = r4
         t5 = r5

         beta = -r5 / r4         
         
         e1 = min(rk * tau, 500._wp) 
         e2 = min(tau / xmu, 500._wp) 
         
         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._wp-ef2))/((rk+gamma1)*(1._wp-beta*ef2))
         endif
           
      end if
  end function two_stream_reflectance 

  ! ########################################################################################
  pure subroutine adding_doubling (npts,Refl, Tran, Refl_tot, Tran_tot)      
    ! Use adding/doubling formulas to compute total reflectance and transmittance from 
    ! layer values
    
    ! INPUTS
    integer,intent(in)                  :: npts
    real(wp),intent(in),dimension(npts) :: Refl,Tran
    ! OUTPUTS
    real(wp),intent(out)                :: Refl_tot, Tran_tot
    ! LOCAL VARIABLES
    integer :: i
    real(wp), dimension(npts) :: Refl_cumulative, Tran_cumulative
    
    Refl_cumulative(1) = Refl(1)
    Tran_cumulative(1) = Tran(1)    
    
    do i=2, npts
       ! 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._wp - Refl_cumulative(i-1) * Refl(i))
       Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1._wp - Refl_cumulative(i-1) * Refl(i))
    end do
    
    Refl_tot = Refl_cumulative(size(Refl))
    Tran_tot = Tran_cumulative(size(Refl))
    
  end subroutine adding_doubling

end module mod_modis_sim