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