Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

mod_modis_sim.f90 @ 5428

Last change on this file since 5428 was 5268, checked in by abarral, 2 months ago
.f90 <-> .F90 depending on cpp key use
File size: 64.2 KB

Rev	Line
[2432]	1	! (c) 2009-2010, Regents of the Unversity of Colorado
	2	! Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences
	3	! All rights reserved.
[2713]	4	! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
[2432]	5	! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $
	6	!
	7	! Redistribution and use in source and binary forms, with or without modification, are permitted
	8	! provided that the following conditions are met:
	9	!
	10	! * Redistributions of source code must retain the above copyright notice, this list
	11	! of conditions and the following disclaimer.
	12	! * Redistributions in binary form must reproduce the above copyright notice, this list
	13	! of conditions and the following disclaimer in the documentation and/or other materials
	14	! provided with the distribution.
	15	! * Neither the name of the Met Office nor the names of its contributors may be used
	16	! to endorse or promote products derived from this software without specific prior written
	17	! permission.
	18	!
	19	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
	20	! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
	21	! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
	22	! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
	23	! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	24	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
	25	! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
	26	! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	27	!
	28
	29	!
	30	! History:
	31	! May 2009 - Robert Pincus - Initial version
	32	! June 2009 - Steve Platnick and Robert Pincus - Simple radiative transfer for size retrievals
	33	! August 2009 - Robert Pincus - Consistency and bug fixes suggested by Rick Hemler (GFDL)
	34	! November 2009 - Robert Pincus - Bux fixes and speed-ups after experience with Rick Hemler using AM2 (GFDL)
	35	! January 2010 - Robert Pincus - Added high, middle, low cloud fractions
	36	!
	37
	38	!
	39	! Notes on using the MODIS simulator:
	40	! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 microns, or
	41	! optical thickness at 0.67 microns and ice- and liquid-water contents (in consistent units of
	42	! your choosing)
	43	! *) Required input also includes the optical thickness and cloud top pressure
	44	! derived from the ISCCP simulator run with parameter top_height = 1.
	45	! *) Cloud particle sizes are specified as radii, measured in meters, though within the module we
	46	! use units of microns. Where particle sizes are outside the bounds used in the MODIS retrieval
	47	! libraries (parameters re_water_min, re_ice_min, etc.) the simulator returns missing values (re_fill)
	48
	49	!
	50	! When error conditions are encountered this code calls the function complain_and_die, supplied at the
	51	! bottom of this module. Users probably want to replace this with something more graceful.
	52	!
	53	module mod_modis_sim
[3241]	54	USE MOD_COSP_TYPES, only: R_UNDEF, ok_debug_cosp
[2432]	55	implicit none
	56	! ------------------------------
	57	! Algorithmic parameters
	58	!
	59
	60	real, parameter :: ice_density = 0.93 ! liquid density is 1.
	61	!
	62	! Retrieval parameters
	63	!
	64	real, parameter :: min_OpticalThickness = 0.3, & ! Minimum detectable optical thickness
	65	CO2Slicing_PressureLimit = 700. * 100., & ! Cloud with higher pressures use thermal methods, units Pa
	66	CO2Slicing_TauLimit = 1., & ! How deep into the cloud does CO2 slicing see?
	67	phase_TauLimit = 1., & ! How deep into the cloud does the phase detection see?
	68	size_TauLimit = 2., & ! Depth of the re retreivals
	69	phaseDiscrimination_Threshold = 0.7 ! What fraction of total extincton needs to be
	70	! in a single category to make phase discrim. work?
	71	real, parameter :: re_fill= -999.
	72	integer, parameter :: phaseIsNone = 0, phaseIsLiquid = 1, phaseIsIce = 2, phaseIsUndetermined = 3
	73
	74	logical, parameter :: useSimpleReScheme = .false.
	75	!
	76	! These are the limits of the libraries for the MODIS collection 5 algorithms
	77	! They are also the limits used in the fits for g and w0
	78	!
	79	real, parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90.
	80	integer, parameter :: num_trial_res = 15 ! increase to make the linear pseudo-retrieval of size more accurate
[2713]	81	! DJS2015: Remove unused parameter
	82	! logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations?
	83	! DJS2015 END
[2432]	84	!
	85	! Precompute near-IR optical params vs size for retrieval scheme
	86	!
	87	integer, private :: i
	88	real, dimension(num_trial_res), parameter :: &
	89	trial_re_w = re_water_min + (re_water_max - re_water_min)/(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /), &
	90	trial_re_i = re_ice_min + (re_ice_max - re_ice_min) /(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /)
	91
	92	! Can't initialze these during compilation, but do in before looping columns in retrievals
	93	real, dimension(num_trial_res) :: g_w, g_i, w0_w, w0_i
	94	! ------------------------------
	95	! Bin boundaries for the joint optical thickness/cloud top pressure histogram
	96	!
	97	integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7
	98
	99	real, private :: dummy_real
	100	real, dimension(numTauHistogramBins + 1), parameter :: &
	101	tauHistogramBoundaries = (/ min_OpticalThickness, 1.3, 3.6, 9.4, 23., 60., 10000. /)
	102	real, dimension(numPressureHistogramBins + 1), parameter :: & ! Units Pa
	103	pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /)
	104	real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100.
	105	!
	106	! For output - nominal bin centers and bin boundaries. On output pressure bins are highest to lowest.
	107	!
	108	integer, private :: k, l
	109	real, parameter, dimension(2, numTauHistogramBins) :: &
	110	nominalTauHistogramBoundaries = &
	111	reshape(source = (/ tauHistogramBoundaries(1), &
	112	((tauHistogramBoundaries(k), l = 1, 2), k = 2, numTauHistogramBins), &
	113	100000. /), &
	114	shape = (/2, numTauHistogramBins /) )
	115	real, parameter, dimension(numTauHistogramBins) :: &
	116	nominalTauHistogramCenters = (nominalTauHistogramBoundaries(1, :) + &
	117	nominalTauHistogramBoundaries(2, :) ) / 2.
	118
	119	real, parameter, dimension(2, numPressureHistogramBins) :: &
	120	nominalPressureHistogramBoundaries = &
	121	reshape(source = (/ 100000., &
	122	((pressureHistogramBoundaries(k), l = 1, 2), k = numPressureHistogramBins, 2, -1), &
	123	0. /), &
	124	shape = (/2, numPressureHistogramBins /) )
	125	real, parameter, dimension(numPressureHistogramBins) :: &
	126	nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + &
	127	nominalPressureHistogramBoundaries(2, :) ) / 2.
