Context Navigation

modis_simulator.F90

Last change on this file was 2720, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2664:2719 into testing branch
File size: 63.2 KB

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

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

Download in other formats:

Original Format