Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

modis_simulator.F90 @ 3702

Last change on this file since 3702 was 3491, checked in by idelkadi, 6 years ago

Integration of version 2 of the COSP simulator in LMDZ
This line, and those below, will be ignored--

M makegcm
M makelmdz
M makelmdz_fcm
M libf/phylmd/physiq_mod.F90
A libf/phylmd/cospv2
A libf/phylmd/cospv2/mo_rng.F90
A libf/phylmd/cospv2/quickbeam_optics.F90
A libf/phylmd/cospv2/cosp_cloudsat_interface.F90
A libf/phylmd/cospv2/cosp_config.F90
A libf/phylmd/cospv2/lidar_simulator.F90
A libf/phylmd/cospv2/prec_scops.F90
A libf/phylmd/cospv2/mrgrnk.F90
A libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90
A libf/phylmd/cospv2/cosp_atlid_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90
A libf/phylmd/cospv2/cosp_math_constants.F90
A libf/phylmd/cospv2/MISR_simulator.F90
A libf/phylmd/cospv2/modis_simulator.F90
A libf/phylmd/cospv2/math_lib.F90
A libf/phylmd/cospv2/cosp_grLidar532_interface.F90
A libf/phylmd/cospv2/cosp_errorHandling.F90
A libf/phylmd/cospv2/cosp_stats.F90
A libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90
A libf/phylmd/cospv2/cosp_utils.F90
A libf/phylmd/cospv2/cosp_optics.F90
A libf/phylmd/cospv2/icarus.F90
A libf/phylmd/cospv2/scops.F90
A libf/phylmd/cospv2/optics_lib.F90
A libf/phylmd/cospv2/cosp_kinds.F90
A libf/phylmd/cospv2/cosp_calipso_interface.F90
A libf/phylmd/cospv2/quickbeam.F90
A libf/phylmd/cospv2/parasol.F90
A libf/phylmd/cospv2/cosp_phys_constants.F90
A libf/phylmd/cospv2/cosp.F90
A libf/phylmd/cospv2/array_lib.F90
A libf/phylmd/cospv2/cosp_isccp_interface.F90
A libf/phylmd/cospv2/cosp_parasol_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_output_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_interface.F90
A libf/phylmd/cospv2/cosp_misr_interface.F90
A libf/phylmd/cospv2/cosp_modis_interface.F90

File size: 44.2 KB

