Context Navigation

mod_modis_sim.F90 @ 5441

Last change on this file since 5441 was 5185, checked in by abarral, 3 months ago
Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
File size: 63.6 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	! History:
29	! May 2009 - Robert Pincus - Initial version
30	! June 2009 - Steve Platnick and Robert Pincus - Simple radiative transfer for size retrievals
31	! August 2009 - Robert Pincus - Consistency and bug fixes suggested by Rick Hemler (GFDL)
32	! November 2009 - Robert Pincus - Bux fixes and speed-ups after experience with Rick Hemler using AM2 (GFDL)
33	! January 2010 - Robert Pincus - Added high, middle, low cloud fractions
34
35	! Notes on using the MODIS simulator:
36	! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 microns, or
37	! optical thickness at 0.67 microns and ice- and liquid-water contents (in consistent units of
38	! your choosing)
39	! *) Required input also includes the optical thickness and cloud top pressure
40	! derived from the ISCCP simulator run with parameter top_height = 1.
41	! *) Cloud particle sizes are specified as radii, measured in meters, though within the module we
42	! use units of microns. Where particle sizes are outside the bounds used in the MODIS retrieval
43	! libraries (parameters re_water_min, re_ice_min, etc.) the simulator returns missing values (re_fill)
44
45	! When error conditions are encountered this code calls the function complain_and_die, supplied at the
46	! bottom of this module. Users probably want to replace this with something more graceful.
47
48	module mod_modis_sim
49	USE MOD_COSP_TYPES, only: R_UNDEF, ok_debug_cosp
50	implicit none
51	! ------------------------------
52	! Algorithmic parameters
53
54	real, parameter :: ice_density = 0.93 ! liquid density is 1.
55
56	! Retrieval parameters
57
58	real, parameter :: min_OpticalThickness = 0.3, & ! Minimum detectable optical thickness
59	CO2Slicing_PressureLimit = 700. * 100., & ! Cloud with higher pressures use thermal methods, units Pa
60	CO2Slicing_TauLimit = 1., & ! How deep into the cloud does CO2 slicing see?
61	phase_TauLimit = 1., & ! How deep into the cloud does the phase detection see?
62	size_TauLimit = 2., & ! Depth of the re retreivals
63	phaseDiscrimination_Threshold = 0.7 ! What fraction of total extincton needs to be
64	! in a single category to make phase discrim. work?
65	real, parameter :: re_fill= -999.
66	integer, parameter :: phaseIsNone = 0, phaseIsLiquid = 1, phaseIsIce = 2, phaseIsUndetermined = 3
67
68	logical, parameter :: useSimpleReScheme = .false.
69
70	! These are the limits of the libraries for the MODIS collection 5 algorithms
71	! They are also the limits used in the fits for g and w0
72
73	real, parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90.
74	integer, parameter :: num_trial_res = 15 ! increase to make the linear pseudo-retrieval of size more accurate
75	! DJS2015: Remove unused parameter
76	! logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations?
77	! DJS2015 END
78
79	! Precompute near-IR optical params vs size for retrieval scheme
80
81	integer, private :: i
82	real, dimension(num_trial_res), parameter :: &
83	trial_re_w = re_water_min + (re_water_max - re_water_min)/(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /), &
84	trial_re_i = re_ice_min + (re_ice_max - re_ice_min) /(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /)
85
86	! Can't initialze these during compilation, but do in before looping columns in retrievals
87	real, dimension(num_trial_res) :: g_w, g_i, w0_w, w0_i
88	! ------------------------------
89	! Bin boundaries for the joint optical thickness/cloud top pressure histogram
90
91	integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7
92
93	real, private :: dummy_real
94	real, dimension(numTauHistogramBins + 1), parameter :: &
95	tauHistogramBoundaries = (/ min_OpticalThickness, 1.3, 3.6, 9.4, 23., 60., 10000. /)
96	real, dimension(numPressureHistogramBins + 1), parameter :: & ! Units Pa
97	pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /)
98	real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100.
99
100	! For output - nominal bin centers and bin boundaries. On output pressure bins are highest to lowest.
101
102	integer, private :: k, l
103	real, parameter, dimension(2, numTauHistogramBins) :: &
104	nominalTauHistogramBoundaries = &
105	reshape(source = (/ tauHistogramBoundaries(1), &
106	((tauHistogramBoundaries(k), l = 1, 2), k = 2, numTauHistogramBins), &
107	100000. /), &
108	shape = (/2, numTauHistogramBins /) )
109	real, parameter, dimension(numTauHistogramBins) :: &
110	nominalTauHistogramCenters = (nominalTauHistogramBoundaries(1, :) + &
111	nominalTauHistogramBoundaries(2, :) ) / 2.
112
113	real, parameter, dimension(2, numPressureHistogramBins) :: &
114	nominalPressureHistogramBoundaries = &
115	reshape(source = (/ 100000., &
116	((pressureHistogramBoundaries(k), l = 1, 2), k = numPressureHistogramBins, 2, -1), &
117	0. /), &
118	shape = (/2, numPressureHistogramBins /) )
119	real, parameter, dimension(numPressureHistogramBins) :: &
120	nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + &
121	nominalPressureHistogramBoundaries(2, :) ) / 2.
122	! DJS2015 START: Add bin descriptions for joint-histograms of partice-sizes and optical depth. This is
123	! identical to what is done in COSPv.2.0.0 for histogram bin initialization.
124	integer :: j
125	integer,parameter :: &
126	numMODISReffLiqBins = 6, & ! Number of bins for tau/ReffLiq joint-histogram
127	numMODISReffIceBins = 6 ! Number of bins for tau/ReffICE joint-histogram
128	real,parameter,dimension(numMODISReffLiqBins+1) :: &
129	reffLIQ_binBounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/)
130	real,parameter,dimension(numMODISReffIceBins+1) :: &
131	reffICE_binBounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/)
132	real,parameter,dimension(2,numMODISReffIceBins) :: &
133	reffICE_binEdges = reshape(source=(/reffICE_binBounds(1),((reffICE_binBounds(k), &
134	l=1,2),k=2,numMODISReffIceBins),reffICE_binBounds(numMODISReffIceBins+1)/), &
135	shape = (/2,numMODISReffIceBins/))
136	real,parameter,dimension(2,numMODISReffLiqBins) :: &
137	reffLIQ_binEdges = reshape(source=(/reffLIQ_binBounds(1),((reffLIQ_binBounds(k), &
138	l=1,2),k=2,numMODISReffLiqBins),reffLIQ_binBounds(numMODISReffIceBins+1)/), &
139	shape = (/2,numMODISReffLiqBins/))
140	real,parameter,dimension(numMODISReffIceBins) :: &
141	reffICE_binCenters = (reffICE_binEdges(1,:)+reffICE_binEdges(2,:))/2.
142	real,parameter,dimension(numMODISReffLiqBins) :: &
143	reffLIQ_binCenters = (reffLIQ_binEdges(1,:)+reffLIQ_binEdges(2,:))/2.
144	! DJS2015 END
145
146	! ------------------------------
147	! There are two ways to call the MODIS simulator:
148	! 1) Provide total optical thickness and liquid/ice water content and we'll partition tau in
149	! subroutine modis_L2_simulator_oneTau, or
150	! 2) Provide ice and liquid optical depths in each layer
151
152	interface modis_L2_simulator
153	module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus
154	end interface
155	contains
156	!------------------------------------------------------------------------------------------------
157	! MODIS simulator using specified liquid and ice optical thickness in each layer
158
159	! Note: this simulator operates on all points; to match MODIS itself night-time
160	! points should be excluded
161
162	! Note: the simulator requires as input the optical thickness and cloud top pressure
163	! derived from the ISCCP simulator run with parameter top_height = 1.
164	! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing
165	! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that
166	! alogrithm in this simulator we simply report the values from the ISCCP simulator.
167
168	subroutine modis_L2_simulator_twoTaus( &
169	temp, pressureLayers, pressureLevels, &
170	liquid_opticalThickness, ice_opticalThickness, &
171	waterSize, iceSize, &
172	isccpTau, isccpCloudTopPressure, &
173	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
174
175	! Grid-mean quantities at layer centers, starting at the model top
176	! dimension nLayers
177	real, dimension(:), intent(in ) :: temp, & ! Temperature, K
178	pressureLayers, & ! Pressure, Pa
179	pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1)
180	! Sub-column quantities
181	! dimension nSubcols, nLayers
182	real, dimension(:, :), intent(in ) :: liquid_opticalThickness, & ! Layer optical thickness @ 0.67 microns due to liquid
183	ice_opticalThickness ! ditto, due to ice
184	real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns
185	iceSize ! Cloud ice effective radius, microns
186
187	! Cloud properties retrieved from ISCCP using top_height = 1
188	! dimension nSubcols
189	real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness
190	isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa)
191
192	! Properties retrieved by MODIS
193	! dimension nSubcols
194	integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer, defined in module header
195	real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
196	retrievedTau, & ! unitless
197	retrievedSize ! microns
198	! ---------------------------------------------------
199	! Local variables
200	logical, dimension(size(retrievedTau)) :: cloudMask
201	real, dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal
202	REAL :: integratedLiquidFraction
203	integer :: i, nSubcols, nLevels
204
205	! ---------------------------------------------------
206	nSubcols = size(liquid_opticalThickness, 1)
207	nLevels = size(liquid_opticalThickness, 2)
208
209	! Initial error checks
210
211	if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), &
212	size(isccpTau), size(isccpCloudTopPressure), &
213	size(retrievedPhase), size(retrievedCloudTopPressure), &
214	size(retrievedTau), size(retrievedSize) /) /= nSubcols )) &
215	call complain_and_die("Differing number of subcolumns in one or more arrays")
216
217	if(any((/ size(ice_opticalThickness, 2), size(waterSize, 2), size(iceSize, 2), &
218	size(temp), size(pressureLayers), size(pressureLevels)-1 /) /= nLevels )) &
219	call complain_and_die("Differing number of levels in one or more arrays")
220
221
222	if(any( (/ any(temp <= 0.), any(pressureLayers <= 0.), &
223	any(liquid_opticalThickness < 0.), &
224	any(ice_opticalThickness < 0.), &
225	any(waterSize < 0.), any(iceSize < 0.) /) )) &
226	call complain_and_die("Input values out of bounds")
227
228	! ---------------------------------------------------
229
230	! Compute the total optical thickness and the proportion due to liquid in each cell
231
232	where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.)
233	tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :))
234	elsewhere
235	tauLiquidFraction(:, :) = 0.
236	end where
237	tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :)
238
239	! Optical depth retrieval
240	! This is simply a sum over the optical thickness in each layer
241	! It should agree with the ISCCP values after min values have been excluded
242
243	retrievedTau(:) = sum(tauTotal(:, :), dim = 2)
244
245	! Cloud detection - does optical thickness exceed detection threshold?
246
247	cloudMask = retrievedTau(:) >= min_OpticalThickness
248
249	! Initialize initial estimates for size retrievals
250
251	if(any(cloudMask) .AND. .not. useSimpleReScheme) then
252	g_w(:) = get_g_nir( phaseIsLiquid, trial_re_w(:))
253	w0_w(:) = get_ssa_nir(phaseIsLiquid, trial_re_w(:))
254	g_i(:) = get_g_nir( phaseIsIce, trial_re_i(:))
255	w0_i(:) = get_ssa_nir(phaseIsIce, trial_re_i(:))
256	end if
257
258	DO i = 1, nSubCols
259	if(cloudMask(i)) then
260
261	! Cloud top pressure determination
262	! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds
263	! lower than that.
264	! For CO2 slicing we report the optical-depth weighted pressure, integrating to a specified
265	! optical depth
266	! This assumes linear variation in p between levels. Linear in ln(p) is probably better,
267	! though we'd need to deal with the lowest pressure gracefully.
268
269	retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), &
270	pressureLevels, &
271	CO2Slicing_TauLimit)
272
273	! Phase determination - determine fraction of total tau that's liquid
274	! When ice and water contribute about equally to the extinction we can't tell
275	! what the phase is
276
277	integratedLiquidFraction = weight_by_extinction(tauTotal(i, :), &
278	tauLiquidFraction(i, :), &
279	phase_TauLimit)
280	if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then
281	retrievedPhase(i) = phaseIsLiquid
282	else if (integratedLiquidFraction <= 1.- phaseDiscrimination_Threshold) then
283	retrievedPhase(i) = phaseIsIce
284	else
285	retrievedPhase(i) = phaseIsUndetermined
286	end if
287
288	! Size determination
289
290	if(useSimpleReScheme) then
291	! This is the extinction-weighted size considering only the phase we've chosen
292
293	if(retrievedPhase(i) == phaseIsIce) then
294	retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :), &
295	iceSize(i, :), &
296	(1. - integratedLiquidFraction) * size_TauLimit)
297
298	else if(retrievedPhase(i) == phaseIsLiquid) then
299	retrievedSize(i) = weight_by_extinction(liquid_opticalThickness(i, :), &
300	waterSize(i, :), &
301	integratedLiquidFraction * size_TauLimit)
302
303	else
304	retrievedSize(i) = 0.
305	end if
306	else
307	retrievedSize(i) = 1.0e-06*retrieve_re(retrievedPhase(i), retrievedTau(i), &
308	obs_Refl_nir = compute_nir_reflectance(liquid_opticalThickness(i, :), waterSize(i, :)*1.0e6, &
309	ice_opticalThickness(i, :), iceSize(i, :)*1.0e6))
310	end if
311	else
312
313	! Values when we don't think there's a cloud.
314
315	retrievedCloudTopPressure(i) = R_UNDEF
316	retrievedPhase(i) = phaseIsNone
317	retrievedSize(i) = R_UNDEF
318	retrievedTau(i) = R_UNDEF
319	end if
320	end do
321	where((retrievedSize(:) < 0.).AND.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill
322
323	! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS
324	! mimics what MODIS does to first order.
325	! Of course, ISCCP cloud top pressures are in mb.
326
327	where(cloudMask(:) .AND. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) &
328	retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100.
329
330	end subroutine modis_L2_simulator_twoTaus
331	!------------------------------------------------------------------------------------------------
332
333	! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents;
334	! we'll partition this into ice and liquid optical thickness and call the full MODIS simulator
335
336	subroutine modis_L2_simulator_oneTau( &
337	temp, pressureLayers, pressureLevels, &
338	opticalThickness, cloudWater, cloudIce, &
339	waterSize, iceSize, &
340	isccpTau, isccpCloudTopPressure, &
341	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
342	! Grid-mean quantities at layer centers,
343	! dimension nLayers
344	real, dimension(:), intent(in ) :: temp, & ! Temperature, K
345	pressureLayers, & ! Pressure, Pa
346	pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1)
347	! Sub-column quantities
348	! dimension nLayers, nSubcols
349	real, dimension(:, :), intent(in ) :: opticalThickness, & ! Layer optical thickness @ 0.67 microns
350	cloudWater, & ! Cloud water content, arbitrary units
351	cloudIce ! Cloud water content, same units as cloudWater
352	real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns
353	iceSize ! Cloud ice effective radius, microns
354
355	! Cloud properties retrieved from ISCCP using top_height = 1
356	! dimension nSubcols
357
358	real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness
359	isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa)
360
361	! Properties retrieved by MODIS
362	! dimension nSubcols
363	integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer
364	real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers
365	retrievedTau, & ! unitless
366	retrievedSize ! microns (or whatever units
367	! waterSize and iceSize are supplied in)
368	! ---------------------------------------------------
369	! Local variables
370	real, dimension(size(opticalThickness, 1), size(opticalThickness, 2)) :: &
371	liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction
372
373	REAL :: seuil
374	! ---------------------------------------------------
375
376	if (ok_debug_cosp) then
377	seuil=1.e-9
378	else
379	seuil=0.0
380	endif
381
382	where(cloudIce(:, :) <= 0.)
383	tauLiquidFraction(:, :) = 1.
384	elsewhere
385	where (cloudWater(:, :) <= 0.)
386	tauLiquidFraction(:, :) = 0.
387	elsewhere
388
389	! Geometic optics limit - tau as LWP/re (proportional to LWC/re)
390
391	! Modif AI 02 2018
392	tauLiquidFraction(:, :) = (cloudWater(:, :)/max(waterSize(:, :), seuil) ) / &
393	(cloudWater(:, :)/max(waterSize(:, :), seuil) + &
394	cloudIce(:, :)/(ice_density * max(iceSize(:, :), seuil)) )
395	end where
396	end where
397	liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :)
398	ice_opticalThickness (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :)
399
400	call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels, &
401	liquid_opticalThickness, ice_opticalThickness, &
402	waterSize, iceSize, &
403	isccpTau, isccpCloudTopPressure, &
404	retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
405
406	end subroutine modis_L2_simulator_oneTau
407
408	! ########################################################################################
409	subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, &
410	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
411	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
412	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
413	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
414	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, &
415	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
416	Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
417
418	! INPUTS
419	integer,intent(in) :: &
420	nPoints, & ! Number of horizontal gridpoints
421	nSubCols ! Number of subcolumns
422	integer,intent(in), dimension(:,:) :: &
423	!ds integer,intent(in), dimension(nPoints, nSubCols) :: &
424	phase
425	real,intent(in),dimension(:,:) :: &
426	!ds real,intent(in),dimension(nPoints, nSubCols) :: &
427	cloud_top_pressure, &
428	optical_thickness, &
429	particle_size
430
431	! OUTPUTS
432	real,intent(inout),dimension(:) :: & !
433	!ds real,intent(inout),dimension(nPoints) :: & !
434	Cloud_Fraction_Total_Mean, & !
435	Cloud_Fraction_Water_Mean, & !
436	Cloud_Fraction_Ice_Mean, & !
437	Cloud_Fraction_High_Mean, & !
438	Cloud_Fraction_Mid_Mean, & !
439	Cloud_Fraction_Low_Mean, & !
440	Optical_Thickness_Total_Mean, & !
441	Optical_Thickness_Water_Mean, & !
442	Optical_Thickness_Ice_Mean, & !
443	Optical_Thickness_Total_MeanLog10, & !
444	Optical_Thickness_Water_MeanLog10, & !
445	Optical_Thickness_Ice_MeanLog10, & !
446	Cloud_Particle_Size_Water_Mean, & !
447	Cloud_Particle_Size_Ice_Mean, & !
448	Cloud_Top_Pressure_Total_Mean, & !
449	Liquid_Water_Path_Mean, & !
450	Ice_Water_Path_Mean !
451	real,intent(inout),dimension(:,:,:) :: &
452	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: &
453	Optical_Thickness_vs_Cloud_Top_Pressure
454	real,intent(inout),dimension(:,:,:) :: &
455	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: &
456	Optical_Thickness_vs_ReffIce
457	real,intent(inout),dimension(:,:,:) :: &
458	!ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: &
459	Optical_Thickness_vs_ReffLiq
460
461	! LOCAL VARIABLES
462	real, parameter :: &
463	LWP_conversion = 2./3. * 1000. ! MKS units
464	integer :: i, j
465	logical, dimension(nPoints,nSubCols) :: &
466	cloudMask, &
467	waterCloudMask, &
468	iceCloudMask, &
469	validRetrievalMask
470	real,dimension(nPoints,nSubCols) :: &
471	tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
472
473	! ########################################################################################
474	! Include only those pixels with successful retrievals in the statistics
475	! ########################################################################################
476	validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
477	cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .AND. &
478	validRetrievalMask(1:nPoints,1:nSubCols)
479	waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .AND. &
480	validRetrievalMask(1:nPoints,1:nSubCols)
481	iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .AND. &
482	validRetrievalMask(1:nPoints,1:nSubCols)
483
484	! ########################################################################################
485	! Use these as pixel counts at first
486	! ########################################################################################
487	Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2))
488	Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
489	Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2))
490	Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .AND. cloud_top_pressure <= &
491	highCloudPressureLimit, dim = 2))
492	Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .AND. cloud_top_pressure > &
493	lowCloudPressureLimit, dim = 2))
494	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
495	- Cloud_Fraction_Low_Mean(1:nPoints)
496
497	! ########################################################################################
498	! Compute mean optical thickness.
499	! ########################################################################################
500	where (Cloud_Fraction_Total_Mean == 0)
501	Optical_Thickness_Total_Mean = R_UNDEF
502	elsewhere
503	Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / &
504	Cloud_Fraction_Total_Mean(1:nPoints)
505	endwhere
506	where (Cloud_Fraction_Water_Mean == 0)
507	Optical_Thickness_Water_Mean = R_UNDEF
508	elsewhere
509	Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
510	Cloud_Fraction_Water_Mean(1:nPoints)
511	endwhere
512	where (Cloud_Fraction_Ice_Mean == 0)
513	Optical_Thickness_Ice_Mean = R_UNDEF
514	elsewhere
515	Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / &
516	Cloud_Fraction_Ice_Mean(1:nPoints)
517	endwhere
518
519	! ########################################################################################
520	! We take the absolute value of optical thickness here to satisfy compilers that complains
521	! when we evaluate the logarithm of a negative number, even though it's not included in
522	! the sum.
523	! ########################################################################################
524	where (Cloud_Fraction_Total_Mean == 0)
525	Optical_Thickness_Total_MeanLog10 = R_UNDEF
526	elsewhere
527	Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
528	dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
529	endwhere
530	where (Cloud_Fraction_Water_Mean == 0)
531	Optical_Thickness_Water_MeanLog10 = R_UNDEF
532	elsewhere
533	Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
534	dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
535	endwhere
536	where (Cloud_Fraction_Ice_Mean == 0)
537	Optical_Thickness_Ice_MeanLog10 = R_UNDEF
538	elsewhere
539	Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
540	dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
541	endwhere
542	where (Cloud_Fraction_Water_Mean == 0)
543	Cloud_Particle_Size_Water_Mean = R_UNDEF
544	elsewhere
545	Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
546	Cloud_Fraction_Water_Mean(1:nPoints)
547	endwhere
548	where (Cloud_Fraction_Ice_Mean == 0)
549	Cloud_Particle_Size_Ice_Mean = R_UNDEF
550	elsewhere
551	Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / &
552	Cloud_Fraction_Ice_Mean(1:nPoints)
553	endwhere
554	where (Cloud_Fraction_Total_Mean == 0)
555	Cloud_Top_Pressure_Total_Mean = R_UNDEF
556	elsewhere
557	Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
558	max(1, count(cloudMask, dim = 2))
559	endwhere
560	where (Cloud_Fraction_Water_Mean == 0)
561	Liquid_Water_Path_Mean = R_UNDEF
562	elsewhere
563	Liquid_Water_Path_Mean(1:nPoints) = LWP_conversionsum(particle_sizeoptical_thickness, &
564	mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
565	endwhere
566	where (Cloud_Fraction_Ice_Mean == 0)
567	Ice_Water_Path_Mean = R_UNDEF
568	elsewhere
569	Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_densitysum(particle_sizeoptical_thickness,&
570	mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
571	endwhere
572	! ########################################################################################
573	! Normalize pixel counts to fraction.
574	! ########################################################################################
575	Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
576	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols
577	Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols
578	Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
579	Cloud_Fraction_Ice_Mean(1:nPoints) = Cloud_Fraction_Ice_Mean(1:nPoints) /nSubcols
580	Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
581
582	! ########################################################################################
583	! Set clear-scenes to undefined
584	! ########################################################################################
585	! where (Cloud_Fraction_High_Mean == 0) Cloud_Fraction_High_Mean = R_UNDEF
586	! where (Cloud_Fraction_Mid_Mean == 0) Cloud_Fraction_Mid_Mean = R_UNDEF
587	! where (Cloud_Fraction_Low_Mean == 0) Cloud_Fraction_Low_Mean = R_UNDEF
588
589	! ########################################################################################
590	! Joint histogram
591	! ########################################################################################
592
593	! Loop over all points
594	tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols)
595	ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols)
596	reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
597	reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
598	DO j=1,nPoints
599
600	! Fill clear and optically thin subcolumns with fill
601	where(.not. cloudMask(j,1:nSubCols))
602	tauWRK(j,1:nSubCols) = -999.
603	ctpWRK(j,1:nSubCols) = -999.
604	endwhere
605	! Joint histogram of tau/CTP
606	call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
607	tauHistogramBoundaries,numTauHistogramBins,&
608	pressureHistogramBoundaries,numPressureHistogramBins,&
609	Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins))
610	! Joint histogram of tau/ReffICE
611	call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, &
612	tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds, &
613	numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins))
614	! Joint histogram of tau/ReffLIQ
615	call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, &
616	tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds, &
617	numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins))
618
619	enddo
620	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = &
621	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols
622	Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = &
623	Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols
624	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = &
625	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols
626
627	end subroutine modis_column
628	! ######################################################################################
629	! SUBROUTINE hist2D
630	! ######################################################################################
631	subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
632	implicit none
633
634	! INPUTS
635	integer, intent(in) :: &
636	npts, & ! Number of data points to be sorted
637	nbin1, & ! Number of bins in histogram direction 1
638	nbin2 ! Number of bins in histogram direction 2
639	real,intent(in),dimension(npts) :: &
640	var1, & ! Variable 1 to be sorted into bins
641	var2 ! variable 2 to be sorted into bins
642	real,intent(in),dimension(nbin1+1) :: &
643	bin1 ! Histogram bin 1 boundaries
644	real,intent(in),dimension(nbin2+1) :: &
645	bin2 ! Histogram bin 2 boundaries
646	! OUTPUTS
647	real,intent(out),dimension(nbin1,nbin2) :: &
648	jointHist
649
650	! LOCAL VARIABLES
651	integer :: ij,ik
652
653	DO ij=2,nbin1+1
654	DO ik=2,nbin2+1
655	jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .AND. var1 .lt. bin1(ij) .AND. &
656	var2 .ge. bin2(ik-1) .AND. var2 .lt. bin2(ik))
657	enddo
658	enddo
659	end subroutine hist2D
660
661	!------------------------------------------------------------------------------------------------
662	subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, &
663	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
664	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
665	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
666	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
667	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
668	Cloud_Top_Pressure_Total_Mean, &
669	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
670	Optical_Thickness_vs_Cloud_Top_Pressure)
671
672	! Inputs; dimension nPoints, nSubcols
673
674	integer, dimension(:, :), intent(in) :: phase
675	real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size
676
677	! Outputs; dimension nPoints
678
679	real, dimension(:), intent(out) :: &
680	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
681	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
682	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
683	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
684	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, &
685	Cloud_Top_Pressure_Total_Mean, &
686	Liquid_Water_Path_Mean, Ice_Water_Path_Mean
687	! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins
688	real, dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure
689	! ---------------------------
690	! Local variables
691
692	real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units
693	integer :: i, j
694	integer :: nPoints, nSubcols
695	logical, dimension(size(phase, 1), size(phase, 2)) :: &
696	cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask
697	logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins ) :: tauMask
698	logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask
699	! ---------------------------
700
701	nPoints = size(phase, 1)
702	nSubcols = size(phase, 2)
703
704	! Array conformance checks
705
706	if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), &
707	size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), &
708	size(Cloud_Fraction_High_Mean), size(Cloud_Fraction_Mid_Mean), size(Cloud_Fraction_Low_Mean), &
709	size(Optical_Thickness_Total_Mean), size(Optical_Thickness_Water_Mean), size(Optical_Thickness_Ice_Mean), &
710	size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), &
711	size(Optical_Thickness_Ice_MeanLog10), size(Cloud_Particle_Size_Water_Mean), &
712	size(Cloud_Particle_Size_Ice_Mean), size(Cloud_Top_Pressure_Total_Mean), &
713	size(Liquid_Water_Path_Mean), size(Ice_Water_Path_Mean) /) /= nPoints)) &
714	call complain_and_die("Some L3 arrays have wrong number of grid points")
715	if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) &
716	call complain_and_die("Some L3 arrays have wrong number of subcolumns")
717
718	! Include only those pixels with successful retrievals in the statistics
719
720	validRetrievalMask(:, :) = particle_size(:, :) > 0.
721	cloudMask = phase(:, :) /= phaseIsNone .AND. validRetrievalMask(:, :)
722	waterCloudMask = phase(:, :) == phaseIsLiquid .AND. validRetrievalMask(:, :)
723	iceCloudMask = phase(:, :) == phaseIsIce .AND. validRetrievalMask(:, :)
724
725	! Use these as pixel counts at first
726
727	Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2))
728	Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2))
729	Cloud_Fraction_Ice_Mean(:) = real(count(iceCloudMask, dim = 2))
730
731	Cloud_Fraction_High_Mean(:) = real(count(cloudMask .AND. cloud_top_pressure <= highCloudPressureLimit, dim = 2))
732	Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .AND. cloud_top_pressure > lowCloudPressureLimit, dim = 2))
733	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:)
734
735	! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0.
736
737	where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1.
738	where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1.
739	where (Cloud_Fraction_Ice_Mean == 0) Cloud_Fraction_Ice_Mean = -1.
740
741	Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask, dim = 2) / Cloud_Fraction_Total_Mean(:)
742	Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
743	Optical_Thickness_Ice_Mean = sum(optical_thickness, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
744
745	! We take the absolute value of optical thickness here to satisfy compilers that complains when we
746	! evaluate the logarithm of a negative number, even though it's not included in the sum.
747	Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask, dim = 2) / &
748	Cloud_Fraction_Total_Mean(:)
749	Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / &
750	Cloud_Fraction_Water_Mean(:)
751	Optical_Thickness_Ice_MeanLog10 = sum(log10(abs(optical_thickness)), mask = iceCloudMask, dim = 2) / &
752	Cloud_Fraction_Ice_Mean(:)
753
754	Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
755	Cloud_Particle_Size_Ice_Mean = sum(particle_size, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:)
756
757	Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2))
758
759	Liquid_Water_Path_Mean = LWP_conversion &
760	* sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) &
761	/ Cloud_Fraction_Water_Mean(:)
762	Ice_Water_Path_Mean = LWP_conversion * ice_density &
763	* sum(particle_size * optical_thickness, mask = iceCloudMask, dim = 2) &
764	/ Cloud_Fraction_Ice_Mean(:)
765
766	! Normalize pixel counts to fraction
767	! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0.
768
769	Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols)
770	Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols)
771	Cloud_Fraction_Ice_Mean(:) = max(0., Cloud_Fraction_Ice_Mean(:) /nSubcols)
772
773	Cloud_Fraction_High_Mean(:) = Cloud_Fraction_High_Mean(:) /nSubcols
774	Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Mid_Mean(:) /nSubcols
775	Cloud_Fraction_Low_Mean(:) = Cloud_Fraction_Low_Mean(:) /nSubcols
776
777	! ----
778	! Joint histogram
779
780	DO i = 1, numTauHistogramBins
781	where(cloudMask(:, :))
782	tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .AND. &
783	optical_thickness(:, :) < tauHistogramBoundaries(i+1)
784	elsewhere
785	tauMask(:, :, i) = .false.
786	end where
787	end do
788
789	DO i = 1, numPressureHistogramBins
790	where(cloudMask(:, :))
791	pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .AND. &
792	cloud_top_pressure(:, :) < pressureHistogramBoundaries(i+1)
793	elsewhere
794	pressureMask(:, :, i) = .false.
795	end where
796	end do
797
798	DO i = 1, numPressureHistogramBins
799	DO j = 1, numTauHistogramBins
800	Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = &
801	real(count(tauMask(:, :, j) .AND. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
802	end do
803	end do
804
805	end subroutine modis_L3_simulator
806	!------------------------------------------------------------------------------------------------
807	function cloud_top_pressure(tauIncrement, pressure, tauLimit)
808	real, dimension(:), intent(in) :: tauIncrement, pressure
809	real, intent(in) :: tauLimit
810	REAL :: cloud_top_pressure
811
812	! Find the extinction-weighted pressure. Assume that pressure varies linearly between
813	! layers and use the trapezoidal rule.
814
815	REAL :: deltaX, totalTau, totalProduct
816	integer :: i
817
818	totalTau = 0.; totalProduct = 0.
819	DO i = 2, size(tauIncrement)
820	if(totalTau + tauIncrement(i) > tauLimit) then
821	deltaX = tauLimit - totalTau
822	totalTau = totalTau + deltaX
823
824	! Result for trapezoidal rule when you take less than a full step
825	! tauIncrement is a layer-integrated value
826
827	totalProduct = totalProduct &
828	+ pressure(i-1) * deltaX &
829	+ (pressure(i) - pressure(i-1)) * deltaX*2/(2. tauIncrement(i))
830	else
831	totalTau = totalTau + tauIncrement(i)
832	totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2.
833	end if
834	if(totalTau >= tauLimit) exit
835	end do
836	cloud_top_pressure = totalProduct/totalTau
837	end function cloud_top_pressure
838	!------------------------------------------------------------------------------------------------
839	function weight_by_extinction(tauIncrement, f, tauLimit)
840	real, dimension(:), intent(in) :: tauIncrement, f
841	real, intent(in) :: tauLimit
842	REAL :: weight_by_extinction
843
844	! Find the extinction-weighted value of f(tau), assuming constant f within each layer
845
846	REAL :: deltaX, totalTau, totalProduct
847	integer :: i
848
849	totalTau = 0.; totalProduct = 0.
850	DO i = 1, size(tauIncrement)
851	if(totalTau + tauIncrement(i) > tauLimit) then
852	deltaX = tauLimit - totalTau
853	totalTau = totalTau + deltaX
854	totalProduct = totalProduct + deltaX * f(i)
855	else
856	totalTau = totalTau + tauIncrement(i)
857	totalProduct = totalProduct + tauIncrement(i) * f(i)
858	end if
859	if(totalTau >= tauLimit) exit
860	end do
861	weight_by_extinction = totalProduct/totalTau
862	end function weight_by_extinction
863	!------------------------------------------------------------------------------------------------
864	pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size)
865	real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
866	REAL :: compute_nir_reflectance
867
868	real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
869	tau, g, w0
870	!----------------------------------------
871	water_g(:) = get_g_nir( phaseIsLiquid, water_size)
872	water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size)
873	ice_g(:) = get_g_nir( phaseIsIce, ice_size)
874	ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size)
875
876	! Combine ice and water optical properties
877
878	g(:) = 0; w0(:) = 0.
879	tau(:) = ice_tau(:) + water_tau(:)
880	where (tau(:) > 0)
881	w0(:) = (water_tau(:) * water_w0(:) + ice_tau(:) * ice_w0(:)) / &
882	tau(:)
883	g(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
884	(w0(:) * tau(:))
885	end where
886
887	compute_nir_reflectance = compute_toa_reflectace(tau, g, w0)
888	end function compute_nir_reflectance
889	!------------------------------------------------------------------------------------------------
890	! Retreivals
891	!------------------------------------------------------------------------------------------------
892	elemental function retrieve_re (phase, tau, obs_Refl_nir)
893	integer, intent(in) :: phase
894	real, intent(in) :: tau, obs_Refl_nir
895	REAL :: retrieve_re
896
897	! Finds the re that produces the minimum mis-match between predicted and observed reflectance in
898	! MODIS band 7 (near IR)
899	! Uses
900	! fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables
901	! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0
902	! adding-doubling for total reflectance
903
904
905
906	! Local variables
907
908	real, parameter :: min_distance_to_boundary = 0.01
909	REAL :: re_min, re_max, delta_re
910	integer :: i
911
912	real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir
913	! --------------------------
914
915	if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then
916	if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then
917	re_min = re_water_min
918	re_max = re_water_max
919	trial_re(:) = trial_re_w
920	g(:) = g_w(:)
921	w0(:) = w0_w(:)
922	else
923	re_min = re_ice_min
924	re_max = re_ice_max
925	trial_re(:) = trial_re_i
926	g(:) = g_i(:)
927	w0(:) = w0_i(:)
928	end if
929
930	! 1st attempt at index: w/coarse re resolution
931
932	predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
933	retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
934
935	! If first retrieval works, can try 2nd iteration using greater re resolution
936
937	! DJS2015: Remove unused piece of code
938	! if(use_two_re_iterations .AND. retrieve_re > 0.) then
939	! re_min = retrieve_re - delta_re
940	! re_max = retrieve_re + delta_re
941	! delta_re = (re_max - re_min)/real(num_trial_res-1)
942
943	! trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
944	! g(:) = get_g_nir( phase, trial_re(:))
945	! w0(:) = get_ssa_nir(phase, trial_re(:))
946	! predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
947	! retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
948	! end if
949	! DJS2015 END
950	else
951	retrieve_re = re_fill
952	end if
953
954	end function retrieve_re
955	! --------------------------------------------
956	pure function interpolate_to_min(x, y, yobs)
957	real, dimension(:), intent(in) :: x, y
958	real, intent(in) :: yobs
959	REAL :: interpolate_to_min
960
961	! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
962	! y must be monotonic in x
963
964	real, dimension(size(x)) :: diff
965	integer :: nPoints, minDiffLoc, lowerBound, upperBound
966	! ---------------------------------
967	nPoints = size(y)
968	diff(:) = y(:) - yobs
969	minDiffLoc = minloc(abs(diff), dim = 1)
970
971	if(minDiffLoc == 1) then
972	lowerBound = minDiffLoc
973	upperBound = minDiffLoc + 1
974	else if(minDiffLoc == nPoints) then
975	lowerBound = minDiffLoc - 1
976	upperBound = minDiffLoc
977	else
978	if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
979	lowerBound = minDiffLoc-1
980	upperBound = minDiffLoc
981	else
982	lowerBound = minDiffLoc
983	upperBound = minDiffLoc + 1
984	end if
985	end if
986
987	if(diff(lowerBound) * diff(upperBound) < 0) then
988
989	! Interpolate the root position linearly if we bracket the root
990
991	interpolate_to_min = x(upperBound) - &
992	diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
993	else
994	interpolate_to_min = re_fill
995	end if
996
997
998	end function interpolate_to_min
999	! --------------------------------------------
1000	! Optical properties
1001	! --------------------------------------------
1002	elemental function get_g_nir (phase, re)
1003
1004	! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
1005	! of size for ice and water
1006	! Fits from Steve Platnick
1007
1008	integer, intent(in) :: phase
1009	real, intent(in) :: re
1010	REAL :: get_g_nir
1011
1012	real, dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
1013	small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
1014	real, dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
1015
1016	! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
1017	if(phase == phaseIsLiquid) then
1018	if(re < 7.) then
1019	get_g_nir = fit_to_quadratic(re, small_water_coefficients)
1020	if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
1021	else
1022	get_g_nir = fit_to_cubic(re, big_water_coefficients)
1023	if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
1024	end if
1025	else
1026	get_g_nir = fit_to_quadratic(re, ice_coefficients)
1027	if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
1028	if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
1029	end if
1030
1031	end function get_g_nir
1032
1033	! --------------------------------------------
1034	elemental function get_ssa_nir (phase, re)
1035	integer, intent(in) :: phase
1036	real, intent(in) :: re
1037	REAL :: get_ssa_nir
1038
1039	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
1040	! of size for ice and water
1041	! Fits from Steve Platnick
1042
1043	real, dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
1044	real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
1045
1046	! approx. fits from MODIS Collection 6 LUT scattering calculations
1047	if(phase == phaseIsLiquid) then
1048	get_ssa_nir = fit_to_quadratic(re, water_coefficients)
1049	if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
1050	if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
1051	else
1052	get_ssa_nir = fit_to_cubic(re, ice_coefficients)
1053	if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
1054	if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
1055	end if
1056
1057	end function get_ssa_nir
1058	! --------------------------------------------
1059	pure function fit_to_cubic(x, coefficients)
1060	real, intent(in) :: x
1061	real, dimension(:), intent(in) :: coefficients
1062	REAL :: fit_to_cubic
1063
1064
1065	fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
1066	end function fit_to_cubic
1067	! --------------------------------------------
1068	pure function fit_to_quadratic(x, coefficients)
1069	real, intent(in) :: x
1070	real, dimension(:), intent(in) :: coefficients
1071	REAL :: fit_to_quadratic
1072
1073
1074	fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
1075	end function fit_to_quadratic
1076	! --------------------------------------------
1077	! Radiative transfer
1078	! --------------------------------------------
1079	pure function compute_toa_reflectace(tau, g, w0)
1080	real, dimension(:), intent(in) :: tau, g, w0
1081	REAL :: compute_toa_reflectace
1082
1083	logical, dimension(size(tau)) :: cloudMask
1084	integer, dimension(count(tau(:) > 0)) :: cloudIndicies
1085	real, dimension(count(tau(:) > 0)) :: Refl, Trans
1086	REAL :: Refl_tot, Trans_tot
1087	integer :: i
1088	! ---------------------------------------
1089
1090	! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation
1091
1092	cloudMask = tau(:) > 0.
1093	cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask)
1094	DO i = 1, size(cloudIndicies)
1095	call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
1096	end do
1097
1098	call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot)
1099
1100	compute_toa_reflectace = Refl_tot
1101
1102	end function compute_toa_reflectace
1103	! --------------------------------------------
1104	pure subroutine two_stream(tauint, gint, w0int, ref, tra)
1105	real, intent(in) :: tauint, gint, w0int
1106	real, intent(out) :: ref, tra
1107
1108	! Compute reflectance in a single layer using the two stream approximation
1109	! The code itself is from Lazaros Oreopoulos via Steve Platnick
1110
1111	! ------------------------
1112	! Local variables
1113	! for delta Eddington code
1114	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
1115	integer, parameter :: beam = 2
1116	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
1117	REAL :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
1118	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
1119
1120	! Compute reflectance and transmittance in a single layer using the two stream approximation
1121	! The code itself is from Lazaros Oreopoulos via Steve Platnick
1122
1123	f = gint**2
1124	tau = (1 - w0int * f) * tauint
1125	w0 = (1 - f) * w0int / (1 - w0int * f)
1126	g = (gint - f) / (1 - f)
1127
1128	! delta-Eddington (Joseph et al. 1976)
1129	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
1130	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
1131	gamma3 = (2 - 3gxmu) / 4.0
1132	gamma4 = 1 - gamma3
1133
1134	if (w0int > minConservativeW0) then
1135	! Conservative scattering
1136	if (beam == 1) then
1137	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
1138	ref = rh / (1 + gamma1 * tau)
1139	tra = 1 - ref
1140	else if(beam == 2) then
1141	ref = gamma1tau/(1 + gamma1tau)
1142	tra = 1 - ref
1143	endif
1144	else
1145	! Non-conservative scattering
1146	a1 = gamma1 * gamma4 + gamma2 * gamma3
1147	a2 = gamma1 * gamma3 + gamma2 * gamma4
1148
1149	rk = sqrt(gamma12 - gamma22)
1150
1151	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
1152	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
1153	r3 = 2 * rk (gamma3 - a2 xmu)
1154	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
1155	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
1156
1157	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
1158	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
1159	t3 = 2 * rk * (gamma4 + a1 * xmu)
1160	t4 = r4
1161	t5 = r5
1162
1163	beta = -r5 / r4
1164
1165	e1 = min(rk * tau, 500.)
1166	e2 = min(tau / xmu, 500.)
1167
1168	if (beam == 1) then
1169	den = r4 * exp(e1) + r5 * exp(-e1)
1170	ref = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
1171	den = t4 * exp(e1) + t5 * exp(-e1)
1172	th = exp(-e2)
1173	tra = th-thw0(t1exp(e1)-t2exp(-e1)-t3*exp(e2))/den
1174	elseif (beam == 2) then
1175	ef1 = exp(-e1)
1176	ef2 = exp(-2*e1)
1177	ref = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
1178	tra = (2rkef1)/((rk+gamma1)(1-betaef2))
1179	endif
1180	end if
1181	end subroutine two_stream
1182	! --------------------------------------------------
1183	elemental function two_stream_reflectance(tauint, gint, w0int)
1184	real, intent(in) :: tauint, gint, w0int
1185	REAL :: two_stream_reflectance
1186
1187	! Compute reflectance in a single layer using the two stream approximation
1188	! The code itself is from Lazaros Oreopoulos via Steve Platnick
1189
1190	! ------------------------
1191	! Local variables
1192	! for delta Eddington code
1193	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
1194	integer, parameter :: beam = 2
1195	real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
1196	REAL :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
1197	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
1198	! ------------------------
1199
1200
1201	f = gint**2
1202	tau = (1 - w0int * f) * tauint
1203	w0 = (1 - f) * w0int / (1 - w0int * f)
1204	g = (gint - f) / (1 - f)
1205
1206	! delta-Eddington (Joseph et al. 1976)
1207	gamma1 = (7 - w0* (4 + 3 * g)) / 4.0
1208	gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
1209	gamma3 = (2 - 3gxmu) / 4.0
1210	gamma4 = 1 - gamma3
1211
1212	if (w0int > minConservativeW0) then
1213	! Conservative scattering
1214	if (beam == 1) then
1215	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
1216	two_stream_reflectance = rh / (1 + gamma1 * tau)
1217	elseif (beam == 2) then
1218	two_stream_reflectance = gamma1tau/(1 + gamma1tau)
1219	endif
1220
1221	else !
1222
1223	! Non-conservative scattering
1224	a1 = gamma1 * gamma4 + gamma2 * gamma3
1225	a2 = gamma1 * gamma3 + gamma2 * gamma4
1226
1227	rk = sqrt(gamma12 - gamma22)
1228
1229	r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
1230	r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
1231	r3 = 2 * rk (gamma3 - a2 xmu)
1232	r4 = (1 - (rk * xmu)*2) (rk + gamma1)
1233	r5 = (1 - (rk * xmu)*2) (rk - gamma1)
1234
1235	t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
1236	t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
1237	t3 = 2 * rk * (gamma4 + a1 * xmu)
1238	t4 = r4
1239	t5 = r5
1240
1241	beta = -r5 / r4
1242
1243	e1 = min(rk * tau, 500.)
1244	e2 = min(tau / xmu, 500.)
1245
1246	if (beam == 1) then
1247	den = r4 * exp(e1) + r5 * exp(-e1)
1248	two_stream_reflectance = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
1249	elseif (beam == 2) then
1250	ef1 = exp(-e1)
1251	ef2 = exp(-2*e1)
1252	two_stream_reflectance = (gamma2(1-ef2))/((rk+gamma1)(1-beta*ef2))
1253	endif
1254
1255	end if
1256	end function two_stream_reflectance
1257	! --------------------------------------------
1258	pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot)
1259	real, dimension(:), intent(in) :: Refl, Tran
1260	real, intent(out) :: Refl_tot, Tran_tot
1261
1262	! Use adding/doubling formulas to compute total reflectance and transmittance from layer values
1263
1264	integer :: i
1265	real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative
1266
1267	Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)
1268
1269	DO i=2, size(Refl)
1270	! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
1271	Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)(Tran_cumulative(i-1)2)/(1 - Refl_cumulative(i-1) Refl(i))
1272	Tran_cumulative(i) = (Tran_cumulative(i-1)Tran(i)) / (1 - Refl_cumulative(i-1) Refl(i))
1273	end do
1274
1275	Refl_tot = Refl_cumulative(size(Refl))
1276	Tran_tot = Tran_cumulative(size(Refl))
1277
1278	end subroutine adding_doubling
1279	! --------------------------------------------------
1280	subroutine complain_and_die(message)
1281	character(len = *), intent(in) :: message
1282
1283	write(6, *) "Failure in MODIS simulator"
1284	write(6, *) trim(message)
1285	stop
1286	end subroutine complain_and_die
1287	!------------------------------------------------------------------------------------------------
1288	end module mod_modis_sim

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

Download in other formats:

Original Format