Context Navigation

modis_simulator.F90 @ 3184

Last change on this file since 3184 was 3184, checked in by idelkadi, 6 years ago
Correction de bogue dans le simulateur Modis (mise a jour dans la version Cosp v1.4.2)
File size: 63.3 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.
535	! ########################################################################################
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
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
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	w0(:) = (water_tau(:) * water_w0(:) + ice_tau(:) * ice_w0(:)) / &
866	tau(:)
867	g(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
868	(w0(:) * 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	!
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
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
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
1010	else
1011	get_g_nir = fit_to_quadratic(re, ice_coefficients)
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	!
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 /)
1030
1031	! approx. fits from MODIS Collection 6 LUT scattering calculations
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	!------------------------------------------------------------------------------------------------
1274	end module mod_modis_sim

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

Download in other formats:

Original Format