[2713]	128	! DJS2015 START: Add bin descriptions for joint-histograms of partice-sizes and optical depth. This is
	129	! identical to what is done in COSPv.2.0.0 for histogram bin initialization.
	130	integer :: j
	131	integer,parameter :: &
	132	numMODISReffLiqBins = 6, & ! Number of bins for tau/ReffLiq joint-histogram
	133	numMODISReffIceBins = 6 ! Number of bins for tau/ReffICE joint-histogram
	134	real,parameter,dimension(numMODISReffLiqBins+1) :: &
	135	reffLIQ_binBounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/)
	136	real,parameter,dimension(numMODISReffIceBins+1) :: &
	137	reffICE_binBounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/)
	138	real,parameter,dimension(2,numMODISReffIceBins) :: &
	139	reffICE_binEdges = reshape(source=(/reffICE_binBounds(1),((reffICE_binBounds(k), &
	140	l=1,2),k=2,numMODISReffIceBins),reffICE_binBounds(numMODISReffIceBins+1)/), &
	141	shape = (/2,numMODISReffIceBins/))
	142	real,parameter,dimension(2,numMODISReffLiqBins) :: &
	143	reffLIQ_binEdges = reshape(source=(/reffLIQ_binBounds(1),((reffLIQ_binBounds(k), &
	144	l=1,2),k=2,numMODISReffLiqBins),reffLIQ_binBounds(numMODISReffIceBins+1)/), &
	145	shape = (/2,numMODISReffLiqBins/))
	146	real,parameter,dimension(numMODISReffIceBins) :: &
	147	reffICE_binCenters = (reffICE_binEdges(1,:)+reffICE_binEdges(2,:))/2.
	148	real,parameter,dimension(numMODISReffLiqBins) :: &
	149	reffLIQ_binCenters = (reffLIQ_binEdges(1,:)+reffLIQ_binEdges(2,:))/2.
	150	! DJS2015 END
	151
