; ============================================================================== ; ; latbin.pro ; ; ; ; Modifications/updates: ; ; May 2016 J. D. Lumpe Updated to allow processing of the new ; Version 5 continuous imaging data. ; Jul 2016 J. D. Lumpe Updated for v5.1 L2 algorithm to allow ; processing of both old, pre-CI mode data and ; and new CI mode (NH16 and beyond). ; Apr 2018 J. D. Lumpe New logic added to consistently handle both ; v4 vs v5 cases and CI vs pre-CI. ; AIR data arrays now passed in and output to ; files. ; Jun 2018 J. D. Lumpe Updated to apply the sensitivity arrays to ; data screening. ; Jun - M. Haken Changed to function, processes all latitude ; Jul 2019 bins and all 3 output types (ALL, CLD, ; NOCLD), returning an array of 3 structures. ; No distinction between ascending/descending ; latitudes. ; Input latitudes lat_in and bins lat_grid ; must be as defined in Level 2 CAT data ; (i.e., including co-latitudes for ASC node). ; Bins are now centered on grid latitudes. ; Changed binning logic from where(mask) ; for each latitude bin to calculating ; bin = grid latitudes from lat_in and sorting ; all input by bin sort; binning loop is ; sequentially over unique sort indices ; rather than grid bins, indexing by index ; range into sorted arrays ; 9 Oct 2019 P. E. Meade Updated to properly allow for lat_grid ; latitude spacing other than 1 degree. ; Also, extensively documented based on ; on reverse-engineering analysis of previous ; coding. Changed 'median()' calls to 'mean()'. ; 25 Oct 2019 P. E. Meade Eliminated frac_min. ; 1 Nov 2019 P. E. Meade Eliminated nlayers_min. Removed explicit ; support for version 4 data (always assume ; version 5 or later). ; ============================================================================== function latbin,ci_flag,lat_grid,sza_min,sza_max,nthresh,threshold,nmin,$ rad_min,rad_unc_max,alb_unc_max,iwc_unc_max,def,$ sza_in,lat_in,lon_in,ltime_in,ut_in,alb_in,$ alb_unc_in,rad_in,rad_unc_in,iwc_in,alb_air_in,alb_air_unc_in,$ iwc_air_in,iwc_air_unc_in,nlayers,cld_presence,sen_30,sen_45,$ sen_60,sen_75 ; ; Input arguments: ; ; ci_flag: integer Continuous Imaging (CI) mode ; flag (1 = CI, 0 = not) ; lat_grid: float array Array of latitudes specifying ; bins for which binned values ; are wanted. These are the MID-POINT ; latitudes for each bin. As used, ; typically from some positive lat ; (e.g., 30) up to 89, then from 91 ; to 180 minus the positive lat, e.g., ; [30,31,...,88,89,91,92,...149,150]. ; Note that there are implict ; assumptions that grid will be ; regularly-spaced such that a ; grid latitude SHOULD fall at ; the pole (90) BUT that the grid ; box at the pole will be omitted ; (i.e., the 90 value will be ; lacking from lat_grid. So for a ; 1 degree grid, the peri-polar ; values would be ; [...,89,91,...], for a 2 degree ; grid [...,88,92,...], etc. This ; reflects the fact that no data ; is acquired in close proximity ; to the pole. ; sza_min: float Specified min solar zenith angle ; sza_max: float Specified max solar zenith angle ; nthresh: integer Number of cloud detection threshold ; levels (in terms of albedo) ; threshold: integer array Cloud detection threshold levels ; nmin: integer Generic minimum number of 'detections' ; to allow calculations to proceed ; rad_min: float Minium radius for cloud pixels to ; be included in radius calculations ; (i.e., very small particles are not ; included in mean radius calculations) ; rad_unc_max: float Upper allowed limit on calculated ; cloud particle radius uncertainty ; alb_unc_max: float Upper allowed limit on calculated ; albedo uncertainty ; iwc_unc_max: float Upper allowed limit on calculated ; ice water column uncertainty ; def: float Default value for calculated arrays. ; Usually -999., indicating that the ; value could not be calculated. ; sza_in: float array CIPS L2 Solar Zenith Angle [xdim,ydim]. ; xdim: along-track, ydim: cross-track. ; lat_in: float array CIPS L2 Latitude [xdim,ydim] ; lon_in: float array CIPS L2 Longitude [xdim,ydim] ; ltime_in: array CIPS L2 Local Time [xdim,ydim] ; ut_in: float array CIPS L2 Universal Time [xdim,ydim] ; alb_in: float array CIPS L2 Cloud Albedo [xdim,ydim] ; alb_unc_in: float array CIPS L2 Albedo Uncertainty [xdim,ydim] ; rad_in: float array CIPS L2 Particle Radius [xdim,ydim] ; rad_unc_in: float array CIPS L2 Radius Uncertainty [xdim,ydim] ; iwc_in: float array CIPS L2 Ice Water Content [xdim,ydim] ; alb_air_in: float array CIPS L2 AIR Albedo [xdim,ydim] ; alb_air_unc_in: float array CIPS L2 AIR Albedo Uncert. [xdim,ydim] ; iwc_air_in: float array CIPS L2 AIR Radius [xdim,ydim] ; iwc_air_unc_in: float array CIPS L2 AIR Radius Uncert. [xdim,ydim] ; nlayers: float array CIPS L2 Number of Observations at each ; each location (each obs at a different ; scattering angle) [xdim,ydim] ; cld_presence: integer array CIPS L2 cloud detection map. Binary ; array, 1: cloud, 2: none [xdim,ydim] ; sen_30: float array ; sen_45: float array ; sen_60: float array ; sen_75: float array ; ; ------------ ; Initialize ; ------------ ; nbins is the number of latitude bins in the output grids nbins=n_elements(lat_grid) ; latinc is the size/spacing of latitude bins in the output grids latinc=abs(lat_grid[1]-lat_grid[0]) l3c_struc={num_cld:intarr(nthresh,nbins), $ ; Output data structure num_obs:intarr(nthresh,nbins), $ lon:fltarr(nthresh,nbins)+def, $ sza:fltarr(nthresh,nbins)+def, $ ltime:fltarr(nthresh,nbins)+def, $ ut:fltarr(nthresh,nbins)+def, $ alb:fltarr(nthresh,nbins)+def, $ rad:fltarr(nthresh,nbins)+def, $ iwc:fltarr(nthresh,nbins)+def, $ alb_std:fltarr(nthresh,nbins)+def, $ rad_std:fltarr(nthresh,nbins)+def, $ iwc_std:fltarr(nthresh,nbins)+def, $ alb_air:fltarr(nthresh,nbins)+def, $ iwc_air:fltarr(nthresh,nbins)+def, $ alb_air_std:fltarr(nthresh,nbins)+def, $ iwc_air_std:fltarr(nthresh,nbins)+def} l3c=replicate(l3c_struc,2) ; [*,0]: CLD, [*,1]: ALL ; lat_in is a 2-d array of the geolocation latitude values associated ; with the input data. Adding latinc/2 to each value and applying floor() ; has the effect mapping data from latitudes halfway less than a ; lat_grid value to halfway more than a lat_grid value to that lat_grid ; value (e.g., for 1-degree bins in lat_grid, lat_in values from 54.500... ; to 55.499... map to the lat_grid 55 bin. 'bins' is therefore a 2D array ; giving, for each lat_in value, the lat_grid value to which it maps: bins=floor(lat_in + latinc/2.) ; sort(bins) yields a vector containing the 1D indices into the ; 'bins' array such that the bins[sort(bins)] yields the elements ; of 'bins' in ascending latitude order, as a 1D array. Similarly ; sort(sbins) can be used to sort any of the input data arrays ; in ascending latitude order. 'sbins', as defined below, contains ALL ; of the values of 'bins' in ascending order. If there are multiple ; occurrences of a latitude value in 'bins' there will be the same ; number of occurrences (sequentially) in 'sbins': sbins=bins[sort(bins)] ; uniq() yields the index of the LAST occurrence of each latitude value ; in 'sbins'. iubins is therefore a 1D array of the indicies into the ; 'sbins' vector of the last occurence of each latitude value in 'sbins': iubins=uniq(sbins) ; 'ubins' picks out each of the unique latitude values in 'sbin', because ; it takes (only) the values at each of the indices in iubins. Latitude ; values in ubins will be in ascending order: ubins=sbins[iubins] ; Sort all input arrays by sort(bins) (i.e., ascending lat order). Prefix ; 's' indicates 'sorted': slat_in = lat_in[sort(bins)] ssza_in = sza_in[sort(bins)] scld_presence = cld_presence[sort(bins)] ssen_30 = sen_30[sort(bins)] ssen_45 = sen_45[sort(bins)] ssen_60 = sen_60[sort(bins)] ssen_75 = sen_75[sort(bins)] salb_in = alb_in[sort(bins)] sltime_in = ltime_in[sort(bins)] sut_in = ut_in[sort(bins)] slon_in = lon_in[sort(bins)] salb_unc_in = alb_unc_in[sort(bins)] srad_in = rad_in[sort(bins)] srad_unc_in = rad_unc_in[sort(bins)] siwc_in = iwc_in[sort(bins)] salb_air_in = alb_air_in[sort(bins)] salb_air_unc_in = alb_air_unc_in[sort(bins)] siwc_air_in = iwc_air_in[sort(bins)] siwc_air_unc_in = iwc_air_unc_in[sort(bins)] snlayers = nlayers[sort(bins)] ; crit_sza is an array of binary values (0 or 1) indicating whether or ; not values in the sza array are between sza_min and sza_max: crit_sza = ssza_in ge sza_min and ssza_in le sza_max ; crit_nocld is an array of binary values (0 or 1) indicating whether ; or not values in the scld_presence array are equal to zero (i.e., ; no clouds): crit_nocld = scld_presence eq 0 ; crit_tot is an array of binary values (0 or 1) indicating whether ; or not values in the scld_presence array are greater than or equal ; to zero (i.e., not flagged as 'bad' values): crit_tot = scld_presence ge 0 ; Prepending an '-1' entry to iubins allows iubins to be used to ; pick out the index ranges (low to high) in the input data arrays ; that correspond to a given latitude bin in lat_grid (see below): iubins=[-1,iubins] ; -1 prepended because iubins is the LAST index for each unique value nubins=n_elements(ubins) ; iubinlo and iubinhi are the lower and upper indices of the iubins ; array which correspond to latitudes in the ubins array which are in ; the range of latitudes allowed in the lat_grid output array. The ; '+1' in the calculation accounts for the fact that iubins was ; shifted by one position relative to ubins when the '-1' element was ; prepended to iubins: iubinlo=min(where(ubins ge min(lat_grid)))+1 ; +1 because of the extra -1 prepended to iubins iubinhi=max(where(ubins le max(lat_grid)))+1 for i = iubinlo,iubinhi do begin ; loop over unique bin indices in sorted arrays ; i is an index into the iubins array ; iubins is an array of indices into the sbins array giving the last ; occurrence of the unique elements of sbins ; sbins is a sorted array of the 'integer-ized' latitudes in the ; input data ; SO: ; index i marches through the unique 'integer-ized' latitudes in the input data ; sbins[iubins[i]] yields the integer latitude value corresponding to i ; ilat is the index into the lat_grid array for the containing the ; latitude given by sbins[iubins[i]]. Note: the ; '(sbins[iubins[i]] gt (90 - fix(latinc)))' term accounts for the ; assumption that there will be no lat_grid bin centered over ; the pole: ilat = fix((sbins[iubins[i]] - fix(lat_grid[0]))/latinc) - $ (sbins[iubins[i]] gt (90 - fix(latinc))) ; The '_bin' arrays contain the subrange of each of the ; corresponding sorted input data arrays that fall into the ; output array latitude bin that is being processed: sza_bin=ssza_in[iubins[i-1]+1:iubins[i]] ltime_bin=sltime_in[iubins[i-1]+1:iubins[i]] ut_bin=sut_in[iubins[i-1]+1:iubins[i]] lon_bin=slon_in[iubins[i-1]+1:iubins[i]] alb_bin=salb_in[iubins[i-1]+1:iubins[i]] alb_unc_bin=salb_unc_in[iubins[i-1]+1:iubins[i]] rad_bin=srad_in[iubins[i-1]+1:iubins[i]] rad_unc_bin=srad_unc_in[iubins[i-1]+1:iubins[i]] iwc_bin=siwc_in[iubins[i-1]+1:iubins[i]] alb_air_bin=salb_air_in[iubins[i-1]+1:iubins[i]] alb_air_unc_bin=salb_air_unc_in[iubins[i-1]+1:iubins[i]] iwc_air_bin=siwc_air_in[iubins[i-1]+1:iubins[i]] iwc_air_unc_bin=siwc_air_unc_in[iubins[i-1]+1:iubins[i]] nlayers_bin=snlayers[iubins[i-1]+1:iubins[i]] ; ----------------------------------- ; Loop over albedo threshold values ; ----------------------------------- for ith= 0,nthresh-1 do begin ngd_alb=0 ngd_iwc=0 ngd_rad=0 ngd_alb_air=0 ngd_iwc_air=0 alb_thresh=threshold[ith] ; ----------------------------------- ; Define data screen criteria ; ----------------------------------- ; The 'crit_' (for criterion) arrays are binary-valued arrays ; (each element is either 0 or 1) indicating whether the ; corresponding elements of the subrange input arrays meet the ; criterion. ; crit_sen indicates whether each of sensitivity input subrange ; arrays meet the albedo threshold criterion: crit_sen = ssen_30[iubins[i-1]+1:iubins[i]] le alb_thresh and $ ssen_45[iubins[i-1]+1:iubins[i]] le alb_thresh and $ ssen_60[iubins[i-1]+1:iubins[i]] le alb_thresh and $ ssen_75[iubins[i-1]+1:iubins[i]] le alb_thresh ; crit_min convolves the crit_sen array with the relevant ; subrange of the crit_sza array: crit_min = crit_sza[iubins[i-1]+1:iubins[i]] and crit_sen ; The 'nocld', 'cld', 'low', and 'total' arrays provide indices ; into the subrange input arrays where the corresponding criteria ; are met (in addition to meeting the sensitivity and solar ; zenith angle criteria): ; 'nocld' indicates cloud-free elements ; 'cld' indicates cloud elements with albedo exceeding the ; albedo threshold ; 'low' indicates cloud elements with albedo less than ; the albedo threshold ; 'total' indicates all elements, cloud or cloud-free ; For each array, the number of elements in the subrange input ; arrays meeting the criteria is also calculated (nnocld, ncld, ; nlow, and ntot). nocld=where(crit_min and crit_nocld[iubins[i-1]+1:iubins[i]],nnocld) crit_cld=scld_presence[iubins[i-1]+1:iubins[i]] eq 1 and salb_in[iubins[i-1]+1:iubins[i]] ge alb_thresh cld=where(crit_min and crit_cld,ncld) crit_low=scld_presence[iubins[i-1]+1:iubins[i]] eq 1 and salb_in[iubins[i-1]+1:iubins[i]] lt alb_thresh low=where(crit_min and crit_low,nlow) total=where(crit_min and crit_tot[iubins[i-1]+1:iubins[i]],ntot) ; Fill in the output structure arrays for the 2 output types: ; cloud = structure l3c[0], all = l3c[1] ; ------------------- ; Clouded elements case ; ------------------- ; Only consider the element of the output data arrays (indexed ; by ilat) to be clouded if the number of clouded elements ; in the contributing subrange input arrays exceeds the ; generic minimum count criterion. if(ncld ge nmin) then begin ; Since we are looking at ONLY clouded elements, both ; num_cld and num_obs in the output arrays are set to ncld ; from the subrange input arrays: l3c[0].num_cld[ith,ilat]=ncld l3c[0].num_obs[ith,ilat]=ncld ; For sza, ltime, and ut in the output arrays, pick the ; mean value of the corresponding quantities in the ; contributing subrange input arrays: l3c[0].sza[ith,ilat]=mean(sza_bin[cld]) ; Alt: median() l3c[0].ltime[ith,ilat]=mean(ltime_bin[cld]) ; Alt: median() l3c[0].ut[ith,ilat]=mean(ut_bin[cld]) ; Alt: median() ; For lon in the output arrays, pick the mean value of the ; corresponding quantity in the contributing subrange input ; arrays, accounting for values in the inputs being on both ; sides of the prime meridian: dm=lon_bin[cld] if(max(dm)-min(dm) gt 300) then begin ;crossed prime meridian p=where(dm gt 180) dm(p)-=360. lonavg=mean(dm) ; Alt: median() if(lonavg lt 0) then lonavg+=360. endif else $ lonavg=mean(dm) ; Alt: median() l3c[0].lon[ith,ilat]=lonavg ; ---------------------------------------------------------------------------------- ; Calculate mean albedo. Screen using the error bars. ; ---------------------------------------------------------------------------------- ; alb_bin[cld] is an array containing the clouded subset of the ; subrange input albedo array ; crit_alb is an of binary values (1 or 0), of the same ; dimension as alb_bin[cld], indicating whether the ; corresponding element of alb_bin[cld] meets ; to-be-defined albedo critera ; Initially, crit_alb only indicates whether or not the ; albedo values are 'good' (finite: defined, non-NAN): crit_alb = finite(alb_bin[cld]) ; Convolve crit_alb with the maximum albedo uncertainty criterion: crit_alb=crit_alb and 100.*alb_unc_bin[cld]/alb_bin[cld] le alb_unc_max ; gd_alb contains the indices of the non-zero elements of ; the crit_alb array, so alb_bin[crit_alb[gd_alb]] represents ; the clouded elements of the subrange input albedo array ; meeting the albedo criteria contributing to calculation ; crit_alb. ngd_alb is the count of 'good' elements: gd_alb=where(crit_alb,ngd_alb) ; If the number of 'good' albedo elements in the subrange ; input albedo array meets the generic minimum count ; criterion, for the albedo value in the output array pick ; the mean value of the corresponding 'good' values in the ; contributing subrange input array. Set the output albedo ; standard deviation to the standard deviation of the 'good' ; input values: if(ngd_alb ge nmin) then begin l3c[0].alb[ith,ilat]=mean(alb_bin[cld[gd_alb]]) ; Alt: median() l3c[0].alb_std[ith,ilat]=stddev(alb_bin[cld[gd_alb]]) endif ; ---------------------------------------------------------------------------------- ; Calculate mean IWC and RAD values. ; Screen out cloud pixels with unrealistically low radius values. ; Screen for radius uncertainty below threshold and NLAYERS GE 2. ; ---------------------------------------------------------------------------------- ; rad_bin[cld] is an array containing the clouded subset of the ; subrange input particle radius array ; crit_rad is an of binary values (1 or 0), of the same ; dimension as rad_bin[cld], indicating whether the ; corresponding element of rad_bin[cld] meets ; to-be-defined radius critera ; Initially, crit_rad only indicates whether or not the ; radius values are 'good' (finite: defined, non-NAN) and ; larger than the specified minimum: crit_rad=finite(rad_bin[cld]) and rad_bin[cld] ge rad_min ; Convolve the radius uncertainty criterion with crit_rad: crit_rad=crit_rad and rad_unc_bin[cld] le rad_unc_max ; In continuous imaging mode, convolve crit_rad with a ; minimum nlayers (scattering angle observations) criterion: if(ci_flag) then $ crit_rad=crit_rad and nlayers_bin[cld] ge 2 ; gd_rad contains the indices of the non-zero elements of ; the crit_rad array, so rad_bin[cld[gd_rad]] represents ; the clouded elements of the subrange input radius array ; meeting the radius criteria contributing to calculation of ; crit_rad. ngd_rad is the count of 'good' elements: gd_rad=where(crit_rad,ngd_rad) ; If the number of 'good' radius elements in the subrange ; input radius array meets the generic minimum count ; criterion, for the radius value in the output array pick ; the mean value of the corresponding 'good' values in the ; contributing subrange input array. Set the output radius ; standard deviation to the standard deviation of the 'good' ; input values. Similarly set the ice water column and ice ; water column standard deviation: if(ngd_rad ge nmin) then begin l3c[0].rad[ith,ilat]=mean(rad_bin[cld[gd_rad]]) ; Alt: median() l3c[0].rad_std[ith,ilat]=stddev(rad_bin[cld[gd_rad]]) l3c[0].iwc[ith,ilat]=mean(iwc_bin[cld[gd_rad]]) ; Alt: median() l3c[0].iwc_std[ith,ilat]=stddev(iwc_bin[cld[gd_rad]]) endif ; ---------------------------------------------------------------------------------- ; Calculate mean AIR data products. Screen using the error bars. ; ---------------------------------------------------------------------------------- ; alb_air_bin[cld] is an array containing the clouded subset of the ; subrange input AIR albedo array ; crit_alb_air is an of binary values (1 or 0), of the same ; dimension as alb_air_bin[cld], indicating whether the ; corresponding element of alb_air_bin[cld] meets ; to-be-defined albedo critera ; Initially, crit_alb_air only indicates whether or not the ; AIR albedo values are 'good' (finite: defined, non-NAN): crit_alb_air = finite(alb_air_bin[cld]) ; Convolve crit_alb_air with the maximum albedo uncertainty criterion: crit_alb_air=crit_alb_air and $ 100.*alb_air_unc_bin[cld]/alb_air_bin[cld] le alb_unc_max ; gd_alb_air contains the indices of the non-zero elements of ; the crit_alb_air array, so alb_air_bin[crit_alb[gd_alb]] represents ; the clouded elements of the subrange input AIR albedo array ; meeting the AIR albedo criteria contributing to calculation ; crit_alb_air. ngd_alb_air is the count of 'good' elements: gd_alb_air=where(crit_alb_air,ngd_alb_air) ; If the number of 'good' AIR albedo elements in the subrange ; input AIR albedo array meets the generic minimum count ; criterion, for the AIR albedo value in the output array pick ; the mean value of the corresponding 'good' values in the ; contributing subrange input array. Set the output AIR albedo ; standard deviation to the standard deviation of the 'good' ; input values: if(ngd_alb_air ge nmin) then begin l3c[0].alb_air[ith,ilat]=mean(alb_air_bin[cld[gd_alb_air]]) l3c[0].alb_air_std[ith,ilat]=stddev(alb_air_bin[cld[gd_alb_air]]) endif ; iwc_air_bin[cld] is an array containing the clouded subset of the ; subrange input AIR ice water column array ; crit_iwc_air is an of binary values (1 or 0), of the same ; dimension as iwc_air_bin[cld], indicating whether the ; corresponding element of iwc_air_bin[cld] meets ; to-be-defined ice water column critera ; Initially, crit_iwc_air only indicates whether or not the ; AIR albedo values are 'good' (finite: defined, non-NAN): crit_iwc_air = finite(iwc_air_bin[cld]) eq 1 ; Convolve crit_iwc_air with the maximum ice water column uncertainty criterion: crit_iwc_air=crit_iwc_air and $ 100.*iwc_air_unc_bin[cld]/iwc_air_bin[cld] le iwc_unc_max ; gd_iwc_air contains the indices of the non-zero elements of ; the crit_iwc_air array, so iwc_air_bin[crit_alb[gd_alb]] represents ; the clouded elements of the subrange input AIR ice water column array ; meeting the AIR ice water column criteria contributing to calculation ; crit_iwc_air. ngd_iwc_air is the count of 'good' elements: gd_iwc_air=where(crit_iwc_air,ngd_iwc_air) ; If the number of 'good' AIR ice water column elements in the subrange ; input AIR ice water column array meets the generic minimum count ; criterion, for the AIR ice water column value in the output array pick ; the mean value of the corresponding 'good' values in the ; contributing subrange input array. Set the output AIR ice ; water column standard deviation to the standard deviation of the 'good' ; input values: if(ngd_iwc_air ge nmin) then begin l3c[0].iwc_air[ith,ilat]=mean(iwc_air_bin[cld[gd_iwc_air]]) l3c[0].iwc_air_std[ith,ilat]=stddev(iwc_air_bin[cld[gd_iwc_air]]) endif endif ; ------------------- ; All elements (cloudy and unclouded) case ; ; Note that for this case, only the mean values of albdeo and ; ice water column are of interest, so neither the standard ; deviations of the quantities nor the radius quantities are calculated. ; ------------------- if(ntot ge nmin) then begin ; Set the total number of elements and the number of clouded ; elements in the subrange input arrays for the output ; data structure: l3c[1].num_cld[ith,ilat]=ncld l3c[1].num_obs[ith,ilat]=ntot ; For sza, ltime, and ut in the output arrays, pick the ; mean value of the corresponding quantities in the ; contributing subrange input arrays: l3c[1].sza[ith,ilat]=mean(sza_bin[total]) ; Alt: median() l3c[1].ltime[ith,ilat]=mean(ltime_bin[total]) ; Alt: median() l3c[1].ut[ith,ilat]=mean(ut_bin[total]) ; Alt: median() ; For lon in the output arrays, pick the mean value of the ; corresponding quantity in the contributing subrange input ; arrays, accounting for values in the inputs being on both ; sides of the prime meridian: dm=lon_bin[total] if(max(dm)-min(dm) gt 300) then begin ;crossed 0 lon line p=where(dm gt 180) dm(p)-=360. lonavg=mean(dm) ; Alt: median() if(lonavg lt 0) then lonavg+=360. endif else $ lonavg=mean(dm) ; Alt: median() l3c[1].lon[ith,ilat]=lonavg if(ncld eq 0) then begin ; If ncld is zero, then the 'all elements' case is ; degenerate with the unclouded case and the albedo, ice ; water column, AIR albedo, and AIR ice water column are ; set to zero in the output data structure. l3c[1].alb[ith,ilat]=0.0 l3c[1].iwc[ith,ilat]=0.0 l3c[1].alb_air[ith,ilat]=0.0 l3c[1].iwc_air[ith,ilat]=0.0 endif else begin ; ---------------------------------------------------------------------------------- ; Calculate mean albedo. Screen using the error bars. ; ---------------------------------------------------------------------------------- ; alb_bin[cld] is an array containing the clouded subset of the ; subrange input albedo array ; crit_alb is an of binary values (1 or 0), of the same ; dimension as alb_bin[cld], indicating whether the ; corresponding element of alb_bin[cld] meets ; to-be-defined albedo critera ; Initially, crit_alb only indicates whether or not the ; albedo values are 'good' (finite: defined, non-NAN): crit_alb = finite(alb_bin[cld]) ; Convolve crit_alb with the maximum albedo uncertainty criterion: crit_alb=crit_alb and 100.*alb_unc_bin[cld]/alb_bin[cld] le alb_unc_max ; gd_alb contains the indices of the non-zero elements of ; the crit_alb array, so alb_bin[crit_alb[gd_alb]] represents ; the clouded elements of the subrange input albedo array ; meeting the albedo criteria contributing to calculation ; crit_alb. ngd_alb is the count of 'good' elements: gd_alb=where(crit_alb,ngd_alb) ; If the number of 'good' albedo elements in the subrange ; input albedo array is non-zero, fill a temporary array ; with the corresponding 'good' values in the ; contributing subrange input array. Otherwise, just set ; the temporary 'array' to scalar zero: if(ngd_alb ne 0) then dum_alb=alb_bin[cld[gd_alb]] else dum_alb=0. ; ----------------------------------------------------------------------------------------- ; Calculate mean IWC. Screen out low radius pixels. ; Screen for radius uncertainty below threshold and NLAYERS GE 2 (CI data only). ; ----------------------------------------------------------------------------------------- ; iwc_bin[cld] is an array containing the clouded subset of the ; subrange input ice water column array ; crit_iwc is an of binary values (1 or 0), of the same ; dimension as iwc_bin[cld], indicating whether the ; corresponding element of iwc_bin[cld] meets ; to-be-defined critera ; Initially, crit_iwc only indicates whether or not the ; ice water column values are 'good' (finite: defined, ; non-NAN) and whether the corresponding radius values ; are larger than the specified minimum: crit_iwc=finite(iwc_bin[cld]) and rad_bin[cld] ge rad_min ; Convolve the radius uncertainty criterion with crit_iwc: crit_iwc=crit_iwc and rad_unc_bin[cld] le rad_unc_max ; In continuous integration mode, convolve crit_iwc with a ; minimum nlayers (scattering angle observations) criterion: if(ci_flag) then $ crit_iwc=crit_iwc and nlayers_bin[cld] ge 2 ; gd_iwc contains the indices of the non-zero elements of ; the crit_iwc array, so iwc_bin[cld[gd_rad]] represents ; the clouded elements of the subrange input radius array ; meeting the radius criteria contributing to calculation of ; crit_iwc. ngd_iwc is the count of 'good' elements: gd_iwc=where(crit_iwc,ngd_iwc) ; If the number of 'good' ice water column elements in ; the subrange input ice water column array is non-zero, ; fill a temporary array with the corresponding 'good' ; values in the contributing subrange input array. ; Otherwise, just set the temporary 'array' to scalar zero: if(ngd_iwc ne 0) then dum_iwc=iwc_bin[cld[gd_iwc]] else dum_iwc=0. ; ---------------------------------------------------------------------------------- ; Calculate mean AIR data products. Screen using the error bars. ; ---------------------------------------------------------------------------------- crit_alb_air = finite(alb_air_bin[cld]) crit_alb_air=crit_alb_air and $ 100.*alb_air_unc_bin[cld]/alb_air_bin[cld] le alb_unc_max gd_alb_air=where(crit_alb_air,ngd_alb_air) if(ngd_alb_air ne 0) then $ dum_alb_air=alb_air_bin[cld[gd_alb_air]] else dum_alb_air=0. ; alb_air_bin[cld] is an array containing the clouded subset of the ; subrange input AIR albedo array ; crit_alb_air is an of binary values (1 or 0), of the same ; dimension as alb_air_bin[cld], indicating whether the ; corresponding element of alb_air_bin[cld] meets ; to-be-defined albedo critera ; Initially, crit_alb_air only indicates whether or not the ; AIR albedo values are 'good' (finite: defined, non-NAN): crit_alb_air = finite(alb_air_bin[cld]) ; Convolve crit_alb_air with the maximum albedo uncertainty criterion: crit_alb_air=crit_alb_air and $ 100.*alb_air_unc_bin[cld]/alb_air_bin[cld] le alb_unc_max ; gd_alb_air contains the indices of the non-zero elements of ; the crit_alb_air array, so alb_air_bin[crit_alb[gd_alb]] represents ; the clouded elements of the subrange input AIR albedo array ; meeting the AIR albedo criteria contributing to calculation ; crit_alb_air. ngd_alb_air is the count of 'good' elements: gd_alb_air=where(crit_alb_air,ngd_alb_air) ; If the number of 'good' AIR albedo elements in the subrange ; input AIR albedo array is non-zero, fill a temporary array ; with the corresponding 'good' values in the contributing ; subrange input array. Otherwise, just set the temporary ; 'array' to scalar zero: if(ngd_alb_air ne 0) then $ dum_alb_air=alb_air_bin[cld[gd_alb_air]] else dum_alb_air=0. ; iwc_air_bin[cld] is an array containing the clouded subset of the ; subrange input AIR ice water column array ; crit_iwc_air is an of binary values (1 or 0), of the same ; dimension as iwc_air_bin[cld], indicating whether the ; corresponding element of iwc_air_bin[cld] meets ; to-be-defined ice water column critera ; Initially, crit_iwc_air only indicates whether or not the ; AIR albedo values are 'good' (finite: defined, non-NAN): crit_iwc_air = finite(iwc_air_bin[cld]) eq 1 ; Convolve crit_iwc_air with the maximum ice water column uncertainty criterion: crit_iwc_air=crit_iwc_air and $ 100.*iwc_air_unc_bin[cld]/iwc_air_bin[cld] le iwc_unc_max ; gd_iwc_air contains the indices of the non-zero elements of ; the crit_iwc_air array, so iwc_air_bin[crit_alb[gd_alb]] represents ; the clouded elements of the subrange input AIR ice water column array ; meeting the AIR ice water column criteria contributing to calculation ; crit_iwc_air. ngd_iwc_air is the count of 'good' elements: gd_iwc_air=where(crit_iwc_air,ngd_iwc_air) ; If the number of 'good' AIR ice water column elements in ; the subrange input AIR ice water column array is non-zero, ; fill a temporary array with corresponding 'good' values ; in the contributing subrange input array. Otherwise, ; just set the temporary 'array' to scalar zero: if(ngd_iwc_air ne 0) then $ dum_iwc_air=iwc_air_bin[cld[gd_iwc_air]] else dum_iwc_air=0. ; ---------------------------------------- ; Now add in the non-cloud/zero pixels ; ---------------------------------------- ; Note: This just adds the needed number of zero elements ; to the ends of the temporary arrays. if(nnocld+nlow gt 0) then begin zeros=fltarr(nnocld+nlow) dum_alb=[dum_alb,zeros] dum_iwc=[dum_iwc,zeros] dum_alb_air=[dum_alb_air,zeros] dum_iwc_air=[dum_iwc_air,zeros] endif ; Set the values of the output data arrays to the mean ; values of the temporary arrays (which come from the ; input data arrays): if(n_elements(dum_alb) ne 0) then l3c[1].alb[ith,ilat]=mean(dum_alb) if(n_elements(dum_iwc) ne 0) then l3c[1].iwc[ith,ilat]=mean(dum_iwc) if(n_elements(dum_alb_air) ne 0) then l3c[1].alb_air[ith,ilat]=mean(dum_alb_air) if(n_elements(dum_iwc_air) ne 0) then l3c[1].iwc_air[ith,ilat]=mean(dum_iwc_air) endelse endif endfor ; loop over threshold values endfor ; end loop over unique bin indices return, l3c end