Context Navigation

modis_simulator.F90 @ 3603

Last change on this file since 3603 was 2822, checked in by idelkadi, 8 years ago

Ecriture des champs de sorties Cosp 4D dans les fichiers Netcdf en utilisant les fonctionalites Xios.

ajout d'une routine histwrite_cosp4D appelee si les sorties sont gerees par Xios (dans cosp_output_write.F90)
definitions des differents axes et grilles 4D pour Cosp (dans context_lmdz.xml)

Traitement des valeurs indefinies par Xios et mise a zero si l'on utilise Ioipsl.

Correction d'erreur sur les sorties Modis

File size: 63.3 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
	54	USE MOD_COSP_TYPES, only: R_UNDEF
	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
	386	! ---------------------------------------------------
	387
	388	where(cloudIce(:, :) <= 0.)
	389	tauLiquidFraction(:, :) = 1.
	390	elsewhere
	391	where (cloudWater(:, :) <= 0.)
	392	tauLiquidFraction(:, :) = 0.
	393	elsewhere
	394	!
	395	! Geometic optics limit - tau as LWP/re (proportional to LWC/re)
	396	!
	397	tauLiquidFraction(:, :) = (cloudWater(:, :)/waterSize(:, :)) / &
	398	(cloudWater(:, :)/waterSize(:, :) + cloudIce(:, :)/(ice_density * iceSize(:, :)) )
	399	end where
	400	end where
	401	liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :)
	402	ice_opticalThickness (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :)
	403
	404	call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels, &
	405	liquid_opticalThickness, ice_opticalThickness, &
	406	waterSize, iceSize, &
	407	isccpTau, isccpCloudTopPressure, &
	408	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
	409
	410	end subroutine modis_L2_simulator_oneTau