[2432]	152	! ------------------------------
	153	! There are two ways to call the MODIS simulator:
	154	! 1) Provide total optical thickness and liquid/ice water content and we'll partition tau in
	155	! subroutine modis_L2_simulator_oneTau, or
	156	! 2) Provide ice and liquid optical depths in each layer
	157	!
	158	interface modis_L2_simulator
	159	module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus
	160	end interface
	161	contains
	162	!------------------------------------------------------------------------------------------------
	163	! MODIS simulator using specified liquid and ice optical thickness in each layer
	164	!
	165	! Note: this simulator operates on all points; to match MODIS itself night-time
	166	! points should be excluded
	167	!
	168	! Note: the simulator requires as input the optical thickness and cloud top pressure
	169	! derived from the ISCCP simulator run with parameter top_height = 1.
	170	! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing
	171	! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that
	172	! alogrithm in this simulator we simply report the values from the ISCCP simulator.
	173	!
	174	subroutine modis_L2_simulator_twoTaus( &
	175	temp, pressureLayers, pressureLevels, &
	176	liquid_opticalThickness, ice_opticalThickness, &
	177	waterSize, iceSize, &
	178	isccpTau, isccpCloudTopPressure, &
	179	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
	180
	181	! Grid-mean quantities at layer centers, starting at the model top
	182	! dimension nLayers
	183	real, dimension(:), intent(in ) :: temp, & ! Temperature, K
	184	pressureLayers, & ! Pressure, Pa
	185	pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1)
	186	! Sub-column quantities
	187	! dimension nSubcols, nLayers
	188	real, dimension(:, :), intent(in ) :: liquid_opticalThickness, & ! Layer optical thickness @ 0.67 microns due to liquid
	189	ice_opticalThickness ! ditto, due to ice
	190	real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns
	191	iceSize ! Cloud ice effective radius, microns
	192
	193	! Cloud properties retrieved from ISCCP using top_height = 1
	194	! dimension nSubcols
	195	real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness
	196	isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa)
	197
	198	! Properties retrieved by MODIS
	199	! dimension nSubcols
	200	integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer, defined in module header
	201	real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
	202	retrievedTau, & ! unitless
	203	retrievedSize ! microns
	204	! ---------------------------------------------------
	205	! Local variables
	206	logical, dimension(size(retrievedTau)) :: cloudMask
	207	real, dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal
	208	real :: integratedLiquidFraction
	209	integer :: i, nSubcols, nLevels
	210
	211	! ---------------------------------------------------
	212	nSubcols = size(liquid_opticalThickness, 1)
	213	nLevels = size(liquid_opticalThickness, 2)
	214
	215	!
	216	! Initial error checks
	217	!
	218	if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), &
	219	size(isccpTau), size(isccpCloudTopPressure), &
	220	size(retrievedPhase), size(retrievedCloudTopPressure), &
	221	size(retrievedTau), size(retrievedSize) /) /= nSubcols )) &
	222	call complain_and_die("Differing number of subcolumns in one or more arrays")
	223
	224	if(any((/ size(ice_opticalThickness, 2), size(waterSize, 2), size(iceSize, 2), &
	225	size(temp), size(pressureLayers), size(pressureLevels)-1 /) /= nLevels )) &
	226	call complain_and_die("Differing number of levels in one or more arrays")
	227
	228
	229	if(any( (/ any(temp <= 0.), any(pressureLayers <= 0.), &
	230	any(liquid_opticalThickness < 0.), &
	231	any(ice_opticalThickness < 0.), &
	232	any(waterSize < 0.), any(iceSize < 0.) /) )) &
	233	call complain_and_die("Input values out of bounds")
	234
	235	! ---------------------------------------------------
	236	!
	237	! Compute the total optical thickness and the proportion due to liquid in each cell
	238	!
	239	where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.)
	240	tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :))
	241	elsewhere
	242	tauLiquidFraction(:, :) = 0.
	243	end where
	244	tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :)
	245
	246	!
	247	! Optical depth retrieval
	248	! This is simply a sum over the optical thickness in each layer
	249	! It should agree with the ISCCP values after min values have been excluded
	250	!
	251	retrievedTau(:) = sum(tauTotal(:, :), dim = 2)
	252
	253	!
	254	! Cloud detection - does optical thickness exceed detection threshold?
	255	!
	256	cloudMask = retrievedTau(:) >= min_OpticalThickness
	257
	258	!
	259	! Initialize initial estimates for size retrievals
	260	!
	261	if(any(cloudMask) .and. .not. useSimpleReScheme) then
	262	g_w(:) = get_g_nir( phaseIsLiquid, trial_re_w(:))
	263	w0_w(:) = get_ssa_nir(phaseIsLiquid, trial_re_w(:))
	264	g_i(:) = get_g_nir( phaseIsIce, trial_re_i(:))
	265	w0_i(:) = get_ssa_nir(phaseIsIce, trial_re_i(:))
	266	end if
	267
	268	do i = 1, nSubCols
	269	if(cloudMask(i)) then
	270	!
	271	! Cloud top pressure determination
	272	! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds
	273	! lower than that.
	274	! For CO2 slicing we report the optical-depth weighted pressure, integrating to a specified
	275	! optical depth
	276	! This assumes linear variation in p between levels. Linear in ln(p) is probably better,
	277	! though we'd need to deal with the lowest pressure gracefully.
	278	!
	279	retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), &
	280	pressureLevels, &
	281	CO2Slicing_TauLimit)
	282
	283
	284	!
	285	! Phase determination - determine fraction of total tau that's liquid
	286	! When ice and water contribute about equally to the extinction we can't tell
	287	! what the phase is
	288	!
	289	integratedLiquidFraction = weight_by_extinction(tauTotal(i, :), &
	290	tauLiquidFraction(i, :), &
	291	phase_TauLimit)
	292	if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then
	293	retrievedPhase(i) = phaseIsLiquid
	294	else if (integratedLiquidFraction <= 1.- phaseDiscrimination_Threshold) then
	295	retrievedPhase(i) = phaseIsIce
	296	else
	297	retrievedPhase(i) = phaseIsUndetermined
	298	end if
	299
	300	!
	301	! Size determination
	302	!
	303	if(useSimpleReScheme) then
	304	! This is the extinction-weighted size considering only the phase we've chosen
	305	!
	306	if(retrievedPhase(i) == phaseIsIce) then
	307	retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :), &
	308	iceSize(i, :), &
	309	(1. - integratedLiquidFraction) * size_TauLimit)
	310
	311	else if(retrievedPhase(i) == phaseIsLiquid) then
	312	retrievedSize(i) = weight_by_extinction(liquid_opticalThickness(i, :), &
	313	waterSize(i, :), &
	314	integratedLiquidFraction * size_TauLimit)
	315
	316	else
	317	retrievedSize(i) = 0.
	318	end if
	319	else
	320	retrievedSize(i) = 1.0e-06*retrieve_re(retrievedPhase(i), retrievedTau(i), &
	321	obs_Refl_nir = compute_nir_reflectance(liquid_opticalThickness(i, :), waterSize(i, :)*1.0e6, &
	322	ice_opticalThickness(i, :), iceSize(i, :)*1.0e6))
	323	end if
	324	else
	325	!
	326	! Values when we don't think there's a cloud.
	327	!
	328	retrievedCloudTopPressure(i) = R_UNDEF
	329	retrievedPhase(i) = phaseIsNone
	330	retrievedSize(i) = R_UNDEF
	331	retrievedTau(i) = R_UNDEF
	332	end if
	333	end do
	334	where((retrievedSize(:) < 0.).and.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill
	335
	336	! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS
	337	! mimics what MODIS does to first order.
	338	! Of course, ISCCP cloud top pressures are in mb.
	339	!
	340	where(cloudMask(:) .and. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) &
	341	retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100.
	342
	343	end subroutine modis_L2_simulator_twoTaus
	344	!------------------------------------------------------------------------------------------------
	345	!
	346	! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents;
	347	! we'll partition this into ice and liquid optical thickness and call the full MODIS simulator
	348	!
	349	subroutine modis_L2_simulator_oneTau( &
	350	temp, pressureLayers, pressureLevels, &
	351	opticalThickness, cloudWater, cloudIce, &
	352	waterSize, iceSize, &
	353	isccpTau, isccpCloudTopPressure, &
	354	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
	355	! Grid-mean quantities at layer centers,
	356	! dimension nLayers
	357	real, dimension(:), intent(in ) :: temp, & ! Temperature, K
	358	pressureLayers, & ! Pressure, Pa
	359	pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1)
	360	! Sub-column quantities
	361	! dimension nLayers, nSubcols
	362	real, dimension(:, :), intent(in ) :: opticalThickness, & ! Layer optical thickness @ 0.67 microns
	363	cloudWater, & ! Cloud water content, arbitrary units
	364	cloudIce ! Cloud water content, same units as cloudWater
	365	real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns
	366	iceSize ! Cloud ice effective radius, microns
	367
	368	! Cloud properties retrieved from ISCCP using top_height = 1
	369	! dimension nSubcols
	370
	371	real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness
	372	isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa)
	373
	374	! Properties retrieved by MODIS
	375	! dimension nSubcols
	376	integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer
	377	real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
	378	retrievedTau, & ! unitless
	379	retrievedSize ! microns (or whatever units
	380	! waterSize and iceSize are supplied in)
	381	! ---------------------------------------------------
	382	! Local variables
	383	real, dimension(size(opticalThickness, 1), size(opticalThickness, 2)) :: &
	384	liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction
	385
[3241]	386	real :: seuil
[2432]	387	! ---------------------------------------------------
	388
[3241]	389	if (ok_debug_cosp) then
	390	seuil=1.e-9
	391	else
	392	seuil=0.0
	393	endif
	394
[2432]	395	where(cloudIce(:, :) <= 0.)
	396	tauLiquidFraction(:, :) = 1.
	397	elsewhere
	398	where (cloudWater(:, :) <= 0.)
	399	tauLiquidFraction(:, :) = 0.
	400	elsewhere
	401	!
	402	! Geometic optics limit - tau as LWP/re (proportional to LWC/re)
	403	!
[3241]	404	! Modif AI 02 2018
	405	tauLiquidFraction(:, :) = (cloudWater(:, :)/max(waterSize(:, :), seuil) ) / &
	406	(cloudWater(:, :)/max(waterSize(:, :), seuil) + &
	407	cloudIce(:, :)/(ice_density * max(iceSize(:, :), seuil)) )
[2432]	408	end where
	409	end where
	410	liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :)
	411	ice_opticalThickness (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :)
	412
	413	call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels, &
	414	liquid_opticalThickness, ice_opticalThickness, &
	415	waterSize, iceSize, &
	416	isccpTau, isccpCloudTopPressure, &
	417	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
	418
	419	end subroutine modis_L2_simulator_oneTau
