Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mod_modis_sim.f90 @ 5445

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

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

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

Download in other formats:

Original Format