[2713]	411
	412	! ########################################################################################
	413	subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, &
	414	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
	415	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
	416	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
	417	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
	418	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, &
	419	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
	420	Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
	421
	422	! INPUTS
	423	integer,intent(in) :: &
	424	nPoints, & ! Number of horizontal gridpoints
	425	nSubCols ! Number of subcolumns
	426	integer,intent(in), dimension(:,:) :: &
	427	!ds integer,intent(in), dimension(nPoints, nSubCols) :: &
	428	phase
	429	real,intent(in),dimension(:,:) :: &
	430	!ds real,intent(in),dimension(nPoints, nSubCols) :: &
	431	cloud_top_pressure, &
	432	optical_thickness, &
	433	particle_size
	434
	435	! OUTPUTS
	436	real,intent(inout),dimension(:) :: & !
	437	!ds real,intent(inout),dimension(nPoints) :: & !
	438	Cloud_Fraction_Total_Mean, & !
	439	Cloud_Fraction_Water_Mean, & !
	440	Cloud_Fraction_Ice_Mean, & !
	441	Cloud_Fraction_High_Mean, & !
	442	Cloud_Fraction_Mid_Mean, & !
	443	Cloud_Fraction_Low_Mean, & !
	444	Optical_Thickness_Total_Mean, & !
	445	Optical_Thickness_Water_Mean, & !
	446	Optical_Thickness_Ice_Mean, & !
	447	Optical_Thickness_Total_MeanLog10, & !
	448	Optical_Thickness_Water_MeanLog10, & !
	449	Optical_Thickness_Ice_MeanLog10, & !
	450	Cloud_Particle_Size_Water_Mean, & !
	451	Cloud_Particle_Size_Ice_Mean, & !
	452	Cloud_Top_Pressure_Total_Mean, & !
	453	Liquid_Water_Path_Mean, & !
	454	Ice_Water_Path_Mean !
	455	real,intent(inout),dimension(:,:,:) :: &
	456	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: &
	457	Optical_Thickness_vs_Cloud_Top_Pressure
	458	real,intent(inout),dimension(:,:,:) :: &
	459	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: &
	460	Optical_Thickness_vs_ReffIce
	461	real,intent(inout),dimension(:,:,:) :: &
	462	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: &
	463	Optical_Thickness_vs_ReffLiq
	464
	465	! LOCAL VARIABLES
	466	real, parameter :: &
	467	LWP_conversion = 2./3. * 1000. ! MKS units
	468	integer :: i, j
	469	logical, dimension(nPoints,nSubCols) :: &
	470	cloudMask, &
	471	waterCloudMask, &
	472	iceCloudMask, &
	473	validRetrievalMask
	474	real,dimension(nPoints,nSubCols) :: &
	475	tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
	476
	477	! ########################################################################################
	478	! Include only those pixels with successful retrievals in the statistics
	479	! ########################################################################################
	480	validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
	481	cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and. &
	482	validRetrievalMask(1:nPoints,1:nSubCols)
	483	waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. &
	484	validRetrievalMask(1:nPoints,1:nSubCols)
	485	iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .and. &
	486	validRetrievalMask(1:nPoints,1:nSubCols)
	487
	488	! ########################################################################################
	489	! Use these as pixel counts at first
	490	! ########################################################################################
	491	Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2))
	492	Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
	493	Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2))
	494	Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure <= &
	495	highCloudPressureLimit, dim = 2))
	496	Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure > &
	497	lowCloudPressureLimit, dim = 2))
	498	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
	499	- Cloud_Fraction_Low_Mean(1:nPoints)
	500
	501	! ########################################################################################
	502	! Compute mean optical thickness.
	503	! ########################################################################################
	504	Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / &
	505	Cloud_Fraction_Total_Mean(1:nPoints)
	506	Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
	507	Cloud_Fraction_Water_Mean(1:nPoints)
	508	Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / &
	509	Cloud_Fraction_Ice_Mean(1:nPoints)
	510
	511	! ########################################################################################
	512	! We take the absolute value of optical thickness here to satisfy compilers that complains
	513	! when we evaluate the logarithm of a negative number, even though it's not included in
	514	! the sum.
	515	! ########################################################################################
	516	Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
	517	dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
	518	Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
	519	dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
	520	Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
	521	dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
	522	Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
	523	Cloud_Fraction_Water_Mean(1:nPoints)
	524	Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / &
	525	Cloud_Fraction_Ice_Mean(1:nPoints)
	526	Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
	527	max(1, count(cloudMask, dim = 2))
	528	Liquid_Water_Path_Mean(1:nPoints) = LWP_conversionsum(particle_sizeoptical_thickness, &
	529	mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
	530	Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_densitysum(particle_sizeoptical_thickness,&
	531	mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
	532
	533	! ########################################################################################
[2822]	534	! Normalize pixel counts to fraction.
[2713]	535	! ########################################################################################
[2822]	536	Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
	537	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols
	538	Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols
	539	Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
	540	Cloud_Fraction_Ice_Mean(1:nPoints) = Cloud_Fraction_Ice_Mean(1:nPoints) /nSubcols
	541	Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
	542
[2713]	543	! ########################################################################################
	544	! Set clear-scenes to undefined
	545	! ########################################################################################
	546	where (Cloud_Fraction_Total_Mean == 0)
	547	Optical_Thickness_Total_Mean = R_UNDEF
	548	Optical_Thickness_Total_MeanLog10 = R_UNDEF
	549	Cloud_Top_Pressure_Total_Mean = R_UNDEF
	550	endwhere
	551	where (Cloud_Fraction_Water_Mean == 0)
	552	Optical_Thickness_Water_Mean = R_UNDEF
	553	Optical_Thickness_Water_MeanLog10 = R_UNDEF
	554	Cloud_Particle_Size_Water_Mean = R_UNDEF
	555	Liquid_Water_Path_Mean = R_UNDEF
	556	endwhere
	557	where (Cloud_Fraction_Ice_Mean == 0)
	558	Optical_Thickness_Ice_Mean = R_UNDEF
	559	Optical_Thickness_Ice_MeanLog10 = R_UNDEF
	560	Cloud_Particle_Size_Ice_Mean = R_UNDEF
	561	Ice_Water_Path_Mean = R_UNDEF
	562	endwhere
	563	where (Cloud_Fraction_High_Mean == 0) Cloud_Fraction_High_Mean = R_UNDEF
	564	where (Cloud_Fraction_Mid_Mean == 0) Cloud_Fraction_Mid_Mean = R_UNDEF
	565	where (Cloud_Fraction_Low_Mean == 0) Cloud_Fraction_Low_Mean = R_UNDEF
	566
	567	! ########################################################################################
	568	! Joint histogram
	569	! ########################################################################################
	570
	571	! Loop over all points
	572	tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols)
	573	ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols)
	574	reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
	575	reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
	576	do j=1,nPoints
	577
	578	! Fill clear and optically thin subcolumns with fill
	579	where(.not. cloudMask(j,1:nSubCols))
	580	tauWRK(j,1:nSubCols) = -999.
	581	ctpWRK(j,1:nSubCols) = -999.
	582	endwhere
	583	! Joint histogram of tau/CTP
	584	call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
	585	tauHistogramBoundaries,numTauHistogramBins,&
	586	pressureHistogramBoundaries,numPressureHistogramBins,&
	587	Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins))
	588	! Joint histogram of tau/ReffICE
	589	call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, &
	590	tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds, &
	591	numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins))
	592	! Joint histogram of tau/ReffLIQ
	593	call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, &
	594	tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds, &
	595	numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins))
	596
	597	enddo
	598	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = &
	599	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols
	600	Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = &
	601	Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols
	602	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = &
	603	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols
	604
	605	end subroutine modis_column
	606	! ######################################################################################
	607	! SUBROUTINE hist2D
	608	! ######################################################################################
	609	subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
	610	implicit none
	611
	612	! INPUTS
	613	integer, intent(in) :: &
	614	npts, & ! Number of data points to be sorted
	615	nbin1, & ! Number of bins in histogram direction 1
	616	nbin2 ! Number of bins in histogram direction 2
	617	real,intent(in),dimension(npts) :: &
	618	var1, & ! Variable 1 to be sorted into bins
	619	var2 ! variable 2 to be sorted into bins
	620	real,intent(in),dimension(nbin1+1) :: &
	621	bin1 ! Histogram bin 1 boundaries
	622	real,intent(in),dimension(nbin2+1) :: &
	623	bin2 ! Histogram bin 2 boundaries
	624	! OUTPUTS
	625	real,intent(out),dimension(nbin1,nbin2) :: &
	626	jointHist
	627
	628	! LOCAL VARIABLES
	629	integer :: ij,ik
	630
	631	do ij=2,nbin1+1
	632	do ik=2,nbin2+1
	633	jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
	634	var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))
	635	enddo
	636	enddo
	637	end subroutine hist2D
	638