[2713]	420
	421	! ########################################################################################
	422	subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, &
	423	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
	424	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
	425	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
	426	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
	427	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, &
	428	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
	429	Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
	430
	431	! INPUTS
	432	integer,intent(in) :: &
	433	nPoints, & ! Number of horizontal gridpoints
	434	nSubCols ! Number of subcolumns
	435	integer,intent(in), dimension(:,:) :: &
	436	!ds integer,intent(in), dimension(nPoints, nSubCols) :: &
	437	phase
	438	real,intent(in),dimension(:,:) :: &
	439	!ds real,intent(in),dimension(nPoints, nSubCols) :: &
	440	cloud_top_pressure, &
	441	optical_thickness, &
	442	particle_size
	443
	444	! OUTPUTS
	445	real,intent(inout),dimension(:) :: & !
	446	!ds real,intent(inout),dimension(nPoints) :: & !
	447	Cloud_Fraction_Total_Mean, & !
	448	Cloud_Fraction_Water_Mean, & !
	449	Cloud_Fraction_Ice_Mean, & !
	450	Cloud_Fraction_High_Mean, & !
	451	Cloud_Fraction_Mid_Mean, & !
	452	Cloud_Fraction_Low_Mean, & !
	453	Optical_Thickness_Total_Mean, & !
	454	Optical_Thickness_Water_Mean, & !
	455	Optical_Thickness_Ice_Mean, & !
	456	Optical_Thickness_Total_MeanLog10, & !
	457	Optical_Thickness_Water_MeanLog10, & !
	458	Optical_Thickness_Ice_MeanLog10, & !
	459	Cloud_Particle_Size_Water_Mean, & !
	460	Cloud_Particle_Size_Ice_Mean, & !
	461	Cloud_Top_Pressure_Total_Mean, & !
	462	Liquid_Water_Path_Mean, & !
	463	Ice_Water_Path_Mean !
	464	real,intent(inout),dimension(:,:,:) :: &
	465	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: &
	466	Optical_Thickness_vs_Cloud_Top_Pressure
	467	real,intent(inout),dimension(:,:,:) :: &
	468	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: &
	469	Optical_Thickness_vs_ReffIce
	470	real,intent(inout),dimension(:,:,:) :: &
	471	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: &
	472	Optical_Thickness_vs_ReffLiq
	473
	474	! LOCAL VARIABLES
	475	real, parameter :: &
	476	LWP_conversion = 2./3. * 1000. ! MKS units
	477	integer :: i, j
	478	logical, dimension(nPoints,nSubCols) :: &
	479	cloudMask, &
	480	waterCloudMask, &
	481	iceCloudMask, &
	482	validRetrievalMask
	483	real,dimension(nPoints,nSubCols) :: &
	484	tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
	485
	486	! ########################################################################################
	487	! Include only those pixels with successful retrievals in the statistics
	488	! ########################################################################################
	489	validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
	490	cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and. &
	491	validRetrievalMask(1:nPoints,1:nSubCols)
	492	waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. &
	493	validRetrievalMask(1:nPoints,1:nSubCols)
	494	iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .and. &
	495	validRetrievalMask(1:nPoints,1:nSubCols)
	496
	497	! ########################################################################################
	498	! Use these as pixel counts at first
	499	! ########################################################################################
	500	Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2))
	501	Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
	502	Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2))
	503	Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure <= &
	504	highCloudPressureLimit, dim = 2))
	505	Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure > &
	506	lowCloudPressureLimit, dim = 2))
	507	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
	508	- Cloud_Fraction_Low_Mean(1:nPoints)
	509
	510	! ########################################################################################
	511	! Compute mean optical thickness.
	512	! ########################################################################################
[3241]	513	where (Cloud_Fraction_Total_Mean == 0)
	514	Optical_Thickness_Total_Mean = R_UNDEF
	515	elsewhere
	516	Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / &
	517	Cloud_Fraction_Total_Mean(1:nPoints)
	518	endwhere
	519	where (Cloud_Fraction_Water_Mean == 0)
	520	Optical_Thickness_Water_Mean = R_UNDEF
	521	elsewhere
	522	Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
	523	Cloud_Fraction_Water_Mean(1:nPoints)
	524	endwhere
	525	where (Cloud_Fraction_Ice_Mean == 0)
	526	Optical_Thickness_Ice_Mean = R_UNDEF
	527	elsewhere
	528	Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / &
	529	Cloud_Fraction_Ice_Mean(1:nPoints)
	530	endwhere
	531
[2713]	532	! ########################################################################################
	533	! We take the absolute value of optical thickness here to satisfy compilers that complains
	534	! when we evaluate the logarithm of a negative number, even though it's not included in
	535	! the sum.
	536	! ########################################################################################