Line
1	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2	! Copyright (c) 2015, Regents of the University of Colorado
3	! All rights reserved.
4	!
5	! Redistribution and use in source and binary forms, with or without modification, are
6	! permitted provided that the following conditions are met:
7	!
8	! 1. Redistributions of source code must retain the above copyright notice, this list of
9	! conditions and the following disclaimer.
10	!
11	! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12	! of conditions and the following disclaimer in the documentation and/or other
13	! materials provided with the distribution.
14	!
15	! 3. Neither the name of the copyright holder nor the names of its contributors may be
16	! used to endorse or promote products derived from this software without specific prior
17	! written permission.
18	!
19	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20	! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21	! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22	! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24	! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26	! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28	!
29	! History
30	! May 2009: Robert Pincus - Initial version
31	! June 2009: Steve Platnick and Robert Pincus - Simple radiative transfer for size
32	! 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
35	! using AM2 (GFDL)
36	! January 2010: Robert Pincus - Added high, middle, low cloud fractions
37	! May 2015: Dustin Swales - Modified for COSPv2.0
38	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39	!
40	! Notes on using the MODIS simulator:
41	! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1
42	! microns, or optical thickness at 0.67 microns and ice- and liquid-water contents
43	! (in consistent units of your choosing)
44	! *) Required input also includes the optical thickness and cloud top pressure
45	! derived from the ISCCP simulator run with parameter top_height = 1.
46	! *) Cloud particle sizes are specified as radii, measured in meters, though within the
47	! module we use units of microns. Where particle sizes are outside the bounds used in
48	! the MODIS retrieval libraries (parameters re_water_min, re_ice_min, etc.) the
49	! simulator returns missing values (re_fill)
50	!
51	! When error conditions are encountered this code calls the function complain_and_die,
52	! supplied at the bottom of this module. Users probably want to replace this with
53	! something more graceful.
54	!
55	module mod_modis_sim
56	USE MOD_COSP_CONFIG, only: R_UNDEF,modis_histTau,modis_histPres,numMODISTauBins, &
57	numMODISPresBins,numMODISReffIceBins,numMODISReffLiqBins, &
58	modis_histReffIce,modis_histReffLiq
59	USE COSP_KINDS, ONLY: wp
60	use MOD_COSP_STATS, ONLY: hist2D
61
62	implicit none
63	! ##########################################################################
64	! Retrieval parameters
65	integer, parameter :: &
66	num_trial_res = 15 ! Increase to make the linear pseudo-retrieval of size more accurate
67
68	real(wp) :: &
69	min_OpticalThickness, & ! Minimum detectable optical thickness
70	CO2Slicing_PressureLimit, & ! Cloud with higher pressures use thermal methods, units Pa
71	CO2Slicing_TauLimit, & ! How deep into the cloud does CO2 slicing see?
72	phase_TauLimit, & ! How deep into the cloud does the phase detection see?
73	size_TauLimit, & ! Depth of the re retreivals
74	phaseDiscrimination_Threshold, & ! What fraction of total extincton needs to be in a single
75	! category to make phase discrim. work?
76	re_fill, & !
77	re_water_min, & ! Minimum effective radius (liquid)
78	re_water_max, & ! Maximum effective radius (liquid)
79	re_ice_min, & ! Minimum effective radius (ice)
80	re_ice_max, & ! Minimum effective radius (ice)
81	highCloudPressureLimit, & ! High cloud pressure limit (Pa)
82	lowCloudPressureLimit ! Low cloud pressure limit (Pa)
83	integer :: &
84	phaseIsNone, & !
85	phaseIsLiquid, & !
86	phaseIsIce, & !
87	phaseIsUndetermined !
88
89	real(wp),dimension(num_trial_res) :: &
90	trial_re_w, & ! Near-IR optical params vs size for retrieval scheme (liquid)
91	trial_re_i ! Near-IR optical params vs size for retrieval scheme (ice)
92	real(wp),dimension(num_trial_res) :: &
93	g_w, & ! Assymettry parameter for size retrieval (liquid)
94	g_i, & ! Assymettry parameter for size retrieval (ice)
95	w0_w, & ! Single-scattering albedo for size retrieval (liquid)
96	w0_i ! Single-scattering albedo for size retrieval (ice)
97	! Algorithmic parameters
98	real(wp),parameter :: &
99	ice_density = 0.93_wp ! Liquid density is 1.
100
101	contains
102	! ########################################################################################
103	! MODIS simulator using specified liquid and ice optical thickness in each layer
104	!
105	! Note: this simulator operates on all points; to match MODIS itself night-time
106	! points should be excluded
107	!
108	! Note: the simulator requires as input the optical thickness and cloud top pressure
109	! derived from the ISCCP simulator run with parameter top_height = 1.
110	! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing
111	! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that
112	! alogrithm in this simulator we simply report the values from the ISCCP simulator.
113	! ########################################################################################
114	subroutine modis_subcolumn(nSubCols, nLevels, pressureLevels, optical_thickness, &
115	tauLiquidFraction, g, w0,isccpCloudTopPressure, &
116	retrievedPhase, retrievedCloudTopPressure, &
117	retrievedTau, retrievedSize)
118
119	! INPUTS
120	integer,intent(in) :: &
121	nSubCols, & ! Number of subcolumns
122	nLevels ! Number of levels
123	real(wp),dimension(nLevels+1),intent(in) :: &
124	pressureLevels ! Gridmean pressure at layer edges (Pa)
125	real(wp),dimension(nSubCols,nLevels),intent(in) :: &
126	optical_thickness, & ! Subcolumn optical thickness @ 0.67 microns.
127	tauLiquidFraction, & ! Liquid water fraction
128	g, & ! Subcolumn assymetry parameter
129	w0 ! Subcolumn single-scattering albedo
130	real(wp),dimension(nSubCols),intent(in) :: &
131	isccpCloudTopPressure ! ISCCP retrieved cloud top pressure (Pa)
132
133	! OUTPUTS
134	integer, dimension(nSubCols), intent(inout) :: &
135	retrievedPhase ! MODIS retrieved phase (liquid/ice/other)
136	real(wp),dimension(nSubCols), intent(inout) :: &
137	retrievedCloudTopPressure, & ! MODIS retrieved CTP (Pa)
138	retrievedTau, & ! MODIS retrieved optical depth (unitless)
139	retrievedSize ! MODIS retrieved particle size (microns)
140
141	! LOCAL VARIABLES
142	logical, dimension(nSubCols) :: &
143	cloudMask
144	real(wp) :: &
145	integratedLiquidFraction, &
146	obs_Refl_nir
147	real(wp),dimension(num_trial_res) :: &
148	predicted_Refl_nir
149	integer :: &
150	i
151
152	! ########################################################################################
153	! Optical depth retrieval
154	! This is simply a sum over the optical thickness in each layer.
155	! It should agree with the ISCCP values after min values have been excluded.
156	! ########################################################################################
157	retrievedTau(1:nSubCols) = sum(optical_thickness(1:nSubCols,1:nLevels), dim = 2)
158
159	! ########################################################################################
160	! Cloud detection
161	! does optical thickness exceed detection threshold?
162	! ########################################################################################
163	cloudMask = retrievedTau(1:nSubCols) >= min_OpticalThickness
164
165	do i = 1, nSubCols
166	if(cloudMask(i)) then
167	! ##################################################################################
168	! Cloud top pressure determination
169	! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal
170	! methods for clouds lower than that. For CO2 slicing we report the optical-depth
171	! weighted pressure, integrating to a specified optical depth.
172	! This assumes linear variation in p between levels. Linear in ln(p) is probably
173	! better, though we'd need to deal with the lowest pressure gracefully.
174	! ##################################################################################
175	retrievedCloudTopPressure(i) = cloud_top_pressure(nLevels,(/ 0._wp, optical_thickness(i,1:nLevels) /), &
176	pressureLevels(1:nLevels),CO2Slicing_TauLimit)
177
178	! ##################################################################################
179	! Phase determination
180	! Determine fraction of total tau that's liquid when ice and water contribute about
181	! equally to the extinction we can't tell what the phase is.
182	! ##################################################################################
183	integratedLiquidFraction = weight_by_extinction(nLevels,optical_thickness(i,1:nLevels), &
184	tauLiquidFraction(i, 1:nLevels), &
185	phase_TauLimit)
186	if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then
187	retrievedPhase(i) = phaseIsLiquid
188	else if (integratedLiquidFraction <= 1._wp- phaseDiscrimination_Threshold) then
189	retrievedPhase(i) = phaseIsIce
190	else
191	retrievedPhase(i) = phaseIsUndetermined
192	end if
193
194	! ##################################################################################
195	! Size determination
196	! ##################################################################################
197
198	! Compute observed reflectance
199	obs_Refl_nir = compute_toa_reflectace(nLevels,optical_thickness(i,1:nLevels), g(i,1:nLevels), w0(i,1:nLevels))
200
201	! Compute predicted reflectance
202	if(any(retrievedPhase(i) == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then
203	if (retrievedPhase(i) == phaseIsLiquid .OR. retrievedPhase(i) == phaseIsUndetermined) then
204	predicted_Refl_nir(1:num_trial_res) = two_stream_reflectance(retrievedTau(i), &
205	g_w(1:num_trial_res), w0_w(1:num_trial_res))
206	retrievedSize(i) = 1.0e-06_wp*interpolate_to_min(trial_re_w(1:num_trial_res), &
207	predicted_Refl_nir(1:num_trial_res), obs_Refl_nir)
208	else
209	predicted_Refl_nir(1:num_trial_res) = two_stream_reflectance(retrievedTau(i), &
210	g_i(1:num_trial_res), w0_i(1:num_trial_res))
211	retrievedSize(i) = 1.0e-06_wp*interpolate_to_min(trial_re_i(1:num_trial_res), &
212	predicted_Refl_nir(1:num_trial_res), obs_Refl_nir)
213	endif
214	else
215	retrievedSize(i) = re_fill
216	endif
217	else
218	! Values when we don't think there's a cloud.
219	retrievedCloudTopPressure(i) = R_UNDEF
220	retrievedPhase(i) = phaseIsNone
221	retrievedSize(i) = R_UNDEF
222	retrievedTau(i) = R_UNDEF
223	end if
224	end do
225	where((retrievedSize(1:nSubCols) < 0.).and.(retrievedSize(1:nSubCols) /= R_UNDEF)) &
226	retrievedSize(1:nSubCols) = 1.0e-06_wp*re_fill
227
228	! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS
229	! mimics what MODIS does to first order.
230	! Of course, ISCCP cloud top pressures are in mb.
231	where(cloudMask(1:nSubCols) .and. retrievedCloudTopPressure(1:nSubCols) > CO2Slicing_PressureLimit) &
232	retrievedCloudTopPressure(1:nSubCols) = isccpCloudTopPressure! * 100._wp
233
234	end subroutine modis_subcolumn
235
236	! ########################################################################################
237	subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, &
238	Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, &
239	Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, &
240	Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, &
241	Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
242	Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, &
243	Liquid_Water_Path_Mean, Ice_Water_Path_Mean, &
244	Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
245
246	! INPUTS
247	integer,intent(in) :: &
248	nPoints, & ! Number of horizontal gridpoints
249	nSubCols ! Number of subcolumns
250	integer,intent(in), dimension(nPoints, nSubCols) :: &
251	phase
252	real(wp),intent(in),dimension(nPoints, nSubCols) :: &
253	cloud_top_pressure, &
254	optical_thickness, &
255	particle_size
256
257	! OUTPUTS
258	real(wp),intent(inout),dimension(nPoints) :: & !
259	Cloud_Fraction_Total_Mean, & !
260	Cloud_Fraction_Water_Mean, & !
261	Cloud_Fraction_Ice_Mean, & !
262	Cloud_Fraction_High_Mean, & !
263	Cloud_Fraction_Mid_Mean, & !
264	Cloud_Fraction_Low_Mean, & !
265	Optical_Thickness_Total_Mean, & !
266	Optical_Thickness_Water_Mean, & !
267	Optical_Thickness_Ice_Mean, & !
268	Optical_Thickness_Total_MeanLog10, & !
269	Optical_Thickness_Water_MeanLog10, & !
270	Optical_Thickness_Ice_MeanLog10, & !
271	Cloud_Particle_Size_Water_Mean, & !
272	Cloud_Particle_Size_Ice_Mean, & !
273	Cloud_Top_Pressure_Total_Mean, & !
274	Liquid_Water_Path_Mean, & !
275	Ice_Water_Path_Mean !
276	real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISPresBins) :: &
277	Optical_Thickness_vs_Cloud_Top_Pressure
278	real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISReffIceBins) :: &
279	Optical_Thickness_vs_ReffIce
280	real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISReffLiqBins) :: &
281	Optical_Thickness_vs_ReffLiq
282
283	! LOCAL VARIABLES
284	real(wp), parameter :: &
285	LWP_conversion = 2._wp/3._wp * 1000._wp ! MKS units
286	integer :: j
287	logical, dimension(nPoints,nSubCols) :: &
288	cloudMask, &
289	waterCloudMask, &
290	iceCloudMask, &
291	validRetrievalMask
292	real(wp),dimension(nPoints,nSubCols) :: &
293	tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
294
295	! ########################################################################################
296	! Include only those pixels with successful retrievals in the statistics
297	! ########################################################################################
298	validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
299	cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and. &
300	validRetrievalMask(1:nPoints,1:nSubCols)
301	waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. &
302	validRetrievalMask(1:nPoints,1:nSubCols)
303	iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .and. &
304	validRetrievalMask(1:nPoints,1:nSubCols)
305
306	! ########################################################################################
307	! Use these as pixel counts at first
308	! ########################################################################################
309	Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2))
310	Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
311	Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2))
312	Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure <= &
313	highCloudPressureLimit, dim = 2))
314	Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure > &
315	lowCloudPressureLimit, dim = 2))
316	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
317	- Cloud_Fraction_Low_Mean(1:nPoints)
318
319	! ########################################################################################
320	! Compute column amounts.
321	! ########################################################################################
322	where(Cloud_Fraction_Total_Mean(1:nPoints) > 0)
323	Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / &
324	Cloud_Fraction_Total_Mean(1:nPoints)
325	Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
326	dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
327	elsewhere
328	Optical_Thickness_Total_Mean = R_UNDEF
329	Optical_Thickness_Total_MeanLog10 = R_UNDEF
330	endwhere
331	where(Cloud_Fraction_Water_Mean(1:nPoints) > 0)
332	Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
333	Cloud_Fraction_Water_Mean(1:nPoints)
334	Liquid_Water_Path_Mean(1:nPoints) = LWP_conversionsum(particle_sizeoptical_thickness, &
335	mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
336	Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
337	dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
338	Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
339	Cloud_Fraction_Water_Mean(1:nPoints)
340	elsewhere
341	Optical_Thickness_Water_Mean = R_UNDEF
342	Optical_Thickness_Water_MeanLog10 = R_UNDEF
343	Cloud_Particle_Size_Water_Mean = R_UNDEF
344	Liquid_Water_Path_Mean = R_UNDEF
345	endwhere
346	where(Cloud_Fraction_Ice_Mean(1:nPoints) > 0)
347	Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / &
348	Cloud_Fraction_Ice_Mean(1:nPoints)
349	Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_densitysum(particle_sizeoptical_thickness,&
350	mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
351	Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
352	dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
353	Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / &
354	Cloud_Fraction_Ice_Mean(1:nPoints)
355	elsewhere
356	Optical_Thickness_Ice_Mean = R_UNDEF
357	Optical_Thickness_Ice_MeanLog10 = R_UNDEF
358	Cloud_Particle_Size_Ice_Mean = R_UNDEF
359	Ice_Water_Path_Mean = R_UNDEF
360	endwhere
361	Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
362	max(1, count(cloudMask, dim = 2))
363
364	! ########################################################################################
365	! Normalize pixel counts to fraction.
366	! ########################################################################################
367	Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
368	Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols
369	Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols
370	Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
371	Cloud_Fraction_Ice_Mean(1:nPoints) = Cloud_Fraction_Ice_Mean(1:nPoints) /nSubcols
372	Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
373
374	! ########################################################################################
375	! Joint histograms
376	! ########################################################################################
377	! Loop over all points
378	tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols)
379	ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols)
380	reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
381	reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
382	do j=1,nPoints
383
384	! Fill clear and optically thin subcolumns with fill
385	where(.not. cloudMask(j,1:nSubCols))
386	tauWRK(j,1:nSubCols) = -999._wp
387	ctpWRK(j,1:nSubCols) = -999._wp
388	endwhere
389	! Joint histogram of tau/CTP
390	call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
391	modis_histTau,numMODISTauBins,&
392	modis_histPres,numMODISPresBins,&
393	Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numMODISTauBins,1:numMODISPresBins))
394	! Joint histogram of tau/ReffICE
395	call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, &
396	modis_histTau,numMODISTauBins,modis_histReffIce, &
397	numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numMODISTauBins,1:numMODISReffIceBins))
398	! Joint histogram of tau/ReffLIQ
399	call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, &
400	modis_histTau,numMODISTauBins,modis_histReffLiq, &
401	numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numMODISTauBins,1:numMODISReffLiqBins))
402
403	enddo
404	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numMODISTauBins,1:numMODISPresBins) = &
405	Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numMODISTauBins,1:numMODISPresBins)/nSubCols
406	Optical_Thickness_vs_ReffIce(1:nPoints,1:numMODISTauBins,1:numMODISReffIceBins) = &
407	Optical_Thickness_vs_ReffIce(1:nPoints,1:numMODISTauBins,1:numMODISReffIceBins)/nSubCols
408	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numMODISTauBins,1:numMODISReffLiqBins) = &
409	Optical_Thickness_vs_ReffLiq(1:nPoints,1:numMODISTauBins,1:numMODISReffLiqBins)/nSubCols
410
411
412	! Unit conversion
413	where(Optical_Thickness_vs_Cloud_Top_Pressure /= R_UNDEF) &
414	Optical_Thickness_vs_Cloud_Top_Pressure = Optical_Thickness_vs_Cloud_Top_Pressure*100._wp
415	where(Optical_Thickness_vs_ReffIce /= R_UNDEF) Optical_Thickness_vs_ReffIce = Optical_Thickness_vs_ReffIce*100._wp
416	where(Optical_Thickness_vs_ReffLiq /= R_UNDEF) Optical_Thickness_vs_ReffLiq = Optical_Thickness_vs_ReffLiq*100._wp
417	where(Cloud_Fraction_Total_Mean /= R_UNDEF) Cloud_Fraction_Total_Mean = Cloud_Fraction_Total_Mean*100._wp
418	where(Cloud_Fraction_Water_Mean /= R_UNDEF) Cloud_Fraction_Water_Mean = Cloud_Fraction_Water_Mean*100._wp
419	where(Cloud_Fraction_Ice_Mean /= R_UNDEF) Cloud_Fraction_Ice_Mean = Cloud_Fraction_Ice_Mean*100._wp
420	where(Cloud_Fraction_High_Mean /= R_UNDEF) Cloud_Fraction_High_Mean = Cloud_Fraction_High_Mean*100._wp
421	where(Cloud_Fraction_Mid_Mean /= R_UNDEF) Cloud_Fraction_Mid_Mean = Cloud_Fraction_Mid_Mean*100._wp
422	where(Cloud_Fraction_Low_Mean /= R_UNDEF) Cloud_Fraction_Low_Mean = Cloud_Fraction_Low_Mean*100._wp
423
424	end subroutine modis_column
425
426	! ########################################################################################
427	function cloud_top_pressure(nLevels,tauIncrement, pressure, tauLimit)
428	! INPUTS
429	integer, intent(in) :: nLevels
430	real(wp),intent(in),dimension(nLevels) :: tauIncrement, pressure
431	real(wp),intent(in) :: tauLimit
432	! OUTPUTS
433	real(wp) :: cloud_top_pressure
434	! LOCAL VARIABLES
435	real(wp) :: deltaX, totalTau, totalProduct
436	integer :: i
437
438	! Find the extinction-weighted pressure. Assume that pressure varies linearly between
439	! layers and use the trapezoidal rule.
440	totalTau = 0._wp; totalProduct = 0._wp
441	do i = 2, size(tauIncrement)
442	if(totalTau + tauIncrement(i) > tauLimit) then
443	deltaX = tauLimit - totalTau
444	totalTau = totalTau + deltaX
445	!
446	! Result for trapezoidal rule when you take less than a full step
447	! tauIncrement is a layer-integrated value
448	!
449	totalProduct = totalProduct &
450	+ pressure(i-1) * deltaX &
451	+ (pressure(i) - pressure(i-1)) * deltaX*2/(2._wp tauIncrement(i))
452	else
453	totalTau = totalTau + tauIncrement(i)
454	totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2._wp
455	end if
456	if(totalTau >= tauLimit) exit
457	end do
458
459	if (totalTau > 0._wp) then
460	cloud_top_pressure = totalProduct/totalTau
461	else
462	cloud_top_pressure = 0._wp
463	endif
464
465	end function cloud_top_pressure
466
467	! ########################################################################################
468	function weight_by_extinction(nLevels,tauIncrement, f, tauLimit)
469	! INPUTS
470	integer, intent(in) :: nLevels
471	real(wp),intent(in),dimension(nLevels) :: tauIncrement, f
472	real(wp),intent(in) :: tauLimit
473	! OUTPUTS
474	real(wp) :: weight_by_extinction
475	! LOCAL VARIABLES
476	real(wp) :: deltaX, totalTau, totalProduct
477	integer :: i
478
479	! Find the extinction-weighted value of f(tau), assuming constant f within each layer
480	totalTau = 0._wp; totalProduct = 0._wp
481	do i = 1, size(tauIncrement)
482	if(totalTau + tauIncrement(i) > tauLimit) then
483	deltaX = tauLimit - totalTau
484	totalTau = totalTau + deltaX
485	totalProduct = totalProduct + deltaX * f(i)
486	else
487	totalTau = totalTau + tauIncrement(i)
488	totalProduct = totalProduct + tauIncrement(i) * f(i)
489	end if
490	if(totalTau >= tauLimit) exit
491	end do
492
493	if (totalTau > 0._wp) then
494	weight_by_extinction = totalProduct/totalTau
495	else
496	weight_by_extinction = 0._wp
497	endif
498
499	end function weight_by_extinction
500
501	! ########################################################################################
502	pure function interpolate_to_min(x, y, yobs)
503	! INPUTS
504	real(wp),intent(in),dimension(num_trial_res) :: x, y
505	real(wp),intent(in) :: yobs
506	! OUTPUTS
507	real(wp) :: interpolate_to_min
508	! LOCAL VARIABLES
509	real(wp), dimension(num_trial_res) :: diff
510	integer :: nPoints, minDiffLoc, lowerBound, upperBound
511
512	! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
513	! y must be monotonic in x
514
515	nPoints = size(y)
516	diff(1:num_trial_res) = y(1:num_trial_res) - yobs
517	minDiffLoc = minloc(abs(diff), dim = 1)
518
519	if(minDiffLoc == 1) then
520	lowerBound = minDiffLoc
521	upperBound = minDiffLoc + 1
522	else if(minDiffLoc == nPoints) then
523	lowerBound = minDiffLoc - 1
524	upperBound = minDiffLoc
525	else
526	if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
527	lowerBound = minDiffLoc-1
528	upperBound = minDiffLoc
529	else
530	lowerBound = minDiffLoc
531	upperBound = minDiffLoc + 1
532	end if
533	end if
534
535	if(diff(lowerBound) * diff(upperBound) < 0) then
536	!
537	! Interpolate the root position linearly if we bracket the root
538	!
539	interpolate_to_min = x(upperBound) - &
540	diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
541	else
542	interpolate_to_min = re_fill
543	end if
544
545
546	end function interpolate_to_min
547
548	! ########################################################################################
549	! Optical properties
550	! ########################################################################################
551	elemental function get_g_nir_old (phase, re)
552	! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
553	! of size for ice and water
554	! Fits from Steve Platnick
555
556	! INPUTS
557	integer, intent(in) :: phase
558	real(wp),intent(in) :: re
559	! OUTPUTS
560	real(wp) :: get_g_nir_old
561	! LOCAL VARIABLES(parameters)
562	real(wp), dimension(3), parameter :: &
563	ice_coefficients = (/ 0.7432, 4.5563e-3, -2.8697e-5 /), &
564	small_water_coefficients = (/ 0.8027, -1.0496e-2, 1.7071e-3 /), &
565	big_water_coefficients = (/ 0.7931, 5.3087e-3, -7.4995e-5 /)
566
567	! approx. fits from MODIS Collection 5 LUT scattering calculations
568	if(phase == phaseIsLiquid) then
569	if(re < 8.) then
570	get_g_nir_old = fit_to_quadratic(re, small_water_coefficients)
571	if(re < re_water_min) get_g_nir_old = fit_to_quadratic(re_water_min, small_water_coefficients)
572	else
573	get_g_nir_old = fit_to_quadratic(re, big_water_coefficients)
574	if(re > re_water_max) get_g_nir_old = fit_to_quadratic(re_water_max, big_water_coefficients)
575	end if
576	else
577	get_g_nir_old = fit_to_quadratic(re, ice_coefficients)
578	if(re < re_ice_min) get_g_nir_old = fit_to_quadratic(re_ice_min, ice_coefficients)
579	if(re > re_ice_max) get_g_nir_old = fit_to_quadratic(re_ice_max, ice_coefficients)
580	end if
581
582	end function get_g_nir_old
583
584	! ########################################################################################
585	elemental function get_ssa_nir_old (phase, re)
586	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
587	! of size for ice and water
588	! Fits from Steve Platnick
589
590	! INPUTS
591	integer, intent(in) :: phase
592	real(wp),intent(in) :: re
593	! OUTPUTS
594	real(wp) :: get_ssa_nir_old
595	! LOCAL VARIABLES (parameters)
596	real(wp), dimension(4), parameter :: ice_coefficients = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /)
597	real(wp), dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /)
598
599	! approx. fits from MODIS Collection 5 LUT scattering calculations
600	if(phase == phaseIsLiquid) then
601	get_ssa_nir_old = fit_to_quadratic(re, water_coefficients)
602	if(re < re_water_min) get_ssa_nir_old = fit_to_quadratic(re_water_min, water_coefficients)
603	if(re > re_water_max) get_ssa_nir_old = fit_to_quadratic(re_water_max, water_coefficients)
604	else
605	get_ssa_nir_old = fit_to_cubic(re, ice_coefficients)
606	if(re < re_ice_min) get_ssa_nir_old = fit_to_cubic(re_ice_min, ice_coefficients)
607	if(re > re_ice_max) get_ssa_nir_old = fit_to_cubic(re_ice_max, ice_coefficients)
608	end if
609
610	end function get_ssa_nir_old
611
612	elemental function get_g_nir (phase, re)
613	!
614	! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
615	! of size for ice and water
616	! Fits from Steve Platnick
617	!
618
619	integer, intent(in) :: phase
620	real(wp), intent(in) :: re
621	real(wp) :: get_g_nir
622
623	real(wp), dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
624	small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
625	real(wp), dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
626
627	! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
628	if(phase == phaseIsLiquid) then
629	if(re < 7.) then
630	get_g_nir = fit_to_quadratic(re, small_water_coefficients)
631	if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
632	else
633	get_g_nir = fit_to_cubic(re, big_water_coefficients)
634	if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
635	end if
636	else
637	get_g_nir = fit_to_quadratic(re, ice_coefficients)
638	if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
639	if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
640	end if
641
642	end function get_g_nir
643
644	! --------------------------------------------
645	elemental function get_ssa_nir (phase, re)
646	integer, intent(in) :: phase
647	real(wp), intent(in) :: re
648	real(wp) :: get_ssa_nir
649	!
650	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
651	! of size for ice and water
652	! Fits from Steve Platnick
653	!
654	real(wp), dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
655	real(wp), dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
656
657	! approx. fits from MODIS Collection 6 LUT scattering calculations
658	if(phase == phaseIsLiquid) then
659	get_ssa_nir = fit_to_quadratic(re, water_coefficients)
660	if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
661	if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
662	else
663	get_ssa_nir = fit_to_cubic(re, ice_coefficients)
664	if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
665	if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
666	end if
667
668	end function get_ssa_nir
669
670
671
672	! ########################################################################################
673	pure function fit_to_cubic(x, coefficients)
674	! INPUTS
675	real(wp), intent(in) :: x
676	real(wp), dimension(4), intent(in) :: coefficients
677	! OUTPUTS
678	real(wp) :: fit_to_cubic
679
680	fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
681	end function fit_to_cubic
682
683	! ########################################################################################
684	pure function fit_to_quadratic(x, coefficients)
685	! INPUTS
686	real(wp), intent(in) :: x
687	real(wp), dimension(3), intent(in) :: coefficients
688	! OUTPUTS
689	real(wp) :: fit_to_quadratic
690
691	fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
692	end function fit_to_quadratic
693
694	! ########################################################################################
695	! Radiative transfer
696	! ########################################################################################
697	pure function compute_toa_reflectace(nLevels,tau, g, w0)
698	! This wrapper reports reflectance only and strips out non-cloudy elements from the
699	! calculation
700
701	! INPUTS
702	integer,intent(in) :: nLevels
703	real(wp),intent(in),dimension(nLevels) :: tau, g, w0
704	! OUTPUTS
705	real(wp) :: compute_toa_reflectace
706	! LOCAL VARIABLES
707	logical, dimension(nLevels) :: cloudMask
708	integer, dimension(count(tau(1:nLevels) > 0)) :: cloudIndicies
709	real(wp),dimension(count(tau(1:nLevels) > 0)) :: Refl,Trans
710	real(wp) :: Refl_tot, Trans_tot
711	integer :: i
712
713	cloudMask(1:nLevels) = tau(1:nLevels) > 0.
714	cloudIndicies = pack((/ (i, i = 1, nLevels) /), mask = cloudMask)
715	do i = 1, size(cloudIndicies)
716	call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
717	end do
718
719	call adding_doubling(count(tau(1:nLevels) > 0),Refl(:), Trans(:), Refl_tot, Trans_tot)
720
721	compute_toa_reflectace = Refl_tot
722
723	end function compute_toa_reflectace
724
725	! ########################################################################################
726	pure subroutine two_stream(tauint, gint, w0int, ref, tra)
727	! Compute reflectance in a single layer using the two stream approximation
728	! The code itself is from Lazaros Oreopoulos via Steve Platnick
729	! INPUTS
730	real(wp), intent(in) :: tauint, gint, w0int
731	! OUTPUTS
732	real(wp), intent(out) :: ref, tra
733	! LOCAL VARIABLES
734	! for delta Eddington code
735	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
736	integer, parameter :: beam = 2
737	real(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
738	real(wp) :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
739	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
740
741	! Compute reflectance and transmittance in a single layer using the two stream approximation
742	! The code itself is from Lazaros Oreopoulos via Steve Platnick
743	f = gint**2
744	tau = (1._wp - w0int * f) * tauint
745	w0 = (1._wp - f) * w0int / (1._wp - w0int * f)
746	g = (gint - f) / (1._wp - f)
747
748	! delta-Eddington (Joseph et al. 1976)
749	gamma1 = (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp
750	gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp
751	gamma3 = (2._wp - 3._wpgxmu) / 4._wp
752	gamma4 = 1._wp - gamma3
753
754	if (w0int > minConservativeW0) then
755	! Conservative scattering
756	if (beam == 1) then
757	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
758
759	ref = rh / (1._wp + gamma1 * tau)
760	tra = 1._wp - ref
761	else if(beam == 2) then
762	ref = gamma1tau/(1._wp + gamma1tau)
763	tra = 1._wp - ref
764	endif
765	else
766	! Non-conservative scattering
767	a1 = gamma1 * gamma4 + gamma2 * gamma3
768	a2 = gamma1 * gamma3 + gamma2 * gamma4
769
770	rk = sqrt(gamma12 - gamma22)
771
772	r1 = (1._wp - rk * xmu) * (a2 + rk * gamma3)
773	r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3)
774	r3 = 2._wp * rk (gamma3 - a2 xmu)
775	r4 = (1._wp - (rk * xmu)*2) (rk + gamma1)
776	r5 = (1._wp - (rk * xmu)*2) (rk - gamma1)
777
778	t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4)
779	t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4)
780	t3 = 2._wp * rk * (gamma4 + a1 * xmu)
781	t4 = r4
782	t5 = r5
783
784	beta = -r5 / r4
785
786	e1 = min(rk * tau, 500._wp)
787	e2 = min(tau / xmu, 500._wp)
788
789	if (beam == 1) then
790	den = r4 * exp(e1) + r5 * exp(-e1)
791	ref = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
792	den = t4 * exp(e1) + t5 * exp(-e1)
793	th = exp(-e2)
794	tra = th-thw0(t1exp(e1)-t2exp(-e1)-t3*exp(e2))/den
795	elseif (beam == 2) then
796	ef1 = exp(-e1)
797	ef2 = exp(-2*e1)
798	ref = (gamma2(1._wp-ef2))/((rk+gamma1)(1._wp-beta*ef2))
799	tra = (2._wprkef1)/((rk+gamma1)(1._wp-betaef2))
800	endif
801	end if
802	end subroutine two_stream
803
804	! ########################################################################################
805	elemental function two_stream_reflectance(tauint, gint, w0int)
806	! Compute reflectance in a single layer using the two stream approximation
807	! The code itself is from Lazaros Oreopoulos via Steve Platnick
808
809	! INPUTS
810	real(wp), intent(in) :: tauint, gint, w0int
811	! OUTPUTS
812	real(wp) :: two_stream_reflectance
813	! LOCAL VARIABLES
814	! for delta Eddington code
815	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
816	integer, parameter :: beam = 2
817	real(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
818	real(wp) :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
819	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
820
821	f = gint**2
822	tau = (1._wp - w0int * f) * tauint
823	w0 = (1._wp - f) * w0int / (1._wp - w0int * f)
824	g = (gint - f) / (1._wp - f)
825
826	! delta-Eddington (Joseph et al. 1976)
827	gamma1 = (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp
828	gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp
829	gamma3 = (2._wp - 3._wpgxmu) / 4._wp
830	gamma4 = 1._wp - gamma3
831
832	if (w0int > minConservativeW0) then
833	! Conservative scattering
834	if (beam == 1) then
835	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
836	two_stream_reflectance = rh / (1._wp + gamma1 * tau)
837	elseif (beam == 2) then
838	two_stream_reflectance = gamma1tau/(1._wp + gamma1tau)
839	endif
840
841	else !
842
843	! Non-conservative scattering
844	a1 = gamma1 * gamma4 + gamma2 * gamma3
845	a2 = gamma1 * gamma3 + gamma2 * gamma4
846
847	rk = sqrt(gamma12 - gamma22)
848
849	r1 = (1._wp - rk * xmu) * (a2 + rk * gamma3)
850	r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3)
851	r3 = 2._wp * rk (gamma3 - a2 xmu)
852	r4 = (1._wp - (rk * xmu)*2) (rk + gamma1)
853	r5 = (1._wp - (rk * xmu)*2) (rk - gamma1)
854
855	t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4)
856	t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4)
857	t3 = 2._wp * rk * (gamma4 + a1 * xmu)
858	t4 = r4
859	t5 = r5
860
861	beta = -r5 / r4
862
863	e1 = min(rk * tau, 500._wp)
864	e2 = min(tau / xmu, 500._wp)
865
866	if (beam == 1) then
867	den = r4 * exp(e1) + r5 * exp(-e1)
868	two_stream_reflectance = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
869	elseif (beam == 2) then
870	ef1 = exp(-e1)
871	ef2 = exp(-2*e1)
872	two_stream_reflectance = (gamma2(1._wp-ef2))/((rk+gamma1)(1._wp-beta*ef2))
873	endif
874
875	end if
876	end function two_stream_reflectance
877
878	! ########################################################################################
879	pure subroutine adding_doubling (npts,Refl, Tran, Refl_tot, Tran_tot)
880	! Use adding/doubling formulas to compute total reflectance and transmittance from
881	! layer values
882
883	! INPUTS
884	integer,intent(in) :: npts
885	real(wp),intent(in),dimension(npts) :: Refl,Tran
886	! OUTPUTS
887	real(wp),intent(out) :: Refl_tot, Tran_tot
888	! LOCAL VARIABLES
889	integer :: i
890	real(wp), dimension(npts) :: Refl_cumulative, Tran_cumulative
891
892	Refl_cumulative(1) = Refl(1)
893	Tran_cumulative(1) = Tran(1)
894
895	do i=2, npts
896	! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
897	Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)(Tran_cumulative(i-1)2)/(1._wp - Refl_cumulative(i-1) Refl(i))
898	Tran_cumulative(i) = (Tran_cumulative(i-1)Tran(i)) / (1._wp - Refl_cumulative(i-1) Refl(i))
899	end do
900
901	Refl_tot = Refl_cumulative(size(Refl))
902	Tran_tot = Tran_cumulative(size(Refl))
903
904	end subroutine adding_doubling
905
906	end module mod_modis_sim

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

Download in other formats:

Original Format