[2432]	639	!------------------------------------------------------------------------------------------------
	640	subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, &
	641	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
	642	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
	643	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
	644	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
	645	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
	646	Cloud_Top_Pressure_Total_Mean, &
	647	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
	648	Optical_Thickness_vs_Cloud_Top_Pressure)
	649	!
	650	! Inputs; dimension nPoints, nSubcols
	651	!
	652	integer, dimension(:, :), intent(in) :: phase
	653	real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size
	654	!
	655	! Outputs; dimension nPoints
	656	!
	657	real, dimension(:), intent(out) :: &
	658	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
	659	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
	660	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
	661	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
	662	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
	663	Cloud_Top_Pressure_Total_Mean, &
	664	Liquid_Water_Path_Mean, Ice_Water_Path_Mean
	665	! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins
	666	real, dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure
	667	! ---------------------------
	668	! Local variables
	669	!
	670	real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units
	671	integer :: i, j
	672	integer :: nPoints, nSubcols
	673	logical, dimension(size(phase, 1), size(phase, 2)) :: &
	674	cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask
	675	logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins ) :: tauMask
	676	logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask
	677	! ---------------------------
	678
	679	nPoints = size(phase, 1)
	680	nSubcols = size(phase, 2)
	681	!
	682	! Array conformance checks
	683	!
	684	if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), &
	685	size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), &
	686	size(Cloud_Fraction_High_Mean), size(Cloud_Fraction_Mid_Mean), size(Cloud_Fraction_Low_Mean), &
	687	size(Optical_Thickness_Total_Mean), size(Optical_Thickness_Water_Mean), size(Optical_Thickness_Ice_Mean), &
	688	size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), &
	689	size(Optical_Thickness_Ice_MeanLog10), size(Cloud_Particle_Size_Water_Mean), &
	690	size(Cloud_Particle_Size_Ice_Mean), size(Cloud_Top_Pressure_Total_Mean), &
	691	size(Liquid_Water_Path_Mean), size(Ice_Water_Path_Mean) /) /= nPoints)) &
	692	call complain_and_die("Some L3 arrays have wrong number of grid points")
	693	if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) &
	694	call complain_and_die("Some L3 arrays have wrong number of subcolumns")
	695
	696
	697	!
	698	! Include only those pixels with successful retrievals in the statistics
	699	!
	700	validRetrievalMask(:, :) = particle_size(:, :) > 0.
	701	cloudMask = phase(:, :) /= phaseIsNone .and. validRetrievalMask(:, :)
	702	waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :)
	703	iceCloudMask = phase(:, :) == phaseIsIce .and. validRetrievalMask(:, :)
	704	!
	705	! Use these as pixel counts at first
	706	!
	707	Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2))
	708	Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2))
	709	Cloud_Fraction_Ice_Mean(:) = real(count(iceCloudMask, dim = 2))
	710
	711	Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2))
	712	Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .and. cloud_top_pressure > lowCloudPressureLimit, dim = 2))
	713	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:)
	714
	715	!
	716	! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0.
	717	!
	718	where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1.
	719	where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1.
	720	where (Cloud_Fraction_Ice_Mean == 0) Cloud_Fraction_Ice_Mean = -1.
	721
	722	Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask, dim = 2) / Cloud_Fraction_Total_Mean(:)
	723	Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
	724	Optical_Thickness_Ice_Mean = sum(optical_thickness, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
	725
	726	! We take the absolute value of optical thickness here to satisfy compilers that complains when we
	727	! evaluate the logarithm of a negative number, even though it's not included in the sum.
	728	Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask, dim = 2) / &
	729	Cloud_Fraction_Total_Mean(:)
	730	Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / &
	731	Cloud_Fraction_Water_Mean(:)
	732	Optical_Thickness_Ice_MeanLog10 = sum(log10(abs(optical_thickness)), mask = iceCloudMask, dim = 2) / &
	733	Cloud_Fraction_Ice_Mean(:)
	734
	735	Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
	736	Cloud_Particle_Size_Ice_Mean = sum(particle_size, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
	737
	738	Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2))
	739
	740	Liquid_Water_Path_Mean = LWP_conversion &
	741	* sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) &
	742	/ Cloud_Fraction_Water_Mean(:)
	743	Ice_Water_Path_Mean = LWP_conversion * ice_density &
	744	* sum(particle_size * optical_thickness, mask = iceCloudMask, dim = 2) &
	745	/ Cloud_Fraction_Ice_Mean(:)
	746
	747	!
	748	! Normalize pixel counts to fraction
	749	! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0.
	750	!
	751	Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols)
	752	Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols)
	753	Cloud_Fraction_Ice_Mean(:) = max(0., Cloud_Fraction_Ice_Mean(:) /nSubcols)
	754
	755	Cloud_Fraction_High_Mean(:) = Cloud_Fraction_High_Mean(:) /nSubcols
	756	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Mid_Mean(:) /nSubcols
	757	Cloud_Fraction_Low_Mean(:) = Cloud_Fraction_Low_Mean(:) /nSubcols
	758
	759	! ----
	760	! Joint histogram
	761	!
	762	do i = 1, numTauHistogramBins
	763	where(cloudMask(:, :))
	764	tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. &
	765	optical_thickness(:, :) < tauHistogramBoundaries(i+1)
	766	elsewhere
	767	tauMask(:, :, i) = .false.
	768	end where
	769	end do
	770
	771	do i = 1, numPressureHistogramBins
	772	where(cloudMask(:, :))
	773	pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. &
	774	cloud_top_pressure(:, :) < pressureHistogramBoundaries(i+1)
	775	elsewhere
	776	pressureMask(:, :, i) = .false.
	777	end where
	778	end do
	779
	780	do i = 1, numPressureHistogramBins
	781	do j = 1, numTauHistogramBins
	782	Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = &
	783	real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
	784	end do
	785	end do
	786
	787	end subroutine modis_L3_simulator
	788	!------------------------------------------------------------------------------------------------
	789	function cloud_top_pressure(tauIncrement, pressure, tauLimit)
	790	real, dimension(:), intent(in) :: tauIncrement, pressure
	791	real, intent(in) :: tauLimit
	792	real :: cloud_top_pressure
	793	!
	794	! Find the extinction-weighted pressure. Assume that pressure varies linearly between
	795	! layers and use the trapezoidal rule.
	796	!
	797
	798	real :: deltaX, totalTau, totalProduct
	799	integer :: i
	800
	801	totalTau = 0.; totalProduct = 0.
	802	do i = 2, size(tauIncrement)
	803	if(totalTau + tauIncrement(i) > tauLimit) then
	804	deltaX = tauLimit - totalTau
	805	totalTau = totalTau + deltaX
	806	!
	807	! Result for trapezoidal rule when you take less than a full step
	808	! tauIncrement is a layer-integrated value
	809	!
	810	totalProduct = totalProduct &
	811	+ pressure(i-1) * deltaX &
	812	+ (pressure(i) - pressure(i-1)) * deltaX*2/(2. tauIncrement(i))
	813	else
	814	totalTau = totalTau + tauIncrement(i)
	815	totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2.
	816	end if
	817	if(totalTau >= tauLimit) exit
	818	end do
	819	cloud_top_pressure = totalProduct/totalTau
	820	end function cloud_top_pressure
	821	!------------------------------------------------------------------------------------------------
	822	function weight_by_extinction(tauIncrement, f, tauLimit)
	823	real, dimension(:), intent(in) :: tauIncrement, f
	824	real, intent(in) :: tauLimit
	825	real :: weight_by_extinction
	826	!
	827	! Find the extinction-weighted value of f(tau), assuming constant f within each layer
	828	!
	829
	830	real :: deltaX, totalTau, totalProduct
	831	integer :: i
	832
	833	totalTau = 0.; totalProduct = 0.
	834	do i = 1, size(tauIncrement)
	835	if(totalTau + tauIncrement(i) > tauLimit) then
	836	deltaX = tauLimit - totalTau
	837	totalTau = totalTau + deltaX
	838	totalProduct = totalProduct + deltaX * f(i)
	839	else
	840	totalTau = totalTau + tauIncrement(i)
	841	totalProduct = totalProduct + tauIncrement(i) * f(i)
	842	end if
	843	if(totalTau >= tauLimit) exit
	844	end do
	845	weight_by_extinction = totalProduct/totalTau
	846	end function weight_by_extinction
	847	!------------------------------------------------------------------------------------------------
	848	pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size)
	849	real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
	850	real :: compute_nir_reflectance
	851
	852	real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
	853	tau, g, w0
	854	!----------------------------------------
	855	water_g(:) = get_g_nir( phaseIsLiquid, water_size)
	856	water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size)
	857	ice_g(:) = get_g_nir( phaseIsIce, ice_size)
	858	ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size)
	859	!
	860	! Combine ice and water optical properties
	861	!
	862	g(:) = 0; w0(:) = 0.
	863	tau(:) = ice_tau(:) + water_tau(:)
	864	where (tau(:) > 0)
	865	g(:) = (water_tau(:) * water_g(:) + ice_tau(:) * ice_g(:) ) / &
	866	tau(:)
	867	w0(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
	868	(g(:) * tau(:))
	869	end where
	870
	871	compute_nir_reflectance = compute_toa_reflectace(tau, g, w0)
	872	end function compute_nir_reflectance
	873	!------------------------------------------------------------------------------------------------
	874	! Retreivals
	875	!------------------------------------------------------------------------------------------------
	876	elemental function retrieve_re (phase, tau, obs_Refl_nir)
	877	integer, intent(in) :: phase
	878	real, intent(in) :: tau, obs_Refl_nir
	879	real :: retrieve_re
	880	!
	881	! Finds the re that produces the minimum mis-match between predicted and observed reflectance in
	882	! MODIS band 7 (near IR)
	883	! Uses
	884	! fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables
	885	! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0
	886	! adding-doubling for total reflectance
	887	!
	888	!
	889	!
	890	! Local variables
	891	!
	892	real, parameter :: min_distance_to_boundary = 0.01
	893	real :: re_min, re_max, delta_re
	894	integer :: i
	895
	896	real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir
	897	! --------------------------
	898
	899	if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then
	900	if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then
	901	re_min = re_water_min
	902	re_max = re_water_max
	903	trial_re(:) = trial_re_w
	904	g(:) = g_w(:)
	905	w0(:) = w0_w(:)
	906	else
	907	re_min = re_ice_min
	908	re_max = re_ice_max
	909	trial_re(:) = trial_re_i
	910	g(:) = g_i(:)
	911	w0(:) = w0_i(:)
	912	end if
	913	!
	914	! 1st attempt at index: w/coarse re resolution
	915	!
	916	predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
	917	retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
	918	!
	919	! If first retrieval works, can try 2nd iteration using greater re resolution
	920	!
[2713]	921	! DJS2015: Remove unused piece of code
	922	! if(use_two_re_iterations .and. retrieve_re > 0.) then
	923	! re_min = retrieve_re - delta_re
	924	! re_max = retrieve_re + delta_re
	925	! delta_re = (re_max - re_min)/real(num_trial_res-1)
	926	!
	927	! trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
	928	! g(:) = get_g_nir( phase, trial_re(:))
	929	! w0(:) = get_ssa_nir(phase, trial_re(:))
	930	! predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
	931	! retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
	932	! end if
	933	! DJS2015 END
[2432]	934	else
	935	retrieve_re = re_fill
	936	end if
	937
	938	end function retrieve_re
	939	! --------------------------------------------
	940	pure function interpolate_to_min(x, y, yobs)
	941	real, dimension(:), intent(in) :: x, y
	942	real, intent(in) :: yobs
	943	real :: interpolate_to_min
	944	!
	945	! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
	946	! y must be monotonic in x
	947	!
	948	real, dimension(size(x)) :: diff
	949	integer :: nPoints, minDiffLoc, lowerBound, upperBound
	950	! ---------------------------------
	951	nPoints = size(y)
	952	diff(:) = y(:) - yobs
	953	minDiffLoc = minloc(abs(diff), dim = 1)
	954
	955	if(minDiffLoc == 1) then
	956	lowerBound = minDiffLoc
	957	upperBound = minDiffLoc + 1
	958	else if(minDiffLoc == nPoints) then
	959	lowerBound = minDiffLoc - 1
	960	upperBound = minDiffLoc
	961	else
	962	if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
	963	lowerBound = minDiffLoc-1
	964	upperBound = minDiffLoc
	965	else
	966	lowerBound = minDiffLoc
	967	upperBound = minDiffLoc + 1
	968	end if
	969	end if
	970
	971	if(diff(lowerBound) * diff(upperBound) < 0) then
	972	!
	973	! Interpolate the root position linearly if we bracket the root
	974	!
	975	interpolate_to_min = x(upperBound) - &
	976	diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
	977	else
	978	interpolate_to_min = re_fill
	979	end if
	980
	981
	982	end function interpolate_to_min
	983	! --------------------------------------------
	984	! Optical properties
	985	! --------------------------------------------
	986	elemental function get_g_nir (phase, re)
	987	!
	988	! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
	989	! of size for ice and water
	990	! Fits from Steve Platnick
	991	!
	992
	993	integer, intent(in) :: phase
	994	real, intent(in) :: re
	995	real :: get_g_nir
[2713]	996
	997	real, dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
	998	small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
	999	real, dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
	1000
	1001	! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
	1002	if(phase == phaseIsLiquid) then
	1003	if(re < 7.) then
	1004	get_g_nir = fit_to_quadratic(re, small_water_coefficients)
	1005	if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
	1006	else
	1007	get_g_nir = fit_to_cubic(re, big_water_coefficients)
	1008	if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
	1009	end if
[2432]	1010	else
[2713]	1011	get_g_nir = fit_to_quadratic(re, ice_coefficients)
[2432]	1012	if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
	1013	if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
	1014	end if
	1015
	1016	end function get_g_nir
	1017
	1018	! --------------------------------------------
	1019	elemental function get_ssa_nir (phase, re)
	1020	integer, intent(in) :: phase
	1021	real, intent(in) :: re
	1022	real :: get_ssa_nir
	1023	!
	1024	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
	1025	! of size for ice and water
	1026	! Fits from Steve Platnick
	1027	!
[2713]	1028	real, dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
	1029	real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
[2432]	1030
[2713]	1031	! approx. fits from MODIS Collection 6 LUT scattering calculations
[2432]	1032	if(phase == phaseIsLiquid) then
	1033	get_ssa_nir = fit_to_quadratic(re, water_coefficients)
	1034	if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
	1035	if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
	1036	else
	1037	get_ssa_nir = fit_to_cubic(re, ice_coefficients)
	1038	if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
	1039	if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
	1040	end if
	1041
	1042	end function get_ssa_nir
	1043	! --------------------------------------------
	1044	pure function fit_to_cubic(x, coefficients)
	1045	real, intent(in) :: x
	1046	real, dimension(:), intent(in) :: coefficients
	1047	real :: fit_to_cubic
	1048
	1049
	1050	fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
	1051	end function fit_to_cubic
	1052	! --------------------------------------------
	1053	pure function fit_to_quadratic(x, coefficients)
	1054	real, intent(in) :: x
	1055	real, dimension(:), intent(in) :: coefficients
	1056	real :: fit_to_quadratic
	1057
	1058
	1059	fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
	1060	end function fit_to_quadratic
	1061	! --------------------------------------------
	1062	! Radiative transfer
	1063	! --------------------------------------------
	1064	pure function compute_toa_reflectace(tau, g, w0)
	1065	real, dimension(:), intent(in) :: tau, g, w0
	1066	real :: compute_toa_reflectace
	1067
	1068	logical, dimension(size(tau)) :: cloudMask
	1069	integer, dimension(count(tau(:) > 0)) :: cloudIndicies
	1070	real, dimension(count(tau(:) > 0)) :: Refl, Trans
	1071	real :: Refl_tot, Trans_tot
	1072	integer :: i
	1073	! ---------------------------------------
	1074	!
	1075	! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation
	1076	!
	1077	cloudMask = tau(:) > 0.
	1078	cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask)
	1079	do i = 1, size(cloudIndicies)
	1080	call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
	1081	end do
	1082
	1083	call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot)
	1084
	1085	compute_toa_reflectace = Refl_tot
	1086
	1087	end function compute_toa_reflectace
	1088	! --------------------------------------------
	1089	pure subroutine two_stream(tauint, gint, w0int, ref, tra)
	1090	real, intent(in) :: tauint, gint, w0int
	1091	real, intent(out) :: ref, tra
	1092	!
	1093	! Compute reflectance in a single layer using the two stream approximation
	1094	! The code itself is from Lazaros Oreopoulos via Steve Platnick
	1095	!
	1096	! ------------------------
	1097	! Local variables
	1098	! for delta Eddington code
	1099	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
	1100	integer, parameter :: beam = 2
	1101	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
	1102	real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
	1103	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
	1104	!
	1105	! Compute reflectance and transmittance in a single layer using the two stream approximation
	1106	! The code itself is from Lazaros Oreopoulos via Steve Platnick
	1107	!
	1108	f = gint**2
	1109	tau = (1 - w0int * f) * tauint
	1110	w0 = (1 - f) * w0int / (1 - w0int * f)
	1111	g = (gint - f) / (1 - f)
	1112
	1113	! delta-Eddington (Joseph et al. 1976)
	1114	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
	1115	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
	1116	gamma3 = (2 - 3gxmu) / 4.0
	1117	gamma4 = 1 - gamma3
	1118
	1119	if (w0int > minConservativeW0) then
	1120	! Conservative scattering
	1121	if (beam == 1) then
	1122	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
	1123	ref = rh / (1 + gamma1 * tau)
	1124	tra = 1 - ref
	1125	else if(beam == 2) then
	1126	ref = gamma1tau/(1 + gamma1tau)
	1127	tra = 1 - ref
	1128	endif
	1129	else
	1130	! Non-conservative scattering
	1131	a1 = gamma1 * gamma4 + gamma2 * gamma3
	1132	a2 = gamma1 * gamma3 + gamma2 * gamma4
	1133
	1134	rk = sqrt(gamma12 - gamma22)
	1135
	1136	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
	1137	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
	1138	r3 = 2 * rk (gamma3 - a2 xmu)
	1139	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
	1140	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
	1141
	1142	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
	1143	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
	1144	t3 = 2 * rk * (gamma4 + a1 * xmu)
	1145	t4 = r4
	1146	t5 = r5
	1147
	1148	beta = -r5 / r4
	1149
	1150	e1 = min(rk * tau, 500.)
	1151	e2 = min(tau / xmu, 500.)
	1152
	1153	if (beam == 1) then
	1154	den = r4 * exp(e1) + r5 * exp(-e1)
	1155	ref = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
	1156	den = t4 * exp(e1) + t5 * exp(-e1)
	1157	th = exp(-e2)
	1158	tra = th-thw0(t1exp(e1)-t2exp(-e1)-t3*exp(e2))/den
	1159	elseif (beam == 2) then
	1160	ef1 = exp(-e1)
	1161	ef2 = exp(-2*e1)
	1162	ref = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
	1163	tra = (2rkef1)/((rk+gamma1)(1-betaef2))
	1164	endif
	1165	end if
	1166	end subroutine two_stream
	1167	! --------------------------------------------------
	1168	elemental function two_stream_reflectance(tauint, gint, w0int)
	1169	real, intent(in) :: tauint, gint, w0int
	1170	real :: two_stream_reflectance
	1171	!
	1172	! Compute reflectance in a single layer using the two stream approximation
	1173	! The code itself is from Lazaros Oreopoulos via Steve Platnick
	1174	!
	1175	! ------------------------
	1176	! Local variables
	1177	! for delta Eddington code
	1178	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
	1179	integer, parameter :: beam = 2
	1180	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
	1181	real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
	1182	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
	1183	! ------------------------
	1184
	1185
	1186	f = gint**2
	1187	tau = (1 - w0int * f) * tauint
	1188	w0 = (1 - f) * w0int / (1 - w0int * f)
	1189	g = (gint - f) / (1 - f)
	1190
	1191	! delta-Eddington (Joseph et al. 1976)
	1192	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
	1193	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
	1194	gamma3 = (2 - 3gxmu) / 4.0
	1195	gamma4 = 1 - gamma3
	1196
	1197	if (w0int > minConservativeW0) then
	1198	! Conservative scattering
	1199	if (beam == 1) then
	1200	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
	1201	two_stream_reflectance = rh / (1 + gamma1 * tau)
	1202	elseif (beam == 2) then
	1203	two_stream_reflectance = gamma1tau/(1 + gamma1tau)
	1204	endif
	1205
	1206	else !
	1207
	1208	! Non-conservative scattering
	1209	a1 = gamma1 * gamma4 + gamma2 * gamma3
	1210	a2 = gamma1 * gamma3 + gamma2 * gamma4
	1211
	1212	rk = sqrt(gamma12 - gamma22)
	1213
	1214	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
	1215	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
	1216	r3 = 2 * rk (gamma3 - a2 xmu)
	1217	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
	1218	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
	1219
	1220	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
	1221	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
	1222	t3 = 2 * rk * (gamma4 + a1 * xmu)
	1223	t4 = r4
	1224	t5 = r5
	1225
	1226	beta = -r5 / r4
	1227
	1228	e1 = min(rk * tau, 500.)
	1229	e2 = min(tau / xmu, 500.)
	1230
	1231	if (beam == 1) then
	1232	den = r4 * exp(e1) + r5 * exp(-e1)
	1233	two_stream_reflectance = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
	1234	elseif (beam == 2) then
	1235	ef1 = exp(-e1)
	1236	ef2 = exp(-2*e1)
	1237	two_stream_reflectance = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
	1238	endif
	1239
	1240	end if
	1241	end function two_stream_reflectance
	1242	! --------------------------------------------
	1243	pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot)
	1244	real, dimension(:), intent(in) :: Refl, Tran
	1245	real, intent(out) :: Refl_tot, Tran_tot
	1246	!
	1247	! Use adding/doubling formulas to compute total reflectance and transmittance from layer values
	1248	!
	1249
	1250	integer :: i
	1251	real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative
	1252
	1253	Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)
	1254
	1255	do i=2, size(Refl)
	1256	! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
	1257	Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)(Tran_cumulative(i-1)2)/(1 - Refl_cumulative(i-1) Refl(i))
	1258	Tran_cumulative(i) = (Tran_cumulative(i-1)Tran(i)) / (1 - Refl_cumulative(i-1) Refl(i))
	1259	end do
	1260
	1261	Refl_tot = Refl_cumulative(size(Refl))
	1262	Tran_tot = Tran_cumulative(size(Refl))
	1263
	1264	end subroutine adding_doubling
	1265	! --------------------------------------------------
	1266	subroutine complain_and_die(message)
	1267	character(len = *), intent(in) :: message
	1268
	1269	write(6, *) "Failure in MODIS simulator"
	1270	write(6, *) trim(message)
	1271	stop
	1272	end subroutine complain_and_die
	1273	!------------------------------------------------------------------------------------------------
[2713]	1274	end module mod_modis_sim

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

Download in other formats:

Original Format