[3241]	537	where (Cloud_Fraction_Total_Mean == 0)
	538	Optical_Thickness_Total_MeanLog10 = R_UNDEF
	539	elsewhere
	540	Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
	541	dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
	542	endwhere
	543	where (Cloud_Fraction_Water_Mean == 0)
	544	Optical_Thickness_Water_MeanLog10 = R_UNDEF
	545	elsewhere
	546	Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
	547	dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
	548	endwhere
	549	where (Cloud_Fraction_Ice_Mean == 0)
	550	Optical_Thickness_Ice_MeanLog10 = R_UNDEF
	551	elsewhere
	552	Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
	553	dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
	554	endwhere
	555	where (Cloud_Fraction_Water_Mean == 0)
	556	Cloud_Particle_Size_Water_Mean = R_UNDEF
	557	elsewhere
	558	Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
	559	Cloud_Fraction_Water_Mean(1:nPoints)
	560	endwhere
	561	where (Cloud_Fraction_Ice_Mean == 0)
	562	Cloud_Particle_Size_Ice_Mean = R_UNDEF
	563	elsewhere
	564	Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / &
	565	Cloud_Fraction_Ice_Mean(1:nPoints)
	566	endwhere
	567	where (Cloud_Fraction_Total_Mean == 0)
	568	Cloud_Top_Pressure_Total_Mean = R_UNDEF
	569	elsewhere
	570	Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
	571	max(1, count(cloudMask, dim = 2))
	572	endwhere
	573	where (Cloud_Fraction_Water_Mean == 0)
	574	Liquid_Water_Path_Mean = R_UNDEF
	575	elsewhere
	576	Liquid_Water_Path_Mean(1:nPoints) = LWP_conversionsum(particle_sizeoptical_thickness, &
	577	mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
	578	endwhere
	579	where (Cloud_Fraction_Ice_Mean == 0)
	580	Ice_Water_Path_Mean = R_UNDEF
	581	elsewhere
	582	Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_densitysum(particle_sizeoptical_thickness,&
	583	mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
	584	endwhere
[2713]	585	! ########################################################################################
[2822]	586	! Normalize pixel counts to fraction.
[2713]	587	! ########################################################################################
[2822]	588	Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
	589	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols
	590	Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols
	591	Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
	592	Cloud_Fraction_Ice_Mean(1:nPoints) = Cloud_Fraction_Ice_Mean(1:nPoints) /nSubcols
	593	Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
	594
[2713]	595	! ########################################################################################
	596	! Set clear-scenes to undefined
	597	! ########################################################################################
[3286]	598	! where (Cloud_Fraction_High_Mean == 0) Cloud_Fraction_High_Mean = R_UNDEF
	599	! where (Cloud_Fraction_Mid_Mean == 0) Cloud_Fraction_Mid_Mean = R_UNDEF
	600	! where (Cloud_Fraction_Low_Mean == 0) Cloud_Fraction_Low_Mean = R_UNDEF
[2713]	601
	602	! ########################################################################################
	603	! Joint histogram
	604	! ########################################################################################
	605
	606	! Loop over all points
	607	tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols)
	608	ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols)
	609	reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
	610	reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
	611	do j=1,nPoints
	612
	613	! Fill clear and optically thin subcolumns with fill
	614	where(.not. cloudMask(j,1:nSubCols))
	615	tauWRK(j,1:nSubCols) = -999.
	616	ctpWRK(j,1:nSubCols) = -999.
	617	endwhere
	618	! Joint histogram of tau/CTP
	619	call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
	620	tauHistogramBoundaries,numTauHistogramBins,&
	621	pressureHistogramBoundaries,numPressureHistogramBins,&
	622	Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins))
	623	! Joint histogram of tau/ReffICE
	624	call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, &
	625	tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds, &
	626	numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins))
	627	! Joint histogram of tau/ReffLIQ
	628	call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, &
	629	tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds, &
	630	numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins))
	631
	632	enddo
	633	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = &
	634	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols
	635	Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = &
	636	Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols
	637	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = &
	638	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols
	639
	640	end subroutine modis_column
	641	! ######################################################################################
	642	! SUBROUTINE hist2D
	643	! ######################################################################################
	644	subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
	645	implicit none
	646
	647	! INPUTS
	648	integer, intent(in) :: &
	649	npts, & ! Number of data points to be sorted
	650	nbin1, & ! Number of bins in histogram direction 1
	651	nbin2 ! Number of bins in histogram direction 2
	652	real,intent(in),dimension(npts) :: &
	653	var1, & ! Variable 1 to be sorted into bins
	654	var2 ! variable 2 to be sorted into bins
	655	real,intent(in),dimension(nbin1+1) :: &
	656	bin1 ! Histogram bin 1 boundaries
	657	real,intent(in),dimension(nbin2+1) :: &
	658	bin2 ! Histogram bin 2 boundaries
	659	! OUTPUTS
	660	real,intent(out),dimension(nbin1,nbin2) :: &
	661	jointHist
	662
	663	! LOCAL VARIABLES
	664	integer :: ij,ik
	665
	666	do ij=2,nbin1+1
	667	do ik=2,nbin2+1
	668	jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
	669	var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))
	670	enddo
	671	enddo
	672	end subroutine hist2D
	673
