Context Navigation

modis_simulator.F90 @ 5209

Last change on this file since 5209 was 5185, checked in by abarral, 9 days ago
Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
File size: 44.1 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	integer, intent(in) :: phase
619	real(wp), intent(in) :: re
620	real(wp) :: get_g_nir
621
622	real(wp), dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
623	small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
624	real(wp), dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
625
626	! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
627	if(phase == phaseIsLiquid) then
628	if(re < 7.) then
629	get_g_nir = fit_to_quadratic(re, small_water_coefficients)
630	if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
631	else
632	get_g_nir = fit_to_cubic(re, big_water_coefficients)
633	if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
634	end if
635	else
636	get_g_nir = fit_to_quadratic(re, ice_coefficients)
637	if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
638	if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
639	end if
640
641	end function get_g_nir
642
643	! --------------------------------------------
644	elemental function get_ssa_nir (phase, re)
645	integer, intent(in) :: phase
646	real(wp), intent(in) :: re
647	real(wp) :: get_ssa_nir
648
649	! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
650	! of size for ice and water
651	! Fits from Steve Platnick
652
653	real(wp), dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
654	real(wp), dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
655
656	! approx. fits from MODIS Collection 6 LUT scattering calculations
657	if(phase == phaseIsLiquid) then
658	get_ssa_nir = fit_to_quadratic(re, water_coefficients)
659	if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
660	if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
661	else
662	get_ssa_nir = fit_to_cubic(re, ice_coefficients)
663	if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
664	if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
665	end if
666
667	end function get_ssa_nir
668
669
670
671	! ########################################################################################
672	pure function fit_to_cubic(x, coefficients)
673	! INPUTS
674	real(wp), intent(in) :: x
675	real(wp), dimension(4), intent(in) :: coefficients
676	! OUTPUTS
677	real(wp) :: fit_to_cubic
678
679	fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
680	end function fit_to_cubic
681
682	! ########################################################################################
683	pure function fit_to_quadratic(x, coefficients)
684	! INPUTS
685	real(wp), intent(in) :: x
686	real(wp), dimension(3), intent(in) :: coefficients
687	! OUTPUTS
688	real(wp) :: fit_to_quadratic
689
690	fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
691	end function fit_to_quadratic
692
693	! ########################################################################################
694	! Radiative transfer
695	! ########################################################################################
696	pure function compute_toa_reflectace(nLevels,tau, g, w0)
697	! This wrapper reports reflectance only and strips out non-cloudy elements from the
698	! calculation
699
700	! INPUTS
701	integer,intent(in) :: nLevels
702	real(wp),intent(in),dimension(nLevels) :: tau, g, w0
703	! OUTPUTS
704	real(wp) :: compute_toa_reflectace
705	! LOCAL VARIABLES
706	logical, dimension(nLevels) :: cloudMask
707	integer, dimension(count(tau(1:nLevels) > 0)) :: cloudIndicies
708	real(wp),dimension(count(tau(1:nLevels) > 0)) :: Refl,Trans
709	real(wp) :: Refl_tot, Trans_tot
710	integer :: i
711
712	cloudMask(1:nLevels) = tau(1:nLevels) > 0.
713	cloudIndicies = pack((/ (i, i = 1, nLevels) /), mask = cloudMask)
714	DO i = 1, size(cloudIndicies)
715	call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
716	end do
717
718	call adding_doubling(count(tau(1:nLevels) > 0),Refl(:), Trans(:), Refl_tot, Trans_tot)
719
720	compute_toa_reflectace = Refl_tot
721
722	end function compute_toa_reflectace
723
724	! ########################################################################################
725	pure subroutine two_stream(tauint, gint, w0int, ref, tra)
726	! Compute reflectance in a single layer using the two stream approximation
727	! The code itself is from Lazaros Oreopoulos via Steve Platnick
728	! INPUTS
729	real(wp), intent(in) :: tauint, gint, w0int
730	! OUTPUTS
731	real(wp), intent(out) :: ref, tra
732	! LOCAL VARIABLES
733	! for delta Eddington code
734	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
735	integer, parameter :: beam = 2
736	real(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
737	real(wp) :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
738	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
739
740	! Compute reflectance and transmittance in a single layer using the two stream approximation
741	! The code itself is from Lazaros Oreopoulos via Steve Platnick
742	f = gint**2
743	tau = (1._wp - w0int * f) * tauint
744	w0 = (1._wp - f) * w0int / (1._wp - w0int * f)
745	g = (gint - f) / (1._wp - f)
746
747	! delta-Eddington (Joseph et al. 1976)
748	gamma1 = (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp
749	gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp
750	gamma3 = (2._wp - 3._wpgxmu) / 4._wp
751	gamma4 = 1._wp - gamma3
752
753	if (w0int > minConservativeW0) then
754	! Conservative scattering
755	if (beam == 1) then
756	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
757
758	ref = rh / (1._wp + gamma1 * tau)
759	tra = 1._wp - ref
760	else if(beam == 2) then
761	ref = gamma1tau/(1._wp + gamma1tau)
762	tra = 1._wp - ref
763	endif
764	else
765	! Non-conservative scattering
766	a1 = gamma1 * gamma4 + gamma2 * gamma3
767	a2 = gamma1 * gamma3 + gamma2 * gamma4
768
769	rk = sqrt(gamma12 - gamma22)
770
771	r1 = (1._wp - rk * xmu) * (a2 + rk * gamma3)
772	r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3)
773	r3 = 2._wp * rk (gamma3 - a2 xmu)
774	r4 = (1._wp - (rk * xmu)*2) (rk + gamma1)
775	r5 = (1._wp - (rk * xmu)*2) (rk - gamma1)
776
777	t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4)
778	t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4)
779	t3 = 2._wp * rk * (gamma4 + a1 * xmu)
780	t4 = r4
781	t5 = r5
782
783	beta = -r5 / r4
784
785	e1 = min(rk * tau, 500._wp)
786	e2 = min(tau / xmu, 500._wp)
787
788	if (beam == 1) then
789	den = r4 * exp(e1) + r5 * exp(-e1)
790	ref = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
791	den = t4 * exp(e1) + t5 * exp(-e1)
792	th = exp(-e2)
793	tra = th-thw0(t1exp(e1)-t2exp(-e1)-t3*exp(e2))/den
794	elseif (beam == 2) then
795	ef1 = exp(-e1)
796	ef2 = exp(-2*e1)
797	ref = (gamma2(1._wp-ef2))/((rk+gamma1)(1._wp-beta*ef2))
798	tra = (2._wprkef1)/((rk+gamma1)(1._wp-betaef2))
799	endif
800	end if
801	end subroutine two_stream
802
803	! ########################################################################################
804	elemental function two_stream_reflectance(tauint, gint, w0int)
805	! Compute reflectance in a single layer using the two stream approximation
806	! The code itself is from Lazaros Oreopoulos via Steve Platnick
807
808	! INPUTS
809	real(wp), intent(in) :: tauint, gint, w0int
810	! OUTPUTS
811	real(wp) :: two_stream_reflectance
812	! LOCAL VARIABLES
813	! for delta Eddington code
814	! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
815	integer, parameter :: beam = 2
816	real(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
817	real(wp) :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
818	rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
819
820	f = gint**2
821	tau = (1._wp - w0int * f) * tauint
822	w0 = (1._wp - f) * w0int / (1._wp - w0int * f)
823	g = (gint - f) / (1._wp - f)
824
825	! delta-Eddington (Joseph et al. 1976)
826	gamma1 = (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp
827	gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp
828	gamma3 = (2._wp - 3._wpgxmu) / 4._wp
829	gamma4 = 1._wp - gamma3
830
831	if (w0int > minConservativeW0) then
832	! Conservative scattering
833	if (beam == 1) then
834	rh = (gamma1tau+(gamma3-gamma1xmu)*(1-exp(-tau/xmu)))
835	two_stream_reflectance = rh / (1._wp + gamma1 * tau)
836	elseif (beam == 2) then
837	two_stream_reflectance = gamma1tau/(1._wp + gamma1tau)
838	endif
839
840	else !
841
842	! Non-conservative scattering
843	a1 = gamma1 * gamma4 + gamma2 * gamma3
844	a2 = gamma1 * gamma3 + gamma2 * gamma4
845
846	rk = sqrt(gamma12 - gamma22)
847
848	r1 = (1._wp - rk * xmu) * (a2 + rk * gamma3)
849	r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3)
850	r3 = 2._wp * rk (gamma3 - a2 xmu)
851	r4 = (1._wp - (rk * xmu)*2) (rk + gamma1)
852	r5 = (1._wp - (rk * xmu)*2) (rk - gamma1)
853
854	t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4)
855	t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4)
856	t3 = 2._wp * rk * (gamma4 + a1 * xmu)
857	t4 = r4
858	t5 = r5
859
860	beta = -r5 / r4
861
862	e1 = min(rk * tau, 500._wp)
863	e2 = min(tau / xmu, 500._wp)
864
865	if (beam == 1) then
866	den = r4 * exp(e1) + r5 * exp(-e1)
867	two_stream_reflectance = w0(r1exp(e1)-r2exp(-e1)-r3exp(-e2))/den
868	elseif (beam == 2) then
869	ef1 = exp(-e1)
870	ef2 = exp(-2*e1)
871	two_stream_reflectance = (gamma2(1._wp-ef2))/((rk+gamma1)(1._wp-beta*ef2))
872	endif
873
874	end if
875	end function two_stream_reflectance
876
877	! ########################################################################################
878	pure subroutine adding_doubling (npts,Refl, Tran, Refl_tot, Tran_tot)
879	! Use adding/doubling formulas to compute total reflectance and transmittance from
880	! layer values
881
882	! INPUTS
883	integer,intent(in) :: npts
884	real(wp),intent(in),dimension(npts) :: Refl,Tran
885	! OUTPUTS
886	real(wp),intent(out) :: Refl_tot, Tran_tot
887	! LOCAL VARIABLES
888	integer :: i
889	real(wp), dimension(npts) :: Refl_cumulative, Tran_cumulative
890
891	Refl_cumulative(1) = Refl(1)
892	Tran_cumulative(1) = Tran(1)
893
894	DO i=2, npts
895	! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
896	Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)(Tran_cumulative(i-1)2)/(1._wp - Refl_cumulative(i-1) Refl(i))
897	Tran_cumulative(i) = (Tran_cumulative(i-1)Tran(i)) / (1._wp - Refl_cumulative(i-1) Refl(i))
898	end do
899
900	Refl_tot = Refl_cumulative(size(Refl))
901	Tran_tot = Tran_cumulative(size(Refl))
902
903	end subroutine adding_doubling
904
905	end module mod_modis_sim

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

Download in other formats:

Original Format