Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

modis_simulator.F90 @ 3331

Last change on this file since 3331 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 48.5 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 15:08:38 +0100 (mer. 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	logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations?
82
83	!
84	! Precompute near-IR optical params vs size for retrieval scheme
85	!
86	integer, private :: i
87	real, dimension(num_trial_res), parameter :: &
88	trial_re_w = re_water_min + (re_water_max - re_water_min)/(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /), &
89	trial_re_i = re_ice_min + (re_ice_max - re_ice_min) /(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /)
90
91	! Can't initialze these during compilation, but do in before looping columns in retrievals
92	real, dimension(num_trial_res) :: g_w, g_i, w0_w, w0_i
93	! ------------------------------
94	! Bin boundaries for the joint optical thickness/cloud top pressure histogram
95	!
96	integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7
97
98	real, private :: dummy_real
99	real, dimension(numTauHistogramBins + 1), parameter :: &
100	tauHistogramBoundaries = (/ min_OpticalThickness, 1.3, 3.6, 9.4, 23., 60., 10000. /)
101	real, dimension(numPressureHistogramBins + 1), parameter :: & ! Units Pa
102	pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /)
103	real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100.
104	!
105	! For output - nominal bin centers and bin boundaries. On output pressure bins are highest to lowest.
106	!
107	integer, private :: k, l
108	real, parameter, dimension(2, numTauHistogramBins) :: &
109	nominalTauHistogramBoundaries = &
110	reshape(source = (/ tauHistogramBoundaries(1), &
111	((tauHistogramBoundaries(k), l = 1, 2), k = 2, numTauHistogramBins), &
112	100000. /), &
113	shape = (/2, numTauHistogramBins /) )
114	real, parameter, dimension(numTauHistogramBins) :: &
115	nominalTauHistogramCenters = (nominalTauHistogramBoundaries(1, :) + &
116	nominalTauHistogramBoundaries(2, :) ) / 2.
117
118	real, parameter, dimension(2, numPressureHistogramBins) :: &
119	nominalPressureHistogramBoundaries = &
120	reshape(source = (/ 100000., &
121	((pressureHistogramBoundaries(k), l = 1, 2), k = numPressureHistogramBins, 2, -1), &
122	0. /), &
123	shape = (/2, numPressureHistogramBins /) )
124	real, parameter, dimension(numPressureHistogramBins) :: &
125	nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + &
126	nominalPressureHistogramBoundaries(2, :) ) / 2.
127	! ------------------------------
128	! There are two ways to call the MODIS simulator:
129	! 1) Provide total optical thickness and liquid/ice water content and we'll partition tau in
130	! subroutine modis_L2_simulator_oneTau, or
131	! 2) Provide ice and liquid optical depths in each layer
132	!
133	interface modis_L2_simulator
134	module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus
135	end interface
136	contains
137	!------------------------------------------------------------------------------------------------
138	! MODIS simulator using specified liquid and ice optical thickness in each layer
139	!
140	! Note: this simulator operates on all points; to match MODIS itself night-time
141	! points should be excluded
142	!
143	! Note: the simulator requires as input the optical thickness and cloud top pressure
144	! derived from the ISCCP simulator run with parameter top_height = 1.
145	! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing
146	! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that
147	! alogrithm in this simulator we simply report the values from the ISCCP simulator.
148	!
149	subroutine modis_L2_simulator_twoTaus( &
150	temp, pressureLayers, pressureLevels, &
151	liquid_opticalThickness, ice_opticalThickness, &
152	waterSize, iceSize, &
153	isccpTau, isccpCloudTopPressure, &
154	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
155
156	! Grid-mean quantities at layer centers, starting at the model top
157	! dimension nLayers
158	real, dimension(:), intent(in ) :: temp, & ! Temperature, K
159	pressureLayers, & ! Pressure, Pa
160	pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1)
161	! Sub-column quantities
162	! dimension nSubcols, nLayers
163	real, dimension(:, :), intent(in ) :: liquid_opticalThickness, & ! Layer optical thickness @ 0.67 microns due to liquid
164	ice_opticalThickness ! ditto, due to ice
165	real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns
166	iceSize ! Cloud ice effective radius, microns
167
168	! Cloud properties retrieved from ISCCP using top_height = 1
169	! dimension nSubcols
170	real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness
171	isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa)
172
173	! Properties retrieved by MODIS
174	! dimension nSubcols
175	integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer, defined in module header
176	real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
177	retrievedTau, & ! unitless
178	retrievedSize ! microns
179	! ---------------------------------------------------
180	! Local variables
181	logical, dimension(size(retrievedTau)) :: cloudMask
182	real, dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal
183	real :: integratedLiquidFraction
184	integer :: i, nSubcols, nLevels
185
186	! ---------------------------------------------------
187	nSubcols = size(liquid_opticalThickness, 1)
188	nLevels = size(liquid_opticalThickness, 2)
189
190	!
191	! Initial error checks
192	!
193	if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), &
194	size(isccpTau), size(isccpCloudTopPressure), &
195	size(retrievedPhase), size(retrievedCloudTopPressure), &
196	size(retrievedTau), size(retrievedSize) /) /= nSubcols )) &
197	call complain_and_die("Differing number of subcolumns in one or more arrays")
198
199	if(any((/ size(ice_opticalThickness, 2), size(waterSize, 2), size(iceSize, 2), &
200	size(temp), size(pressureLayers), size(pressureLevels)-1 /) /= nLevels )) &
201	call complain_and_die("Differing number of levels in one or more arrays")
202
203
204	if(any( (/ any(temp <= 0.), any(pressureLayers <= 0.), &
205	any(liquid_opticalThickness < 0.), &
206	any(ice_opticalThickness < 0.), &
207	any(waterSize < 0.), any(iceSize < 0.) /) )) &
208	call complain_and_die("Input values out of bounds")
209
210	! ---------------------------------------------------
211	!
212	! Compute the total optical thickness and the proportion due to liquid in each cell
213	!
214	where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.)
215	tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :))
216	elsewhere
217	tauLiquidFraction(:, :) = 0.
218	end where
219	tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :)
220
221	!
222	! Optical depth retrieval
223	! This is simply a sum over the optical thickness in each layer
224	! It should agree with the ISCCP values after min values have been excluded
225	!
226	retrievedTau(:) = sum(tauTotal(:, :), dim = 2)
227
228	!
229	! Cloud detection - does optical thickness exceed detection threshold?
230	!
231	cloudMask = retrievedTau(:) >= min_OpticalThickness
232
233	!
234	! Initialize initial estimates for size retrievals
235	!
236	if(any(cloudMask) .and. .not. useSimpleReScheme) then
237	g_w(:) = get_g_nir( phaseIsLiquid, trial_re_w(:))
238	w0_w(:) = get_ssa_nir(phaseIsLiquid, trial_re_w(:))
239	g_i(:) = get_g_nir( phaseIsIce, trial_re_i(:))
240	w0_i(:) = get_ssa_nir(phaseIsIce, trial_re_i(:))
241	end if
242
243	do i = 1, nSubCols
244	if(cloudMask(i)) then
245	!
246	! Cloud top pressure determination
247	! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds
248	! lower than that.
249	! For CO2 slicing we report the optical-depth weighted pressure, integrating to a specified
250	! optical depth
251	! This assumes linear variation in p between levels. Linear in ln(p) is probably better,
252	! though we'd need to deal with the lowest pressure gracefully.
253	!
254	retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), &
255	pressureLevels, &
256	CO2Slicing_TauLimit)
257
258
259	!
260	! Phase determination - determine fraction of total tau that's liquid
261	! When ice and water contribute about equally to the extinction we can't tell
262	! what the phase is
263	!
264	integratedLiquidFraction = weight_by_extinction(tauTotal(i, :), &
265	tauLiquidFraction(i, :), &
266	phase_TauLimit)
267	if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then
268	retrievedPhase(i) = phaseIsLiquid
269	else if (integratedLiquidFraction <= 1.- phaseDiscrimination_Threshold) then
270	retrievedPhase(i) = phaseIsIce
271	else
272	retrievedPhase(i) = phaseIsUndetermined
273	end if
274
275	!
276	! Size determination
277	!
278	if(useSimpleReScheme) then
279	! This is the extinction-weighted size considering only the phase we've chosen
280	!
281	if(retrievedPhase(i) == phaseIsIce) then
282	retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :), &
283	iceSize(i, :), &
284	(1. - integratedLiquidFraction) * size_TauLimit)
285
286	else if(retrievedPhase(i) == phaseIsLiquid) then
287	retrievedSize(i) = weight_by_extinction(liquid_opticalThickness(i, :), &
288	waterSize(i, :), &
289	integratedLiquidFraction * size_TauLimit)
290
291	else
292	retrievedSize(i) = 0.
293	end if
294	else
295	retrievedSize(i) = 1.0e-06*retrieve_re(retrievedPhase(i), retrievedTau(i), &
296	obs_Refl_nir = compute_nir_reflectance(liquid_opticalThickness(i, :), waterSize(i, :)*1.0e6, &
297	ice_opticalThickness(i, :), iceSize(i, :)*1.0e6))
298	end if
299	else
300	!
301	! Values when we don't think there's a cloud.
302	!
303	retrievedCloudTopPressure(i) = R_UNDEF
304	retrievedPhase(i) = phaseIsNone
305	retrievedSize(i) = R_UNDEF
306	retrievedTau(i) = R_UNDEF
307	end if
308	end do
309	where((retrievedSize(:) < 0.).and.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill
310
311	! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS
312	! mimics what MODIS does to first order.
313	! Of course, ISCCP cloud top pressures are in mb.
314	!
315	where(cloudMask(:) .and. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) &
316	retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100.
317
318	end subroutine modis_L2_simulator_twoTaus
319	!------------------------------------------------------------------------------------------------
320	!
321	! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents;
322	! we'll partition this into ice and liquid optical thickness and call the full MODIS simulator
323	!
324	subroutine modis_L2_simulator_oneTau( &
325	temp, pressureLayers, pressureLevels, &
326	opticalThickness, cloudWater, cloudIce, &
327	waterSize, iceSize, &
328	isccpTau, isccpCloudTopPressure, &
329	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
330	! Grid-mean quantities at layer centers,
331	! dimension nLayers
332	real, dimension(:), intent(in ) :: temp, & ! Temperature, K
333	pressureLayers, & ! Pressure, Pa
334	pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1)
335	! Sub-column quantities
336	! dimension nLayers, nSubcols
337	real, dimension(:, :), intent(in ) :: opticalThickness, & ! Layer optical thickness @ 0.67 microns
338	cloudWater, & ! Cloud water content, arbitrary units
339	cloudIce ! Cloud water content, same units as cloudWater
340	real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns
341	iceSize ! Cloud ice effective radius, microns
342
343	! Cloud properties retrieved from ISCCP using top_height = 1
344	! dimension nSubcols
345
346	real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness
347	isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa)
348
349	! Properties retrieved by MODIS
350	! dimension nSubcols
351	integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer
352	real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
353	retrievedTau, & ! unitless
354	retrievedSize ! microns (or whatever units
355	! waterSize and iceSize are supplied in)
356	! ---------------------------------------------------
357	! Local variables
358	real, dimension(size(opticalThickness, 1), size(opticalThickness, 2)) :: &
359	liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction
360
361	! ---------------------------------------------------
362
363	where(cloudIce(:, :) <= 0.)
364	tauLiquidFraction(:, :) = 1.
365	elsewhere
366	where (cloudWater(:, :) <= 0.)
367	tauLiquidFraction(:, :) = 0.
368	elsewhere
369	!
370	! Geometic optics limit - tau as LWP/re (proportional to LWC/re)
371	!
372	tauLiquidFraction(:, :) = (cloudWater(:, :)/waterSize(:, :)) / &
373	(cloudWater(:, :)/waterSize(:, :) + cloudIce(:, :)/(ice_density * iceSize(:, :)) )
374	end where
375	end where
376	liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :)
377	ice_opticalThickness (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :)
378
379	call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels, &
380	liquid_opticalThickness, ice_opticalThickness, &
381	waterSize, iceSize, &
382	isccpTau, isccpCloudTopPressure, &
383	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
384
385	end subroutine modis_L2_simulator_oneTau
386	!------------------------------------------------------------------------------------------------
387	subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, &
388	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
389	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
390	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
391	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
392	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
393	Cloud_Top_Pressure_Total_Mean, &
394	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
395	Optical_Thickness_vs_Cloud_Top_Pressure)
396	!
397	! Inputs; dimension nPoints, nSubcols
398	!
399	integer, dimension(:, :), intent(in) :: phase
400	real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size
401	!
402	! Outputs; dimension nPoints
403	!
404	real, dimension(:), intent(out) :: &
405	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
406	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
407	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
408	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
409	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
410	Cloud_Top_Pressure_Total_Mean, &
411	Liquid_Water_Path_Mean, Ice_Water_Path_Mean
412	! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins
413	real, dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure
414	! ---------------------------
415	! Local variables
416	!
417	real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units
418	integer :: i, j
419	integer :: nPoints, nSubcols
420	logical, dimension(size(phase, 1), size(phase, 2)) :: &
421	cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask
422	logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins ) :: tauMask
423	logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask
424	! ---------------------------
425
426	nPoints = size(phase, 1)
427	nSubcols = size(phase, 2)
428	!
429	! Array conformance checks
430	!
431	if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), &
432	size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), &
433	size(Cloud_Fraction_High_Mean), size(Cloud_Fraction_Mid_Mean), size(Cloud_Fraction_Low_Mean), &
434	size(Optical_Thickness_Total_Mean), size(Optical_Thickness_Water_Mean), size(Optical_Thickness_Ice_Mean), &
435	size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), &
436	size(Optical_Thickness_Ice_MeanLog10), size(Cloud_Particle_Size_Water_Mean), &
437	size(Cloud_Particle_Size_Ice_Mean), size(Cloud_Top_Pressure_Total_Mean), &
438	size(Liquid_Water_Path_Mean), size(Ice_Water_Path_Mean) /) /= nPoints)) &
439	call complain_and_die("Some L3 arrays have wrong number of grid points")
440	if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) &
441	call complain_and_die("Some L3 arrays have wrong number of subcolumns")
442
443
444	!
445	! Include only those pixels with successful retrievals in the statistics
446	!
447	validRetrievalMask(:, :) = particle_size(:, :) > 0.
448	cloudMask = phase(:, :) /= phaseIsNone .and. validRetrievalMask(:, :)
449	waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :)
450	iceCloudMask = phase(:, :) == phaseIsIce .and. validRetrievalMask(:, :)
451	!
452	! Use these as pixel counts at first
453	!
454	Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2))
455	Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2))
456	Cloud_Fraction_Ice_Mean(:) = real(count(iceCloudMask, dim = 2))
457
458	Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2))
459	Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .and. cloud_top_pressure > lowCloudPressureLimit, dim = 2))
460	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:)
461
462	!
463	! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0.
464	!
465	where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1.
466	where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1.
467	where (Cloud_Fraction_Ice_Mean == 0) Cloud_Fraction_Ice_Mean = -1.
468
469	Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask, dim = 2) / Cloud_Fraction_Total_Mean(:)
470	Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
471	Optical_Thickness_Ice_Mean = sum(optical_thickness, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
472
473	! We take the absolute value of optical thickness here to satisfy compilers that complains when we
474	! evaluate the logarithm of a negative number, even though it's not included in the sum.
475	Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask, dim = 2) / &
476	Cloud_Fraction_Total_Mean(:)
477	Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / &
478	Cloud_Fraction_Water_Mean(:)
479	Optical_Thickness_Ice_MeanLog10 = sum(log10(abs(optical_thickness)), mask = iceCloudMask, dim = 2) / &
480	Cloud_Fraction_Ice_Mean(:)
481
482	Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
483	Cloud_Particle_Size_Ice_Mean = sum(particle_size, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
484
485	Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2))
486
487	Liquid_Water_Path_Mean = LWP_conversion &
488	* sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) &
489	/ Cloud_Fraction_Water_Mean(:)
490	Ice_Water_Path_Mean = LWP_conversion * ice_density &
491	* sum(particle_size * optical_thickness, mask = iceCloudMask, dim = 2) &
492	/ Cloud_Fraction_Ice_Mean(:)
493
494	!
495	! Normalize pixel counts to fraction
496	! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0.
497	!
498	Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols)
499	Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols)
500	Cloud_Fraction_Ice_Mean(:) = max(0., Cloud_Fraction_Ice_Mean(:) /nSubcols)
501
502	Cloud_Fraction_High_Mean(:) = Cloud_Fraction_High_Mean(:) /nSubcols
503	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Mid_Mean(:) /nSubcols
504	Cloud_Fraction_Low_Mean(:) = Cloud_Fraction_Low_Mean(:) /nSubcols
505
506	! ----
507	! Joint histogram
508	!
509	do i = 1, numTauHistogramBins
510	where(cloudMask(:, :))
511	tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. &
512	optical_thickness(:, :) < tauHistogramBoundaries(i+1)
513	elsewhere
514	tauMask(:, :, i) = .false.
515	end where
516	end do
517
518	do i = 1, numPressureHistogramBins
519	where(cloudMask(:, :))
520	pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. &
521	cloud_top_pressure(:, :) < pressureHistogramBoundaries(i+1)
522	elsewhere
523	pressureMask(:, :, i) = .false.
524	end where
525	end do
526
527	do i = 1, numPressureHistogramBins
528	do j = 1, numTauHistogramBins
529	Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = &
530	real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
531	end do
532	end do
533
534	end subroutine modis_L3_simulator
535	!------------------------------------------------------------------------------------------------
536	function cloud_top_pressure(tauIncrement, pressure, tauLimit)
537	real, dimension(:), intent(in) :: tauIncrement, pressure
538	real, intent(in) :: tauLimit
539	real :: cloud_top_pressure
540	!
541	! Find the extinction-weighted pressure. Assume that pressure varies linearly between
542	! layers and use the trapezoidal rule.
543	!
544
545	real :: deltaX, totalTau, totalProduct
546	integer :: i
547
548	totalTau = 0.; totalProduct = 0.
549	do i = 2, size(tauIncrement)
550	if(totalTau + tauIncrement(i) > tauLimit) then
551	deltaX = tauLimit - totalTau
552	totalTau = totalTau + deltaX
553	!
554	! Result for trapezoidal rule when you take less than a full step
555	! tauIncrement is a layer-integrated value
556	!
557	totalProduct = totalProduct &
558	+ pressure(i-1) * deltaX &
559	+ (pressure(i) - pressure(i-1)) * deltaX*2/(2. tauIncrement(i))
560	else
561	totalTau = totalTau + tauIncrement(i)
562	totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2.
563	end if
564	if(totalTau >= tauLimit) exit
565	end do
566	cloud_top_pressure = totalProduct/totalTau
567	end function cloud_top_pressure
568	!------------------------------------------------------------------------------------------------
569	function weight_by_extinction(tauIncrement, f, tauLimit)
570	real, dimension(:), intent(in) :: tauIncrement, f
571	real, intent(in) :: tauLimit
572	real :: weight_by_extinction
573	!
574	! Find the extinction-weighted value of f(tau), assuming constant f within each layer
575	!
576
577	real :: deltaX, totalTau, totalProduct
578	integer :: i
579
580	totalTau = 0.; totalProduct = 0.
581	do i = 1, size(tauIncrement)
582	if(totalTau + tauIncrement(i) > tauLimit) then
583	deltaX = tauLimit - totalTau
584	totalTau = totalTau + deltaX
585	totalProduct = totalProduct + deltaX * f(i)
586	else
587	totalTau = totalTau + tauIncrement(i)
588	totalProduct = totalProduct + tauIncrement(i) * f(i)
589	end if
590	if(totalTau >= tauLimit) exit
591	end do
592	weight_by_extinction = totalProduct/totalTau
593	end function weight_by_extinction
594	!------------------------------------------------------------------------------------------------
595	pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size)
596	real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
597	real :: compute_nir_reflectance
598
599	real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
600	tau, g, w0
601	!----------------------------------------
602	water_g(:) = get_g_nir( phaseIsLiquid, water_size)
603	water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size)
604	ice_g(:) = get_g_nir( phaseIsIce, ice_size)
605	ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size)
606	!
607	! Combine ice and water optical properties
608	!
609	g(:) = 0; w0(:) = 0.
610	tau(:) = ice_tau(:) + water_tau(:)
611	where (tau(:) > 0)
612	g(:) = (water_tau(:) * water_g(:) + ice_tau(:) * ice_g(:) ) / &
613	tau(:)
614	w0(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
615	(g(:) * tau(:))
616	end where
617
618	compute_nir_reflectance = compute_toa_reflectace(tau, g, w0)
619	end function compute_nir_reflectance
620	!------------------------------------------------------------------------------------------------
621	! Retreivals
622	!------------------------------------------------------------------------------------------------
623	elemental function retrieve_re (phase, tau, obs_Refl_nir)
624	integer, intent(in) :: phase
625	real, intent(in) :: tau, obs_Refl_nir
626	real :: retrieve_re
627	!
628	! Finds the re that produces the minimum mis-match between predicted and observed reflectance in
629	! MODIS band 7 (near IR)
630	! Uses
631	! fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables
632	! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0
633	! adding-doubling for total reflectance
634	!
635	!
636	!
637	! Local variables
638	!
639	real, parameter :: min_distance_to_boundary = 0.01
640	real :: re_min, re_max, delta_re
641	integer :: i
642
643	real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir
644	! --------------------------
645
646	if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then
647	if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then
648	re_min = re_water_min
649	re_max = re_water_max
650	trial_re(:) = trial_re_w
651	g(:) = g_w(:)
652	w0(:) = w0_w(:)
653	else
654	re_min = re_ice_min
655	re_max = re_ice_max
656	trial_re(:) = trial_re_i
657	g(:) = g_i(:)
658	w0(:) = w0_i(:)
659	end if
660	!
661	! 1st attempt at index: w/coarse re resolution
662	!
663	predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
664	retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
665	!
666	! If first retrieval works, can try 2nd iteration using greater re resolution
667	!
668	if(use_two_re_iterations .and. retrieve_re > 0.) then
669	re_min = retrieve_re - delta_re
670	re_max = retrieve_re + delta_re
671	delta_re = (re_max - re_min)/real(num_trial_res-1)
672
673	trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
674	g(:) = get_g_nir( phase, trial_re(:))
675	w0(:) = get_ssa_nir(phase, trial_re(:))
676	predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
677	retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
678	end if
679	else
680	retrieve_re = re_fill
681	end if
682
683	end function retrieve_re
684	! --------------------------------------------
685	pure function interpolate_to_min(x, y, yobs)
686	real, dimension(:), intent(in) :: x, y
687	real, intent(in) :: yobs
688	real :: interpolate_to_min
689	!
690	! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
691	! y must be monotonic in x
692	!
693	real, dimension(size(x)) :: diff
694	integer :: nPoints, minDiffLoc, lowerBound, upperBound
695	! ---------------------------------
696	nPoints = size(y)
697	diff(:) = y(:) - yobs
698	minDiffLoc = minloc(abs(diff), dim = 1)
699
700	if(minDiffLoc == 1) then
701	lowerBound = minDiffLoc
702	upperBound = minDiffLoc + 1
703	else if(minDiffLoc == nPoints) then
704	lowerBound = minDiffLoc - 1
705	upperBound = minDiffLoc
706	else
707	if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
708	lowerBound = minDiffLoc-1
709	upperBound = minDiffLoc
710	else
711	lowerBound = minDiffLoc
712	upperBound = minDiffLoc + 1
713	end if
714	end if
715
716	if(diff(lowerBound) * diff(upperBound) < 0) then
717	!
718	! Interpolate the root position linearly if we bracket the root
719	!
720	interpolate_to_min = x(upperBound) - &
721	diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
722	else
723	interpolate_to_min = re_fill
724	end if
725
726
727	end function interpolate_to_min
728	! --------------------------------------------
729	! Optical properties
730	! --------------------------------------------
731	elemental function get_g_nir (phase, re)
732	!
733	! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
734	! of size for ice and water
735	! Fits from Steve Platnick
736	!
737
738	integer, intent(in) :: phase
739	real, intent(in) :: re
740	real :: get_g_nir
741
742	real, dimension(3), parameter :: ice_coefficients = (/ 0.7432, 4.5563e-3, -2.8697e-5 /), &
743	small_water_coefficients = (/ 0.8027, -1.0496e-2, 1.7071e-3 /), &
744	big_water_coefficients = (/ 0.7931, 5.3087e-3, -7.4995e-5 /)
745
746	! approx. fits from MODIS Collection 5 LUT scattering calculations
747	if(phase == phaseIsLiquid) then
748	if(re < 8.) then
749	get_g_nir = fit_to_quadratic(re, small_water_coefficients)
750	if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
751	else
752	get_g_nir = fit_to_quadratic(re, big_water_coefficients)
753	if(re > re_water_max) get_g_nir = fit_to_quadratic(re_water_max, big_water_coefficients)
754	end if
755	else
756	get_g_nir = fit_to_quadratic(re, ice_coefficients)
757	if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
758	if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
759	end if
760
761	end function get_g_nir
762
763	! --------------------------------------------
764	elemental function get_ssa_nir (phase, re)
765	integer, intent(in) :: phase
766	real, intent(in) :: re
767	real :: get_ssa_nir
768	!
769	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
770	! of size for ice and water
771	! Fits from Steve Platnick
772	!
773
774	real, dimension(4), parameter :: ice_coefficients = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /)
775	real, dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /)
776
777	! approx. fits from MODIS Collection 5 LUT scattering calculations
778	if(phase == phaseIsLiquid) then
779	get_ssa_nir = fit_to_quadratic(re, water_coefficients)
780	if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
781	if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
782	else
783	get_ssa_nir = fit_to_cubic(re, ice_coefficients)
784	if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
785	if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
786	end if
787
788	end function get_ssa_nir
789	! --------------------------------------------
790	pure function fit_to_cubic(x, coefficients)
791	real, intent(in) :: x
792	real, dimension(:), intent(in) :: coefficients
793	real :: fit_to_cubic
794
795
796	fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
797	end function fit_to_cubic
798	! --------------------------------------------
799	pure function fit_to_quadratic(x, coefficients)
800	real, intent(in) :: x
801	real, dimension(:), intent(in) :: coefficients
802	real :: fit_to_quadratic
803
804
805	fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
806	end function fit_to_quadratic
807	! --------------------------------------------
808	! Radiative transfer
809	! --------------------------------------------
810	pure function compute_toa_reflectace(tau, g, w0)
811	real, dimension(:), intent(in) :: tau, g, w0
812	real :: compute_toa_reflectace
813
814	logical, dimension(size(tau)) :: cloudMask
815	integer, dimension(count(tau(:) > 0)) :: cloudIndicies
816	real, dimension(count(tau(:) > 0)) :: Refl, Trans
817	real :: Refl_tot, Trans_tot
818	integer :: i
819	! ---------------------------------------
820	!
821	! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation
822	!
823	cloudMask = tau(:) > 0.
824	cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask)
825	do i = 1, size(cloudIndicies)
826	call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
827	end do
828
829	call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot)
830
831	compute_toa_reflectace = Refl_tot
832
833	end function compute_toa_reflectace
834	! --------------------------------------------
835	pure subroutine two_stream(tauint, gint, w0int, ref, tra)
836	real, intent(in) :: tauint, gint, w0int
837	real, intent(out) :: ref, tra
838	!
839	! Compute reflectance in a single layer using the two stream approximation
840	! The code itself is from Lazaros Oreopoulos via Steve Platnick
841	!
842	! ------------------------
843	! Local variables
844	! for delta Eddington code
845	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
846	integer, parameter :: beam = 2
847	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
848	real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
849	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
850	!
851	! Compute reflectance and transmittance in a single layer using the two stream approximation
852	! The code itself is from Lazaros Oreopoulos via Steve Platnick
853	!
854	f = gint**2
855	tau = (1 - w0int * f) * tauint
856	w0 = (1 - f) * w0int / (1 - w0int * f)
857	g = (gint - f) / (1 - f)
858
859	! delta-Eddington (Joseph et al. 1976)
860	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
861	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
862	gamma3 = (2 - 3gxmu) / 4.0
863	gamma4 = 1 - gamma3
864
865	if (w0int > minConservativeW0) then
866	! Conservative scattering
867	if (beam == 1) then
868	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
869	ref = rh / (1 + gamma1 * tau)
870	tra = 1 - ref
871	else if(beam == 2) then
872	ref = gamma1tau/(1 + gamma1tau)
873	tra = 1 - ref
874	endif
875	else
876	! Non-conservative scattering
877	a1 = gamma1 * gamma4 + gamma2 * gamma3
878	a2 = gamma1 * gamma3 + gamma2 * gamma4
879
880	rk = sqrt(gamma12 - gamma22)
881
882	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
883	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
884	r3 = 2 * rk (gamma3 - a2 xmu)
885	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
886	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
887
888	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
889	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
890	t3 = 2 * rk * (gamma4 + a1 * xmu)
891	t4 = r4
892	t5 = r5
893
894	beta = -r5 / r4
895
896	e1 = min(rk * tau, 500.)
897	e2 = min(tau / xmu, 500.)
898
899	if (beam == 1) then
900	den = r4 * exp(e1) + r5 * exp(-e1)
901	ref = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
902	den = t4 * exp(e1) + t5 * exp(-e1)
903	th = exp(-e2)
904	tra = th-thw0(t1exp(e1)-t2exp(-e1)-t3*exp(e2))/den
905	elseif (beam == 2) then
906	ef1 = exp(-e1)
907	ef2 = exp(-2*e1)
908	ref = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
909	tra = (2rkef1)/((rk+gamma1)(1-betaef2))
910	endif
911	end if
912	end subroutine two_stream
913	! --------------------------------------------------
914	elemental function two_stream_reflectance(tauint, gint, w0int)
915	real, intent(in) :: tauint, gint, w0int
916	real :: two_stream_reflectance
917	!
918	! Compute reflectance in a single layer using the two stream approximation
919	! The code itself is from Lazaros Oreopoulos via Steve Platnick
920	!
921	! ------------------------
922	! Local variables
923	! for delta Eddington code
924	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
925	integer, parameter :: beam = 2
926	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
927	real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
928	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
929	! ------------------------
930
931
932	f = gint**2
933	tau = (1 - w0int * f) * tauint
934	w0 = (1 - f) * w0int / (1 - w0int * f)
935	g = (gint - f) / (1 - f)
936
937	! delta-Eddington (Joseph et al. 1976)
938	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
939	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
940	gamma3 = (2 - 3gxmu) / 4.0
941	gamma4 = 1 - gamma3
942
943	if (w0int > minConservativeW0) then
944	! Conservative scattering
945	if (beam == 1) then
946	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
947	two_stream_reflectance = rh / (1 + gamma1 * tau)
948	elseif (beam == 2) then
949	two_stream_reflectance = gamma1tau/(1 + gamma1tau)
950	endif
951
952	else !
953
954	! Non-conservative scattering
955	a1 = gamma1 * gamma4 + gamma2 * gamma3
956	a2 = gamma1 * gamma3 + gamma2 * gamma4
957
958	rk = sqrt(gamma12 - gamma22)
959
960	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
961	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
962	r3 = 2 * rk (gamma3 - a2 xmu)
963	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
964	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
965
966	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
967	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
968	t3 = 2 * rk * (gamma4 + a1 * xmu)
969	t4 = r4
970	t5 = r5
971
972	beta = -r5 / r4
973
974	e1 = min(rk * tau, 500.)
975	e2 = min(tau / xmu, 500.)
976
977	if (beam == 1) then
978	den = r4 * exp(e1) + r5 * exp(-e1)
979	two_stream_reflectance = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
980	elseif (beam == 2) then
981	ef1 = exp(-e1)
982	ef2 = exp(-2*e1)
983	two_stream_reflectance = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
984	endif
985
986	end if
987	end function two_stream_reflectance
988	! --------------------------------------------
989	pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot)
990	real, dimension(:), intent(in) :: Refl, Tran
991	real, intent(out) :: Refl_tot, Tran_tot
992	!
993	! Use adding/doubling formulas to compute total reflectance and transmittance from layer values
994	!
995
996	integer :: i
997	real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative
998
999	Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)
1000
1001	do i=2, size(Refl)
1002	! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
1003	Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)(Tran_cumulative(i-1)2)/(1 - Refl_cumulative(i-1) Refl(i))
1004	Tran_cumulative(i) = (Tran_cumulative(i-1)Tran(i)) / (1 - Refl_cumulative(i-1) Refl(i))
1005	end do
1006
1007	Refl_tot = Refl_cumulative(size(Refl))
1008	Tran_tot = Tran_cumulative(size(Refl))
1009
1010	end subroutine adding_doubling
1011	! --------------------------------------------------
1012	subroutine complain_and_die(message)
1013	character(len = *), intent(in) :: message
1014
1015	write(6, *) "Failure in MODIS simulator"
1016	write(6, *) trim(message)
1017	stop
1018	end subroutine complain_and_die
1019	!------------------------------------------------------------------------------------------------
1020	end module mod_modis_sim

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

Download in other formats:

Original Format