[2432]	674	!------------------------------------------------------------------------------------------------
	675	subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, &
	676	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
	677	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
	678	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
	679	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
	680	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
	681	Cloud_Top_Pressure_Total_Mean, &
	682	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
	683	Optical_Thickness_vs_Cloud_Top_Pressure)
	684	!
	685	! Inputs; dimension nPoints, nSubcols
	686	!
	687	integer, dimension(:, :), intent(in) :: phase
	688	real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size
	689	!
	690	! Outputs; dimension nPoints
	691	!
	692	real, dimension(:), intent(out) :: &
	693	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
	694	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
	695	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
	696	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
	697	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
	698	Cloud_Top_Pressure_Total_Mean, &
	699	Liquid_Water_Path_Mean, Ice_Water_Path_Mean
	700	! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins
	701	real, dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure
	702	! ---------------------------
	703	! Local variables
	704	!
	705	real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units
	706	integer :: i, j
	707	integer :: nPoints, nSubcols
	708	logical, dimension(size(phase, 1), size(phase, 2)) :: &
	709	cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask
	710	logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins ) :: tauMask
	711	logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask
	712	! ---------------------------
	713
	714	nPoints = size(phase, 1)
	715	nSubcols = size(phase, 2)
	716	!
	717	! Array conformance checks
	718	!
	719	if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), &
	720	size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), &
	721	size(Cloud_Fraction_High_Mean), size(Cloud_Fraction_Mid_Mean), size(Cloud_Fraction_Low_Mean), &
	722	size(Optical_Thickness_Total_Mean), size(Optical_Thickness_Water_Mean), size(Optical_Thickness_Ice_Mean), &
	723	size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), &
	724	size(Optical_Thickness_Ice_MeanLog10), size(Cloud_Particle_Size_Water_Mean), &
	725	size(Cloud_Particle_Size_Ice_Mean), size(Cloud_Top_Pressure_Total_Mean), &
	726	size(Liquid_Water_Path_Mean), size(Ice_Water_Path_Mean) /) /= nPoints)) &
	727	call complain_and_die("Some L3 arrays have wrong number of grid points")
	728	if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) &
	729	call complain_and_die("Some L3 arrays have wrong number of subcolumns")
	730
	731
	732	!
	733	! Include only those pixels with successful retrievals in the statistics
	734	!
	735	validRetrievalMask(:, :) = particle_size(:, :) > 0.
	736	cloudMask = phase(:, :) /= phaseIsNone .and. validRetrievalMask(:, :)
	737	waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :)
	738	iceCloudMask = phase(:, :) == phaseIsIce .and. validRetrievalMask(:, :)
	739	!
	740	! Use these as pixel counts at first
	741	!
	742	Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2))
	743	Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2))
	744	Cloud_Fraction_Ice_Mean(:) = real(count(iceCloudMask, dim = 2))
	745
	746	Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2))
	747	Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .and. cloud_top_pressure > lowCloudPressureLimit, dim = 2))
	748	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:)
	749
	750	!
	751	! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0.
	752	!
	753	where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1.
	754	where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1.
	755	where (Cloud_Fraction_Ice_Mean == 0) Cloud_Fraction_Ice_Mean = -1.
	756
	757	Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask, dim = 2) / Cloud_Fraction_Total_Mean(:)
	758	Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
	759	Optical_Thickness_Ice_Mean = sum(optical_thickness, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
	760
	761	! We take the absolute value of optical thickness here to satisfy compilers that complains when we
	762	! evaluate the logarithm of a negative number, even though it's not included in the sum.
	763	Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask, dim = 2) / &
	764	Cloud_Fraction_Total_Mean(:)
	765	Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / &
	766	Cloud_Fraction_Water_Mean(:)
	767	Optical_Thickness_Ice_MeanLog10 = sum(log10(abs(optical_thickness)), mask = iceCloudMask, dim = 2) / &
	768	Cloud_Fraction_Ice_Mean(:)
	769
	770	Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
	771	Cloud_Particle_Size_Ice_Mean = sum(particle_size, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
	772
	773	Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2))
	774
	775	Liquid_Water_Path_Mean = LWP_conversion &
	776	* sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) &
	777	/ Cloud_Fraction_Water_Mean(:)
	778	Ice_Water_Path_Mean = LWP_conversion * ice_density &
	779	* sum(particle_size * optical_thickness, mask = iceCloudMask, dim = 2) &
	780	/ Cloud_Fraction_Ice_Mean(:)
	781
	782	!
	783	! Normalize pixel counts to fraction
	784	! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0.
	785	!
	786	Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols)
	787	Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols)
	788	Cloud_Fraction_Ice_Mean(:) = max(0., Cloud_Fraction_Ice_Mean(:) /nSubcols)
	789
	790	Cloud_Fraction_High_Mean(:) = Cloud_Fraction_High_Mean(:) /nSubcols
	791	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Mid_Mean(:) /nSubcols
	792	Cloud_Fraction_Low_Mean(:) = Cloud_Fraction_Low_Mean(:) /nSubcols
	793
	794	! ----
	795	! Joint histogram
	796	!
	797	do i = 1, numTauHistogramBins
	798	where(cloudMask(:, :))
	799	tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. &
	800	optical_thickness(:, :) < tauHistogramBoundaries(i+1)
	801	elsewhere
	802	tauMask(:, :, i) = .false.
	803	end where
	804	end do
	805
	806	do i = 1, numPressureHistogramBins
	807	where(cloudMask(:, :))
	808	pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. &
	809	cloud_top_pressure(:, :) < pressureHistogramBoundaries(i+1)
	810	elsewhere
	811	pressureMask(:, :, i) = .false.
	812	end where
	813	end do
	814
	815	do i = 1, numPressureHistogramBins
	816	do j = 1, numTauHistogramBins
	817	Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = &
	818	real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
	819	end do
	820	end do
	821
	822	end subroutine modis_L3_simulator
	823	!------------------------------------------------------------------------------------------------
	824	function cloud_top_pressure(tauIncrement, pressure, tauLimit)
	825	real, dimension(:), intent(in) :: tauIncrement, pressure
	826	real, intent(in) :: tauLimit
	827	real :: cloud_top_pressure
	828	!
	829	! Find the extinction-weighted pressure. Assume that pressure varies linearly between
	830	! layers and use the trapezoidal rule.
	831	!
	832
	833	real :: deltaX, totalTau, totalProduct
	834	integer :: i
	835
	836	totalTau = 0.; totalProduct = 0.
	837	do i = 2, size(tauIncrement)
	838	if(totalTau + tauIncrement(i) > tauLimit) then
	839	deltaX = tauLimit - totalTau
	840	totalTau = totalTau + deltaX
	841	!
	842	! Result for trapezoidal rule when you take less than a full step
	843	! tauIncrement is a layer-integrated value
	844	!
	845	totalProduct = totalProduct &
	846	+ pressure(i-1) * deltaX &
	847	+ (pressure(i) - pressure(i-1)) * deltaX*2/(2. tauIncrement(i))
	848	else
	849	totalTau = totalTau + tauIncrement(i)
	850	totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2.
	851	end if
	852	if(totalTau >= tauLimit) exit
	853	end do
	854	cloud_top_pressure = totalProduct/totalTau
	855	end function cloud_top_pressure
	856	!------------------------------------------------------------------------------------------------
	857	function weight_by_extinction(tauIncrement, f, tauLimit)
	858	real, dimension(:), intent(in) :: tauIncrement, f
	859	real, intent(in) :: tauLimit
	860	real :: weight_by_extinction
	861	!
	862	! Find the extinction-weighted value of f(tau), assuming constant f within each layer
	863	!
	864
	865	real :: deltaX, totalTau, totalProduct
	866	integer :: i
	867
	868	totalTau = 0.; totalProduct = 0.
	869	do i = 1, size(tauIncrement)
	870	if(totalTau + tauIncrement(i) > tauLimit) then
	871	deltaX = tauLimit - totalTau
	872	totalTau = totalTau + deltaX
	873	totalProduct = totalProduct + deltaX * f(i)
	874	else
	875	totalTau = totalTau + tauIncrement(i)
	876	totalProduct = totalProduct + tauIncrement(i) * f(i)
	877	end if
	878	if(totalTau >= tauLimit) exit
	879	end do
	880	weight_by_extinction = totalProduct/totalTau
	881	end function weight_by_extinction
	882	!------------------------------------------------------------------------------------------------
	883	pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size)
	884	real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
	885	real :: compute_nir_reflectance
	886
	887	real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
	888	tau, g, w0
	889	!----------------------------------------
	890	water_g(:) = get_g_nir( phaseIsLiquid, water_size)
	891	water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size)
	892	ice_g(:) = get_g_nir( phaseIsIce, ice_size)
	893	ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size)
	894	!
	895	! Combine ice and water optical properties
	896	!
	897	g(:) = 0; w0(:) = 0.
	898	tau(:) = ice_tau(:) + water_tau(:)
	899	where (tau(:) > 0)
