Context Navigation

modis_simulator.F90 @ 5017

Last change on this file since 5017 was 3358, checked in by idelkadi, 6 years ago
Implementation de la nouvelle version COSPv2 dans LMDZ. Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"
File size: 44.2 KB

Rev	Line
[3358]	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 = 0._wp
	329	Optical_Thickness_Total_MeanLog10 = 0._wp
	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 = 0._wp
	342	Optical_Thickness_Water_MeanLog10 = 0._wp
	343	Cloud_Particle_Size_Water_Mean = 0._wp
	344	Liquid_Water_Path_Mean = 0._wp
	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 = 0._wp
	357	Optical_Thickness_Ice_MeanLog10 = 0._wp
	358	Cloud_Particle_Size_Ice_Mean = 0._wp
	359	Ice_Water_Path_Mean = 0._wp
	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