[3184]	900	w0(:) = (water_tau(:) * water_w0(:) + ice_tau(:) * ice_w0(:)) / &
	901	tau(:)
	902	g(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
	903	(w0(:) * tau(:))
[2432]	904	end where
	905
	906	compute_nir_reflectance = compute_toa_reflectace(tau, g, w0)
	907	end function compute_nir_reflectance
	908	!------------------------------------------------------------------------------------------------
	909	! Retreivals
	910	!------------------------------------------------------------------------------------------------
	911	elemental function retrieve_re (phase, tau, obs_Refl_nir)
	912	integer, intent(in) :: phase
	913	real, intent(in) :: tau, obs_Refl_nir
	914	real :: retrieve_re
	915	!
	916	! Finds the re that produces the minimum mis-match between predicted and observed reflectance in
	917	! MODIS band 7 (near IR)
	918	! Uses
	919	! fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables
	920	! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0
	921	! adding-doubling for total reflectance
	922	!
	923	!
	924	!
	925	! Local variables
	926	!
	927	real, parameter :: min_distance_to_boundary = 0.01
	928	real :: re_min, re_max, delta_re
	929	integer :: i
	930
	931	real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir
	932	! --------------------------
	933
	934	if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then
	935	if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then
	936	re_min = re_water_min
	937	re_max = re_water_max
	938	trial_re(:) = trial_re_w
	939	g(:) = g_w(:)
	940	w0(:) = w0_w(:)
	941	else
	942	re_min = re_ice_min
	943	re_max = re_ice_max
	944	trial_re(:) = trial_re_i
	945	g(:) = g_i(:)
	946	w0(:) = w0_i(:)
	947	end if
	948	!
	949	! 1st attempt at index: w/coarse re resolution
	950	!
	951	predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
	952	retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
	953	!
	954	! If first retrieval works, can try 2nd iteration using greater re resolution
	955	!
[2713]	956	! DJS2015: Remove unused piece of code
	957	! if(use_two_re_iterations .and. retrieve_re > 0.) then
	958	! re_min = retrieve_re - delta_re
	959	! re_max = retrieve_re + delta_re
	960	! delta_re = (re_max - re_min)/real(num_trial_res-1)
	961	!
	962	! trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
	963	! g(:) = get_g_nir( phase, trial_re(:))
	964	! w0(:) = get_ssa_nir(phase, trial_re(:))
	965	! predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
	966	! retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
	967	! end if
	968	! DJS2015 END
[2432]	969	else
	970	retrieve_re = re_fill
	971	end if
	972
	973	end function retrieve_re
	974	! --------------------------------------------
	975	pure function interpolate_to_min(x, y, yobs)
	976	real, dimension(:), intent(in) :: x, y
	977	real, intent(in) :: yobs
	978	real :: interpolate_to_min
	979	!
	980	! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
	981	! y must be monotonic in x
	982	!
	983	real, dimension(size(x)) :: diff
	984	integer :: nPoints, minDiffLoc, lowerBound, upperBound
	985	! ---------------------------------
	986	nPoints = size(y)
	987	diff(:) = y(:) - yobs
	988	minDiffLoc = minloc(abs(diff), dim = 1)
	989
	990	if(minDiffLoc == 1) then
	991	lowerBound = minDiffLoc
	992	upperBound = minDiffLoc + 1
	993	else if(minDiffLoc == nPoints) then
	994	lowerBound = minDiffLoc - 1
	995	upperBound = minDiffLoc
	996	else
	997	if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
	998	lowerBound = minDiffLoc-1
	999	upperBound = minDiffLoc
	1000	else
	1001	lowerBound = minDiffLoc
	1002	upperBound = minDiffLoc + 1
	1003	end if
	1004	end if
	1005
	1006	if(diff(lowerBound) * diff(upperBound) < 0) then
	1007	!
	1008	! Interpolate the root position linearly if we bracket the root
	1009	!
	1010	interpolate_to_min = x(upperBound) - &
	1011	diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
	1012	else
	1013	interpolate_to_min = re_fill
	1014	end if
	1015
	1016
	1017	end function interpolate_to_min
	1018	! --------------------------------------------
	1019	! Optical properties
	1020	! --------------------------------------------
	1021	elemental function get_g_nir (phase, re)
	1022	!
	1023	! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
	1024	! of size for ice and water
	1025	! Fits from Steve Platnick
	1026	!
	1027
	1028	integer, intent(in) :: phase
	1029	real, intent(in) :: re
	1030	real :: get_g_nir
[2713]	1031
	1032	real, dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
	1033	small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
	1034	real, dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
	1035
	1036	! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
	1037	if(phase == phaseIsLiquid) then
	1038	if(re < 7.) then
	1039	get_g_nir = fit_to_quadratic(re, small_water_coefficients)
	1040	if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
	1041	else
	1042	get_g_nir = fit_to_cubic(re, big_water_coefficients)
	1043	if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
	1044	end if
[2432]	1045	else
[2713]	1046	get_g_nir = fit_to_quadratic(re, ice_coefficients)
[2432]	1047	if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
	1048	if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
	1049	end if
	1050
	1051	end function get_g_nir
	1052
	1053	! --------------------------------------------
	1054	elemental function get_ssa_nir (phase, re)
	1055	integer, intent(in) :: phase
	1056	real, intent(in) :: re
	1057	real :: get_ssa_nir
	1058	!
	1059	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
	1060	! of size for ice and water
	1061	! Fits from Steve Platnick
	1062	!
[2713]	1063	real, dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
	1064	real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
[2432]	1065
[2713]	1066	! approx. fits from MODIS Collection 6 LUT scattering calculations
[2432]	1067	if(phase == phaseIsLiquid) then
	1068	get_ssa_nir = fit_to_quadratic(re, water_coefficients)
	1069	if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
	1070	if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
	1071	else
	1072	get_ssa_nir = fit_to_cubic(re, ice_coefficients)
	1073	if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
	1074	if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
	1075	end if
	1076
	1077	end function get_ssa_nir
	1078	! --------------------------------------------
	1079	pure function fit_to_cubic(x, coefficients)
	1080	real, intent(in) :: x
	1081	real, dimension(:), intent(in) :: coefficients
	1082	real :: fit_to_cubic
	1083
	1084
	1085	fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
	1086	end function fit_to_cubic
	1087	! --------------------------------------------
	1088	pure function fit_to_quadratic(x, coefficients)
	1089	real, intent(in) :: x
	1090	real, dimension(:), intent(in) :: coefficients
	1091	real :: fit_to_quadratic
	1092
	1093
	1094	fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
	1095	end function fit_to_quadratic
	1096	! --------------------------------------------
	1097	! Radiative transfer
	1098	! --------------------------------------------
	1099	pure function compute_toa_reflectace(tau, g, w0)
	1100	real, dimension(:), intent(in) :: tau, g, w0
	1101	real :: compute_toa_reflectace
	1102
	1103	logical, dimension(size(tau)) :: cloudMask
	1104	integer, dimension(count(tau(:) > 0)) :: cloudIndicies
	1105	real, dimension(count(tau(:) > 0)) :: Refl, Trans
	1106	real :: Refl_tot, Trans_tot
	1107	integer :: i
	1108	! ---------------------------------------
	1109	!
	1110	! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation
	1111	!
	1112	cloudMask = tau(:) > 0.
	1113	cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask)
	1114	do i = 1, size(cloudIndicies)
	1115	call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
	1116	end do
	1117
	1118	call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot)
	1119
	1120	compute_toa_reflectace = Refl_tot
	1121
	1122	end function compute_toa_reflectace
	1123	! --------------------------------------------
	1124	pure subroutine two_stream(tauint, gint, w0int, ref, tra)
	1125	real, intent(in) :: tauint, gint, w0int
	1126	real, intent(out) :: ref, tra
	1127	!
	1128	! Compute reflectance in a single layer using the two stream approximation
	1129	! The code itself is from Lazaros Oreopoulos via Steve Platnick
	1130	!
	1131	! ------------------------
	1132	! Local variables
	1133	! for delta Eddington code
	1134	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
	1135	integer, parameter :: beam = 2
	1136	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
	1137	real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
	1138	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
	1139	!
	1140	! Compute reflectance and transmittance in a single layer using the two stream approximation
	1141	! The code itself is from Lazaros Oreopoulos via Steve Platnick
	1142	!
	1143	f = gint**2
	1144	tau = (1 - w0int * f) * tauint
	1145	w0 = (1 - f) * w0int / (1 - w0int * f)
	1146	g = (gint - f) / (1 - f)
	1147
	1148	! delta-Eddington (Joseph et al. 1976)
	1149	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
	1150	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
	1151	gamma3 = (2 - 3gxmu) / 4.0
	1152	gamma4 = 1 - gamma3
	1153
	1154	if (w0int > minConservativeW0) then
	1155	! Conservative scattering
	1156	if (beam == 1) then
	1157	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
	1158	ref = rh / (1 + gamma1 * tau)
	1159	tra = 1 - ref
	1160	else if(beam == 2) then
	1161	ref = gamma1tau/(1 + gamma1tau)
	1162	tra = 1 - ref
	1163	endif
	1164	else
	1165	! Non-conservative scattering
	1166	a1 = gamma1 * gamma4 + gamma2 * gamma3
	1167	a2 = gamma1 * gamma3 + gamma2 * gamma4
	1168
	1169	rk = sqrt(gamma12 - gamma22)
	1170
	1171	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
	1172	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
	1173	r3 = 2 * rk (gamma3 - a2 xmu)
	1174	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
	1175	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
	1176
	1177	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
	1178	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
	1179	t3 = 2 * rk * (gamma4 + a1 * xmu)
	1180	t4 = r4
	1181	t5 = r5
	1182
	1183	beta = -r5 / r4
	1184
	1185	e1 = min(rk * tau, 500.)
	1186	e2 = min(tau / xmu, 500.)
	1187
	1188	if (beam == 1) then
	1189	den = r4 * exp(e1) + r5 * exp(-e1)
	1190	ref = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
	1191	den = t4 * exp(e1) + t5 * exp(-e1)
	1192	th = exp(-e2)
	1193	tra = th-thw0(t1exp(e1)-t2exp(-e1)-t3*exp(e2))/den
	1194	elseif (beam == 2) then
	1195	ef1 = exp(-e1)
	1196	ef2 = exp(-2*e1)
	1197	ref = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
	1198	tra = (2rkef1)/((rk+gamma1)(1-betaef2))
	1199	endif
	1200	end if
	1201	end subroutine two_stream
	1202	! --------------------------------------------------
	1203	elemental function two_stream_reflectance(tauint, gint, w0int)
	1204	real, intent(in) :: tauint, gint, w0int
	1205	real :: two_stream_reflectance
	1206	!
	1207	! Compute reflectance in a single layer using the two stream approximation
	1208	! The code itself is from Lazaros Oreopoulos via Steve Platnick
	1209	!
	1210	! ------------------------
	1211	! Local variables
	1212	! for delta Eddington code
	1213	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
	1214	integer, parameter :: beam = 2
	1215	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
	1216	real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
	1217	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
	1218	! ------------------------
	1219
	1220
	1221	f = gint**2
	1222	tau = (1 - w0int * f) * tauint
	1223	w0 = (1 - f) * w0int / (1 - w0int * f)
	1224	g = (gint - f) / (1 - f)
	1225
	1226	! delta-Eddington (Joseph et al. 1976)
	1227	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
	1228	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
	1229	gamma3 = (2 - 3gxmu) / 4.0
	1230	gamma4 = 1 - gamma3
	1231
	1232	if (w0int > minConservativeW0) then
	1233	! Conservative scattering
	1234	if (beam == 1) then
	1235	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
	1236	two_stream_reflectance = rh / (1 + gamma1 * tau)
	1237	elseif (beam == 2) then
	1238	two_stream_reflectance = gamma1tau/(1 + gamma1tau)
	1239	endif
	1240
	1241	else !
	1242
	1243	! Non-conservative scattering
	1244	a1 = gamma1 * gamma4 + gamma2 * gamma3
	1245	a2 = gamma1 * gamma3 + gamma2 * gamma4
	1246
	1247	rk = sqrt(gamma12 - gamma22)
	1248
	1249	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
	1250	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
	1251	r3 = 2 * rk (gamma3 - a2 xmu)
	1252	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
	1253	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
	1254
	1255	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
	1256	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
	1257	t3 = 2 * rk * (gamma4 + a1 * xmu)
	1258	t4 = r4
	1259	t5 = r5
	1260
	1261	beta = -r5 / r4
	1262
	1263	e1 = min(rk * tau, 500.)
	1264	e2 = min(tau / xmu, 500.)
	1265
	1266	if (beam == 1) then
	1267	den = r4 * exp(e1) + r5 * exp(-e1)
	1268	two_stream_reflectance = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
	1269	elseif (beam == 2) then
	1270	ef1 = exp(-e1)
	1271	ef2 = exp(-2*e1)
	1272	two_stream_reflectance = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
	1273	endif
	1274
	1275	end if
	1276	end function two_stream_reflectance
	1277	! --------------------------------------------
	1278	pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot)
	1279	real, dimension(:), intent(in) :: Refl, Tran
	1280	real, intent(out) :: Refl_tot, Tran_tot
	1281	!
	1282	! Use adding/doubling formulas to compute total reflectance and transmittance from layer values
	1283	!
	1284
	1285	integer :: i
	1286	real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative
	1287
	1288	Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)
	1289
	1290	do i=2, size(Refl)
	1291	! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
	1292	Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)(Tran_cumulative(i-1)2)/(1 - Refl_cumulative(i-1) Refl(i))
	1293	Tran_cumulative(i) = (Tran_cumulative(i-1)Tran(i)) / (1 - Refl_cumulative(i-1) Refl(i))
	1294	end do
	1295
	1296	Refl_tot = Refl_cumulative(size(Refl))
	1297	Tran_tot = Tran_cumulative(size(Refl))
	1298
	1299	end subroutine adding_doubling
	1300	! --------------------------------------------------
	1301	subroutine complain_and_die(message)
	1302	character(len = *), intent(in) :: message
	1303
	1304	write(6, *) "Failure in MODIS simulator"
	1305	write(6, *) trim(message)
	1306	stop
	1307	end subroutine complain_and_die
	1308	!------------------------------------------------------------------------------------------------
[2713]	1309	end module mod_modis_sim

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format