;$Id$
@apply_filters
@write_3d_table
@compute_class_cond
@make_true_phase
@ecm_sfc_type
@make_acc_metrics
;--------------------------------------------------------------------------------------------
;  probabilistic phase for 2d naive bayesian
;
;--------------------------------------------------------------------------------------------
pro prob_mask_phase_3d, day = day, night=night, sensor_name = sensor_name, noplot=noplot,  $
                        new=new, old=old, test_name=test_name, sfc_type_on_flag_input = sfc_type_on_flag_input, $
                        path_dir = path_dir, nomod = nomod, Time_Diff_Max = Time_Diff_Max, $
                        cod_clear_ice_thresh = cod_clear_ice_thresh, cod_clear_water_thresh = cod_clear_water_thresh, yesplot=yesplot, $
			isccp=isccp

if (~keyword_set(path_dir)) then path_dir='./'
if (~keyword_set(sensor_name)) then sensor_name='myd02ssh'
if (~keyword_set(Time_Diff_Max)) then Time_Diff_Max = 10.0  ; minutes
if (~keyword_set(cod_clear_ice_thresh)) then cod_clear_ice_thresh = 0.1
if (~keyword_set(cod_clear_water_thresh)) then cod_clear_water_thresh = 0.1
if (sensor_name eq 'myd02ssh') then training_source = 'ecm_training_myd02ssh_train.sav'
if (sensor_name eq 'viirs') then training_source = 'ecm_training_viirs_train_snpp_5+333_2013.sav'
if (sensor_name eq 'abi') then training_source = 'ecm_training_abi_train.sav'
if (sensor_name eq 'abhi') then training_source = 'ecm_training_abhi_train.sav'
if (sensor_name eq 'avhrr') then training_source = 'ecm_training_avhrr_train.sav'
if (sensor_name eq 'mtsat2') then training_source = 'ecm_training_mtsat2_train.sav'
if (sensor_name eq 'vgac') then training_source = 'ecm_training_vgac_train.sav'
if (sensor_name eq 'seviri') then training_source = 'ecm_training_seviri_train.sav'

;--- check for sensor and test inconsistencies
if (strpos(sensor_name,'viirs') lt 0 and strpos(test_name,'dnb') gt 0) then begin
  print, 'DNB classifiers impossible'
  return
endif

;-- control training data modification in make_true_phase routine
mod_flag = 1
if (keyword_set(nomod)) then mod_flag = 0

;-------------------------------
; read in data
;-------------------------------
training_path = './training_data/'
restore, training_path + training_source
f = training_data
training_data = !null
ntrain = f.length

nb_rank = 3

nsfc = 7

if (~keyword_set(sfc_type_on_flag_input)) then sfc_type_on_flag_input = make_array(nsfc,/INT,VALUE=1)

;---- some parameters
laplace_flag = 0
max_r = 1000.0
count_min = 10.0
missing = -999.0
rut_solzen_thresh = 80.0

;--- partial cloud filter
cfrac_min = 0.1
cfrac_max = 0.9

;--- set any geographical limits
lat_min = -90.0
lat_max = 90.0

;---- adjust nasa dnb to look like noaa dnb [REMOVE WHEN ADJUSTMENT in TRAINING DATA]
Dnb_Coef = [-0.373685,0.977945,-0.000261637]
Dnb_Coef = [-0.118767,0.962452,-0.000144502] ; Denis

idx = where(f.refl_lunar_dnb_nom ne missing,cc)

if (cc gt 0) then begin
   f[idx].refl_lunar_dnb_nom = Dnb_Coef[0] + $
                              Dnb_Coef[1]*f[idx].refl_lunar_dnb_nom + $
                              Dnb_Coef[2]*f[idx].refl_lunar_dnb_nom^2
endif

;--- product filtering
cod_cloud_thresh = 1.0
etrop_cloud_thresh = 0.9

;---  default classifier limits - only modify in classifier definition if needed
zen_min = 0.0
zen_max = 90.0
solzen_min = 0.0
solzen_max = 180.0
solzen_day_max = 85.0
solglintzen_min = 0.0
solglintzen_max = 180.0
solscatang_min = 0.0
solscatang_max = 180.0
solscatang_day_min = 75.0
solglint_mask_min = 0
solglint_mask_max = 1
lunzen_min = 0.0
lunzen_max = 180.0
lunglintzen_min = 0.0
lunglintzen_max = 180.0
lunscatang_min = 0.0
lunscatang_max = 180.0
lunglint_mask_min = 0
lunglint_mask_max = 1
coast_mask_min = 0
coast_mask_max = 1
city_mask_min = 0
city_mask_max = 1
moon_illum_frac_min = 0.0
moon_illum_frac_max = 100.0
slant_tpw_min = 0.0
slant_tpw_max = 100.0
tsfc_min = 180.0
tsfc_max = 340.0
zsfc_min = -100.0
zsfc_max = 20000.0
zsfc_std_min = 00.0
zsfc_std_max = 1000.0
snow_class_min = 1
snow_class_max = 3
sfc_type_on_flag = sfc_type_on_flag_input
zsfc_std_max_iruni = 300.0

;----- set solar limits
day_string = 'all'
solzen_min = 0.0
solzen_max = 180.0
if (keyword_set(day)) then begin
   day_string = 'day'
   solzen_min = 0.0
   solzen_max = solzen_day_max
endif
if (keyword_set(night)) then begin
   day_string = 'night'
   solzen_min = 90.0
   solzen_max = 180.0
endif
;------------------------------------------------------------------------------
; small positive (near zero) dnb reflectances should be missing
;------------------------------------------------------------------------------
idx = where(f.refl_lunar_dnb_nom lt 0.1,cc) & if (cc gt 0) then f[idx].refl_lunar_dnb_nom = missing

;------------------------------------------------------------------------------
; loop through classifiers
;------------------------------------------------------------------------------
nc_header = '3D Mask and Phase LUT'

print, '==============================================='
print , "3D Classifier =  ", test_name
print, '==============================================='

case test_name of

"etropo11_topa_logbt11std": begin
   classifier_name ='etropo11_topa_logbt11std'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'topa'
   y = f.cld_temp_opaque
   y_min = 180.0 & y_bin = 5.0 & nbins_y = 32
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 1 & channel_wvl_used = [11000]
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end

"zopa_topa_logbt11std": begin
   classifier_name ='zopa_topa_logbt11std'
   xtitle = 'zopa' & x = f.cld_height_opaque & x_min = 0.0 & x_bin = 1000.0 & nbins_x = 20
   ytitle = 'topa' & y = f.cld_temp_opaque & y_min = 180.0 & y_bin = 5.0 & nbins_y = 32
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 1 & channel_wvl_used = [11000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"etropo11_zopa_logbt11std": begin
   classifier_name ='etropo11_zopa_logbt11std'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'zopa'
   y = f.cld_height_opaque
   y_min = 0.0 & y_bin = 1000.0 & nbins_y = 20
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 1 & channel_wvl_used = [11000]
   coast_mask_max = 0
   ;zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end

"logzopa_refldnb_refl38": begin
   classifier_name = 'logzopa_refldnb_refl38'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'refldnb' & y = f.refl_lunar_dnb_nom & y_min = -5 & y_bin = 2.5 & nbins_y = 50
   ztitle = 'refl38' & z = f.refl_3_75um_nom  & z_min = -20.0 & z_bin = 0.5 & nbins_z = 50.0
   nchan_used = 3 & channel_wvl_used = [650,3750,1100]
   tsfc_min = 240.0
   solzen_min = 90.0
   solzen_max = 180.0
   lunzen_min = 0.0
   lunzen_max = 80.0
   lunglint_mask_min = 0
   lunglint_mask_max = 0
   moon_illum_frac_min = 25.0
   moon_illum_frac_max = 100.0
end

"logzopa_refldnb_btd3811": begin
   classifier_name = 'logzopa_refldnb_btd3811'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'refldnb' & y = f.refl_lunar_dnb_nom & y_min = -5 & y_bin = 2.5 & nbins_y = 50
   ztitle = 'btd3811' & z = f.temp_3_75um_nom - f.temp_11_0um_nom & z_min = -10 & z_bin = 1.0 & nbins_z = 40.0
   nchan_used = 3 & channel_wvl_used = [650,3750,1100]
   tsfc_min = 240.0
   solzen_min = 90.0
   solzen_max = 180.0
   lunzen_min = 0.0
   lunzen_max = 80.0
   lunglint_mask_min = 0
   lunglint_mask_max = 0
   moon_illum_frac_min = 25.0
   moon_illum_frac_max = 100.0
end


"bt73_btd1112_logbt11std": begin
   classifier_name = 'bt73_btd1112_logbt11std'
   xtitle = 'bt73' & x = f.temp_7_3um_nom & x_min = 190.0 & x_bin = 2.0 & nbins_x = 50.0
   ytitle = 'btd1112' & y = f.temp_11_0um_nom - f.temp_12_0um_nom & y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 3 & channel_wvl_used = [7300,11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"logzopa_btd1112_logbt11std": begin
   classifier_name = 'logzopa_btd1112_logbt11std'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'btd1112' & y = f.temp_11_0um_nom - f.temp_12_0um_nom & y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"logzopa_btd11133_logbt11std": begin
   classifier_name = 'logzopa_btd1112_logbt11std'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'btd11133' & y = f.temp_11_0um_nom - f.temp_13_3um_nom & y_min = -10.0 & y_bin = 2.0 & nbins_y = 40
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,13300]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"etropo11_btd1112_logbt11std": begin
   classifier_name = 'etropo11_btd1112_logbt11std'
   xtitle = 'etropo11' & x = f.emiss_tropo_11_0um_nom & x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'btd1112' & y = f.temp_11_0um_nom - f.temp_12_0um_nom & y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt11ratio_btd1112_logbt11std": begin
   classifier_name = 'bt11ratio_btd1112_logbt11std'
   xtitle = 'bt11ratio' 
   x = btratio(f.temp_11_0um_nom,f.surface_temperature_nwp,f.tropopause_temperature_nwp,missing)
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30

   ytitle = 'btd1112' & y = f.temp_11_0um_nom - f.temp_12_0um_nom & y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt11ratio_bt6711covar_logbt11std": begin
   classifier_name = 'bt11ratio_bt6711covar_logbt11std'

   xtitle = 'bt11ratio' 
   x = btratio(f.temp_11_0um_nom,f.surface_temperature_nwp,f.tropopause_temperature_nwp,missing)
   x_min = -0.1 & x_bin = 0.02 & nbins_x = 60

   ytitle = 'bt6711covar' & y = f.temp_11um_vs_67um_covar_5x5 & y_min = -2 & y_bin = 0.2 & nbins_y = 50.0


   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [6700,11000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"logzopa_topa_logbt11std": begin
   classifier_name = 'logzopa_topa_logbt11std'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'topa' & y = f.cld_temp_opaque & y_min = 180.0 & y_bin = 5.0 & nbins_y = 32
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"logzopa_dtsfcopa_logbt11std": begin
   classifier_name = 'logzopa_dtsfcopa_logbt11std'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'dtsfcopa' & y = f.surface_temperature_nwp - f.cld_temp_opaque & y_min = -10.0 & y_bin = 2.0 & nbins_y = 30
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"etropo11_dtsfcopa_logbt11std": begin
   classifier_name = 'etropo11_dtsfcopa_logbt11std'
   xtitle = 'etropo11' & x = f.emiss_tropo_11_0um_nom & x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'dtsfcopa' & y = f.surface_temperature_nwp - f.cld_temp_opaque & y_min = -10.0 & y_bin = 2.0 & nbins_y = 30
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"etropo11_topa_logbt11std": begin
   classifier_name = 'etropo11_topa_logbt11std'
   xtitle = 'etropo11' & x = f.emiss_tropo_11_0um_nom & x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'topa' & y = f.cld_temp_opaque & y_min = 180.0 & y_bin = 5.0 & nbins_y = 32
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"tsfc_btd1112_logbt11std": begin
   classifier_name = 'tsfc_btd1112_logbt11std'
   xtitle = 'tsfc' & x = f.surface_temperature_nwp & x_min = 200.0 & x_bin = 5.0 & nbins_x = 30
   ytitle = 'btd1112' & y = f.temp_11_0um_nom - f.temp_12_0um_nom & y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"tsfc_tpw_btd1112": begin
   classifier_name = 'tsfc_tpw_btd1112'
   xtitle = 'tsfc' & x = f.surface_temperature_nwp & x_min = 200.0 & x_bin = 5.0 & nbins_x = 30
   ytitle = 'tpw' & y = f.total_precipitable_water_nwp & y_min = 0.0 & y_bin = 0.2 & nbins_y = 40
   ztitle = 'btd1112' & z = f.temp_11_0um_nom - f.temp_12_0um_nom & z_min = -2.0 & z_bin = 0.2 & nbins_z = 50
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
end

"tsfc_tpw_btd11133": begin
   classifier_name = 'tsfc_tpw_btd11133'
   xtitle = 'tsfc' & x = f.surface_temperature_nwp & x_min = 200.0 & x_bin = 5.0 & nbins_x = 30
   ytitle = 'tpw' & y = f.total_precipitable_water_nwp & y_min = 0.0 & y_bin = 0.2 & nbins_y = 40
   ztitle = 'btd11133' & z = f.temp_11_0um_nom - f.temp_13_3um_nom & z_min = -10.0 & z_bin = 2.0 & nbins_z = 40
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
end


"logzopa_btd3811clr_logbt11std_day": begin
   classifier_name = 'logzopa_btd3811clr_logbt11std_day'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'dbtd3811clr'
   y = (f.temp_3_75um_nom - f.temp_11_0um_nom) - (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky)
   y_min = -10 & y_bin = 1.0 & nbins_y = 40.0
   ztitle = 'logbt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   tsfc_min = 230.0
   solzen_max=80.0
   coast_mask_max = 0
   zsfc_std_max = zsfc_std_max_iruni
end

"zopa_btd1112_bt11std": begin
   classifier_name = 'zopa_btd1112_bt11std'
   xtitle = 'zopa' & x = f.cld_height_opaque & x_min = 0.0 & x_bin = 250.0 & nbins_x = 60
   ytitle = 'btd1112' & y = f.temp_11_0um_nom - f.temp_12_0um_nom & y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'bt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = 0.0 & z_bin = 0.1 & nbins_z = 60.0
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   coast_mask_max = 0
   ;zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"zopa_btd3811clr_bt11std_day": begin
   classifier_name = 'zopa_btd3811clr_bt11std_day'
   xtitle = 'zopa' & x = f.cld_height_opaque & x_min = 0.0 & x_bin = 250.0 & nbins_x = 60
   ytitle = 'dbtd3811clr'
   y = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $
       (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky)
   y_min = -10 & y_bin = 1.0 & nbins_y = 40.0
   ztitle = 'bt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = 0.0 & z_bin = 0.1 & nbins_z = 60.0
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   tsfc_min = 230.0
   solzen_max=80.0
   coast_mask_max = 0
   ;zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end

"zopa_btd3811clr_bt11std_night": begin
   classifier_name = 'zopa_btd3811clr_bt11std_night'
   xtitle = 'zopa' & x = f.cld_height_opaque & x_min = 0.0 & x_bin = 250.0 & nbins_x = 60
   ytitle = 'dbtd3811clr'
   y = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $
       (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky)
   y_min = -10 & y_bin = 1.0 & nbins_y = 40.0
   ztitle = 'bt11std' & z = f.TEMP_11_0UM_NOM_STDDEV_3X3 & z_min = 0.0 & z_bin = 0.1 & nbins_z = 60.0
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   tsfc_min = 230.0
   solzen_min=90.0
   solzen_max=180.0
   coast_mask_max = 0
   ;zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end

"zopa_btd3811clr_logbt11std_day": begin
   classifier_name = 'zopa_btd3811clr_logbt11std_day'
   xtitle = 'zopa'
   x = f.cld_height_opaque
   x_min = 0.0 & x_bin = 500.0 & nbins_x = 30
   ytitle = 'dbtd3811clr'
   y = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $
       (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky)
   y_min = -10 & y_bin = 0.5 & nbins_y = 80.0
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   tsfc_min = 230.0
   solzen_max=80.0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end

"zopa_btd3811clr_logbt11std_night": begin
   classifier_name = 'zopa_btd3811clr_logbt11std_night'
   xtitle = 'zopa'
   x = f.cld_height_opaque
   x_min = 0.0 & x_bin = 500.0 & nbins_x = 30
   ytitle = 'dbtd3811clr'
   y = (f.temp_3_75um_nom - f.temp_11_0um_nom) - $
       (f.temp_3_75um_nom_clear_sky - f.temp_11_0um_nom_clear_sky)
   y_min = -10 & y_bin = 0.5 & nbins_y = 80.0
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   tsfc_min = 230.0
   solzen_min=90.0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end

"zopa_btd1112_btd8511": begin
   classifier_name = 'zopa_btd1112_btd8511'
   xtitle = 'zopa'
   x = f.cld_height_opaque
   x_min = 0.0 & x_bin = 500.0 & nbins_x = 30
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd8511'
   z = f.temp_8_5um_nom - f.temp_11_0um_nom
   z_min = -5 & z_bin = 0.5 & nbins_z = 30.0
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
end

"logzopa_btd1112_btd8511": begin
   classifier_name = 'logzopa_btd1112_btd8511'

   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])

   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd8511'
   z = f.temp_8_5um_nom - f.temp_11_0um_nom
   z_min = -5 & z_bin = 0.5 & nbins_z = 30.0
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
end

"bt11_btd1112_btd8511": begin
   classifier_name = 'bt11_btd1112_btd8511'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd8511'
   z = f.temp_8_5um_nom - f.temp_11_0um_nom
   z_min = -5 & z_bin = 0.5 & nbins_z = 30.0
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
end

"bt11_btd1112_btd3811_day": begin
   classifier_name = 'bt11_btd1112_btd3811_day'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd3811'
   z = f.temp_3_75um_nom - f.temp_11_0um_nom
   z_min = -5.0 & z_bin = 1.0 & nbins_z = 30
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
end

"bt11_btd1112_btd3811_night": begin
   classifier_name = 'bt11_btd1112_btd3811_night'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd3811'
   z = f.temp_3_75um_nom - f.temp_11_0um_nom
   z_min = -5.0 & z_bin = 1.0 & nbins_z = 30
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
   solzen_min = 90.0
end

"bt11_btd11133_btd3811_day": begin
   classifier_name = 'bt11_btd11133_btd3811_day'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd11133'
   y = f.temp_11_0um_nom - f.temp_13_3um_nom
   y_min = -10.0 & y_bin = 2.0 & nbins_y = 40
   ztitle = 'btd3811'
   z = f.temp_3_75um_nom - f.temp_11_0um_nom
   z_min = -5.0 & z_bin = 1.0 & nbins_z = 30
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
end

"bt11_btd11133_btd3811_night": begin
   classifier_name = 'bt11_btd11133_btd3811_night'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd11133'
   y = f.temp_11_0um_nom - f.temp_13_3um_nom
   y_min = -10.0 & y_bin = 2.0 & nbins_y = 40
   ztitle = 'btd3811'
   z = f.temp_3_75um_nom - f.temp_11_0um_nom
   z_min = -5.0 & z_bin = 1.0 & nbins_z = 30
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
   solzen_min = 90.0
end

"bt11_btd1112_logbt11std": begin
   classifier_name = 'bt11_btd1112_logbt11std'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   ztitle = 'logbt11std'
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [11000,12000]
   tsfc_min = 240.0
   zsfc_std_max = zsfc_std_max_iruni
end

"btd11clr_btd1112_btd8511": begin
   classifier_name = 'btd11clr_btd1112_btd8511'
   xtitle = 'btdclr11'
   x = f.temp_11_0um_nom_clear_sky - f.temp_11_0um_nom
   x_min = -15 & x_bin = 1.0 & nbins_x = 50.0
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd8511'
   z = f.temp_8_5um_nom - f.temp_11_0um_nom
   z_min = -5 & z_bin = 0.5 & nbins_z = 30.0
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
   tsfc_min = 240.0
end

"btd8511_btd1112_btd1110": begin
   classifier_name = 'btd8511_btd1112_btd1110'
   xtitle = 'btd8511'
   x = f.temp_8_5um_nom - f.temp_11_0um_nom
   x_min = -5 & x_bin = 0.5 & nbins_x = 30.0
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd1110'
   z = f.temp_11_0um_nom - f.temp_10_4um_nom
   z_min = -10 & z_bin = 0.5 & nbins_z = 40.0
   nchan_used = 3 & channel_wvl_used = [10400,11000,12000]
   tsfc_min = 240.0
end

"bt11_btd1112_btd1110": begin
   classifier_name = 'bt11_btd1112_btd1110'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd1112'
   y = f.temp_11_0um_nom - f.temp_12_0um_nom
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 50
   ztitle = 'btd1110'
   z = f.temp_11_0um_nom - f.temp_10_4um_nom
   z_min = -10 & z_bin = 0.5 & nbins_z = 40.0
   nchan_used = 3 & channel_wvl_used = [10400,11000,12000]
   tsfc_min = 240.0
end

"bt11_btd3811_logbt11std_all": begin
   classifier_name = 'bt11_btd3811_logbt11std_all'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd3811'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt11_btd3811_logbt11std_day": begin
   classifier_name = 'bt11_btd3811_logbt11std_day'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd3811'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt11_btd3811_logbt11std_night": begin
   classifier_name = 'bt11_btd3811_logbt11std_night'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd3811'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   solzen_min = 90.0
   solzen_max = 180.0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt10_btd3810_logbt10std_all": begin
   classifier_name = 'bt10_btd3810_logbt10std_all'
   xtitle = 'bt10'
   x = f.temp_10_4um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd3810'
   y = f.temp_3_75um_nom - f.temp_10_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt10std'
   z = f.TEMP_10_4UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,10400]
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt10_btd3810_logbt10std_day": begin
   classifier_name = 'bt10_btd3810_logbt10std_day'
   xtitle = 'bt10'
   x = f.temp_10_4um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd3810'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt10std'
   z = f.TEMP_10_4UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,10400]
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"bt10_btd3810_logbt10std_night": begin
   classifier_name = 'bt10_btd3810_logbt10std_night'
   xtitle = 'bt10'
   x = f.temp_10_4um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd3810'
   y = f.temp_3_75um_nom - f.temp_10_4um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt10std'
   z = f.TEMP_10_4UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,10400]
   solzen_min = 90.0
   solzen_max = 180.0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
end

"etropo11_btd8511_btd1112": begin
   classifier_name = 'etropo11_btd8511_btd1112'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'btd8511'
   y = f.temp_8_5um_nom - f.temp_11_0um_nom
   y_min = -5 & y_bin = 0.5 & nbins_y = 30.0
   ztitle = 'btd1112' & z = f.temp_11_0um_nom - f.temp_12_0um_nom & z_min = -2.0 & z_bin = 0.2 & nbins_z = 50
   nchan_used = 3 & channel_wvl_used = [8500,11000,12000]
end

"etropo11_btd8511_logbt11std": begin
   classifier_name = 'etropo11_btd8511_logbt11std'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'btd8511'
   y = f.temp_8_5um_nom - f.temp_11_0um_nom
   y_min = -5 & y_bin = 0.5 & nbins_y = 30.0
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [8500,10400]
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 230.0
end

"bt11_btd8511_logbt11std": begin
   classifier_name = 'bt11_btd8511_logbt11std'
   xtitle = 'bt11'
   x = f.temp_11_0um_nom
   x_min = 180.0 & x_bin = 5.0 & nbins_x = 32
   ytitle = 'btd8511'
   y = f.temp_8_5um_nom - f.temp_11_0um_nom
   y_min = -5 & y_bin = 0.5 & nbins_y = 30.0
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [8500,10400]
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
end


"etropo11_btd3811_logbt11std_all": begin
   classifier_name = 'etropo11_btd3811_logbt11std_all'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'btd3811'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,10400]
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
   tsfc_min = 230.0
end

"etropo11_btd3811_logbt11std_day": begin
   classifier_name = 'etropo11_btd3811_logbt11std_day'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'btd3811'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
   tsfc_min = 230.0
end

"etropo11_btd3811_logbt11std_night": begin
   classifier_name = 'etropo11_btd3811_logbt11std_night'
   xtitle = 'etropo11'
   x = f.emiss_tropo_11_0um_nom
   x_min = -0.4 & x_bin = 0.05 & nbins_x = 30
   ytitle = 'btd3811'
   y = f.temp_3_75um_nom - f.temp_11_0um_nom
   y_min = -5.0 & y_bin = 1.0 & nbins_y = 30
   ztitle = 'logbt11std'
   z = f.TEMP_11_0UM_NOM_STDDEV_3X3
   z_min = -1.5 & z_bin = 0.25 & nbins_z = 15.0
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc) & if (cc gt 0) then z[idx] = z_min
   nchan_used = 2 & channel_wvl_used = [3750,11000]
   solzen_min = 90.0
   solzen_max = 180.0
   coast_mask_max = 0
;  zsfc_max = 2000.0
   zsfc_std_max = zsfc_std_max_iruni
   tsfc_min = 240.0
   tsfc_min = 230.0
end

"logbt11std_logrefl065std_drefl065maxsubclr": begin
    classifier_name ='logbt11std_logrefl065std_drefl065maxsubclr'
    xtitle = 'logbt11std'
    x_min = -1.5 & x_bin = 0.10 & nbins_x = 36.0
    x = f.TEMP_11_0UM_NOM_STDDEV_3X3
    idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
    idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
    ytitle = 'logrefl065std'
    y_min = -1.5 & y_bin = 0.10 & nbins_y = 36.0
    y = f.REFL_0_65UM_NOM_STDDEV_3X3
    idx = where(y gt 0.0,cc) & if (cc gt 0) then y[idx] = alog10(y[idx])
    idx = where(y eq 0.0,cc)
    if (cc gt 0) then y[idx] = y_min
    ztitle = 'drefl065maxsubclr'
    z = f.refl_0_65um_nom_max_sub - f.refl_0_65um_nom_clear_sky
    z_min = -40.0 & z_bin = 2.0 & nbins_z = 60
    nchan_used = 2 & channel_wvl_used = [650,11000]
    coast_mask_min = 0
    coast_mask_max = 0
    solzen_min = 0.0
    solzen_max = solzen_day_max
    solscatang_min = solscatang_day_min
;   zsfc_max = 2000.0
    zsfc_std_max = zsfc_std_max_iruni
    snow_class_max = 1
end

"logbt11std_logrefl065std_drefl065clr": begin
    classifier_name ='logbt11std_logrefl065std_drefl065clr'
    xtitle = 'logbt11std'
    x_min = -1.5 & x_bin = 0.10 & nbins_x = 36.0
    x = f.TEMP_11_0UM_NOM_STDDEV_3X3
    idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
    idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
    ytitle = 'logrefl065std'
    y_min = -1.5 & y_bin = 0.10 & nbins_y = 36.0
    y = f.REFL_0_65UM_NOM_STDDEV_3X3
    idx = where(y gt 0.0,cc) & if (cc gt 0) then y[idx] = alog10(y[idx])
    idx = where(y eq 0.0,cc)
    if (cc gt 0) then y[idx] = y_min
    ztitle = 'drefl065clr'
    z = f.refl_0_65um_nom - f.refl_0_65um_nom_clear_sky
    z_min = -40.0 & z_bin = 2.0 & nbins_z = 60
    nchan_used = 2 & channel_wvl_used = [650,11000]
    coast_mask_min = 0
    coast_mask_max = 0
    solzen_min = 0.0
    solzen_max = solzen_day_max
    solscatang_min = solscatang_day_min
    ;zsfc_max = 2000.0
    zsfc_std_max = zsfc_std_max_iruni
    snow_class_max = 1
 
end

"refl065_refl160_refl138": begin
   classifier_name = 'refl065_refl160_refl138'
   xtitle = 'refl065' & x = f.refl_0_65um_nom & x_min = 0.0 & x_bin = 1.0 & nbins_x = 40.0
   ytitle = 'refl160' & y = f.refl_1_60um_nom & y_min = 0.0 & y_bin = 1.0 & nbins_y = 40.0
   ztitle = 'refl138' & z = f.refl_1_38um_nom & z_min = 0.0 & z_bin = 1.0 & nbins_z = 40.0
   nchan_used = 3 & channel_wvl_used = [650,1600,1380]
   tsfc_min = 180.0
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   slant_tpw_min = 0.5
   snow_class_max = 1
end

"refl065_refl160_logrefl138": begin
   classifier_name = 'refl065_refl160_logrefl138'
   xtitle = 'refl065' & x = f.refl_0_65um_nom & x_min = 0.0 & x_bin = 2.0 & nbins_x = 40.0
   ytitle = 'refl160' & y = f.refl_1_60um_nom & y_min = 0.0 & y_bin = 1.5 & nbins_y = 40.0
   ztitle = 'logrefl138' 
   z = f.refl_1_38um_nom 
   z_min = -2.0 & z_bin = 0.2 & nbins_z = 20
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc)
   if (cc gt 0) then z[idx] = z_min
   nchan_used = 3 & channel_wvl_used = [650,1600,1380]
   tsfc_min = 180.0
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   slant_tpw_min = 0.5
   snow_class_max = 1
end

"logrefl065_logrefl160_logrefl138": begin

   classifier_name = 'logrefl065_logrefl160_logrefl138'

   xtitle = 'logrefl065'
   x = f.refl_0_65um_nom 
   x_min = -2.0 & x_bin = 0.2 & nbins_x = 20
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc)
   if (cc gt 0) then x[idx] = x_min

   ytitle = 'logrefl160' 
   y = f.refl_1_60um_nom 
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 20
   idx = where(y gt 0.0,cc) & if (cc gt 0) then y[idx] = alog10(y[idx])
   idx = where(y eq 0.0,cc)
   if (cc gt 0) then y[idx] = y_min

   ztitle = 'logrefl138' 
   z = f.refl_1_38um_nom 
   z_min = -2.0 & z_bin = 0.2 & nbins_z = 20
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc)
   if (cc gt 0) then z[idx] = z_min

   nchan_used = 3 & channel_wvl_used = [650,1600,1380]
   tsfc_min = 180.0
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   slant_tpw_min = 0.5
   snow_class_max = 1
end

"logzopa_logrefl065_logrefl160": begin

   classifier_name = 'logzopa_logrefl065_logrefl160'

   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min

   ytitle = 'logrefl065'
   y = f.refl_0_65um_nom 
   y_min = -2.0 & y_bin = 0.2 & nbins_y = 20
   idx = where(y gt 0.0,cc) & if (cc gt 0) then y[idx] = alog10(y[idx])
   idx = where(y eq 0.0,cc)
   if (cc gt 0) then y[idx] = y_min

   ztitle = 'logrefl160' 
   z = f.refl_1_60um_nom 
   z_min = -2.0 & z_bin = 0.2 & nbins_z = 20
   idx = where(z gt 0.0,cc) & if (cc gt 0) then z[idx] = alog10(z[idx])
   idx = where(z eq 0.0,cc)
   if (cc gt 0) then z[idx] = z_min

   nchan_used = 3 & channel_wvl_used = [11000,650,1600]
   tsfc_min = 180.0
   solzen_min = 0.0
   solzen_max = solzen_day_max
   solscatang_min = solscatang_day_min
   solglint_mask_min = 0
   solglint_mask_max = 0
   snow_class_max = 1
end

"logzopa_btd1167_btd3811_night": begin
   classifier_name = 'logzopa_btd1167_btd3811_night'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'btd1167' & y = f.temp_11_0um_nom - f.temp_6_7um_nom & y_min = -5 & y_bin = 2.0 & nbins_y = 40.0
   ztitle = 'btd3811' & z = f.temp_3_75um_nom - f.temp_11_0um_nom & z_min = -10 & z_bin = 1.0 & nbins_z = 40
   nchan_used = 3 & channel_wvl_used = [3750,6700,1100]
   tsfc_min = 240.0
   solzen_min = 90.0
   solzen_max = 180.0
end

"logzopa_btd1167_btd3811_day": begin
   classifier_name = 'logzopa_btd1167_btd3811_day'
   xtitle = 'logzopa' & x = f.cld_height_opaque & x_min = 1.0 & x_bin = 0.1 & nbins_x = 30.0
   idx = where(x gt 0.0,cc) & if (cc gt 0) then x[idx] = alog10(x[idx])
   idx = where(x eq 0.0,cc) & if (cc gt 0) then x[idx] = x_min
   ytitle = 'btd1167' & y = f.temp_11_0um_nom - f.temp_6_7um_nom & y_min = -5 & y_bin = 2.0 & nbins_y = 40.0
   ztitle = 'btd3811' & z = f.temp_3_75um_nom - f.temp_11_0um_nom & z_min = -10 & z_bin = 1.0 & nbins_z = 40
   nchan_used = 3 & channel_wvl_used = [3750,6700,1100]
   tsfc_min = 240.0
   solzen_min = 0.0
   solzen_max = solzen_day_max
end

else: begin
      print, "error: unknown test"
      stop
      end

endcase

;------------------------------------------------
; make bin bounds
;------------------------------------------------
x_x = x_min + findgen(nbins_x)*x_bin + x_bin/2
x_y = y_min + findgen(nbins_y)*y_bin + y_bin/2
x_z = z_min + findgen(nbins_z)*z_bin + z_bin/2

;-----------------------------------------------
; make surface type
;----------------------------------------------
sfc_type_mask = f.bayes_mask_sfc_type

if (keyword_set(isccp)) then begin
   make_ecm_surface_type, f, sfc_type_mask_2, /isccp
endif

;-----------------------------------------------
;  apply filters
;----------------------------------------------
cod_lidar = f.CLOSEST_CALIPSO_COD
slant_tpw = f.total_precipitable_water_nwp / cos(f.sensor_zenith_angle * !DTOR)

if (keyword_set(new)) then begin
      phase_lidar = f.CLOSEST_CALIPSO_CLOUD_PHASE_TOP_LAYER
      frac_lidar = f.CLOSEST_CALIPSO_CLOUD_FRACTION_5KM
      tc_lidar = f.CLOSEST_CALIPSO_TOP_MID_TEMPERATURE
;     phase_lidar = f.CLOSEST_CALIPSO_PHASE_TOP_LAYER
;     frac_lidar = f.CLOSEST_CALIPSO_CLOUD_FRACTION
;     tc_lidar = f.CLOSEST_CALIPSO_MID_TEMPERATURE
      codi_lidar = f.CLOSEST_CALIPSO_COD_ICE
endif 
if (keyword_set(old)) then begin
      phase_lidar = f.CLOSEST_CALIPSO_CLOUD_PHASE_TOP_LAYER
      frac_lidar =f.CLOSEST_CALIPSO_CLOUD_FRACTION_x3
      codi_lidar = f.CLOSEST_CALIPSO_COD
      codi_lidar[*] = missing
endif

;--- convert tc_lidar from C to K
idx = where(tc_lidar ne  missing, cc)
if (cc gt 0) then tc_lidar[idx] = tc_lidar[idx] + 273.15

make_true_phase, phase_lidar, frac_lidar, tc_lidar, cod_lidar, codi_lidar,  $
                    cod_clear_water_thresh, cod_clear_ice_thresh, $
                    f.bayes_mask_sfc_type, f.cld_temp_opaque, f.cld_height_opaque, $
                    f.solar_zenith_angle, f.glint_mask, f.coast_mask, $
                    f.surface_elevation, f.surface_temperature_nwp, $
                    f.refl_0_65um_nom_min_3x3, f.refl_0_65um_nom, $
                    f.temp_11_0um_nom_max_3x3, f.temp_11_0um_nom, $
                    f.cld_opd_mask_0_65um_nom, cod_cloud_thresh, $
                    f.emiss_tropo_11_0um_nom, etrop_cloud_thresh, mod_flag, phase_true, missing

;-------------------------------------------------------------------------------
; fill in things that might be missing in training data
;-------------------------------------------------------------------------------
idx = where(tag_names(f) eq "CITY_MASK",cc)
if (cc eq 0) then begin
   city_mask_temp = make_array(ntrain,/FLOAT, VALUE=MISSING)
endif else begin
   city_mask_temp = f.city_mask
endelse

idx = where(tag_names(f) eq "MOON_ILLUM_FRAC",cc)
if (cc eq 0) then begin
   moon_illum_frac_temp = make_array(ntrain,/FLOAT, VALUE=MISSING)
endif else begin
   moon_illum_frac_temp = f.moon_illum_frac
endelse

apply_filters, phase_true, $
               f.CLOSEST_CALIPSO_TIME_DIFFERENCE, Time_Diff_Max, $
               f.sensor_zenith_angle, zen_min, zen_max, $
               f.solar_zenith_angle, solzen_min, solzen_max, $
               f.lunar_zenith_angle, lunzen_min, lunzen_max, $
               f.land_class,  $
	       frac_lidar, cfrac_min, cfrac_max,  $
               f.glint_zenith_angle, solglintzen_min, solglintzen_max,  $
               f.lunar_glint_zenith_angle, lunglintzen_min, lunglintzen_max, $
               f.scattering_angle, solscatang_min, solscatang_max,  $
               f.snow_class, snow_class_min, snow_class_max, $
               f.glint_mask, solglint_mask_min, solglint_mask_max, $
               f.glint_mask_lunar, lunglint_mask_min, lunglint_mask_max, $
               city_mask_temp, city_mask_min, city_mask_max, $
               moon_illum_frac_temp, moon_illum_frac_min, moon_illum_frac_max, $
               f.coast_mask, coast_mask_min, coast_mask_max, $
               slant_tpw, slant_tpw_min, slant_tpw_max, $
               f.surface_temperature_nwp, tsfc_min, tsfc_max, $
               f.surface_elevation, zsfc_min, zsfc_max, $
               f.surface_elevation_stddev_3x3, zsfc_std_min, zsfc_std_max, $
               ;f.latitude, lat_min, lat_max, missing, /phase_uni_flag
               f.latitude, lat_min, lat_max, missing
;-------------------------------
; select valid data within these limits
;-------------------------------
idx = where(phase_true ge 0.0 and x ne missing and y ne missing and z ne missing,cc)

print, 'total number after filtering= ', cc

;--- throw away filtered data
x = x[idx]
y = y[idx]
z = z[idx]
phase_true = phase_true[idx]
sfc_type_mask = sfc_type_mask[idx]

;-------------------------------------------------------------
; make arrays
;-------------------------------------------------------------
ice_table = make_array(nbins_x, nbins_y, nbins_z, nsfc, /FLOAT, VALUE=MISSING)
water_table = make_array(nbins_x, nbins_y, nbins_z, nsfc, /FLOAT, VALUE=MISSING)
clear_table = make_array(nbins_x, nbins_y, nbins_z, nsfc, /FLOAT, VALUE=MISSING)
obs_table = make_array(nbins_x, nbins_y, nbins_z, nsfc, /FLOAT, VALUE=0.0)
cloud_fraction = make_array(nsfc,/FLOAT,VALUE=missing)
clear_fraction = make_array(nsfc,/FLOAT,VALUE=missing)
ice_fraction = make_array(nsfc,/FLOAT,VALUE=missing)
water_fraction = make_array(nsfc,/FLOAT,VALUE=missing)
observation_count = make_array(nsfc,/FLOAT,VALUE=missing)

;--------------------------------------------------------------------
; loop over surface type
;--------------------------------------------------------------------
;for isfc = 0, 0 do begin             ;sfc_loop
for isfc = 0, nsfc - 1 do begin             ;sfc_loop

idx = where(sfc_type_mask eq isfc + 1,cc)
;print, 'isfc, count = ', isfc, cc
if (cc gt count_min) then begin
   x_sfc = x[idx]
   y_sfc = y[idx]
   z_sfc = z[idx]
   phase_true_sfc = phase_true[idx]
   observation_count[isfc] = cc
   ;print, 'isfc, count = ', isfc, cc
endif else begin
   goto, cycle_sfc
endelse

;---  make priors
idx = where(phase_true_sfc ge 0, count_all)
idx = where(phase_true_sfc eq 2, count_yes_ice_all)
idx = where(phase_true_sfc eq 1, count_yes_water_all)
idx = where(phase_true_sfc eq 0, count_yes_clear_all)

cloud_fraction[isfc] = float(count_yes_water_all + count_yes_ice_all) / float(count_all)
ice_fraction[isfc] = float(count_yes_ice_all) / float(count_all)
water_fraction[isfc] = float(count_yes_water_all) / float(count_all)
clear_fraction[isfc] = 1.0 - cloud_fraction[isfc]

cond_ratio_ice = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=MISSING)
cond_ratio_water = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=MISSING)
cond_ratio_clear = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=MISSING)
post_prob_ice = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=MISSING)
post_prob_water = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=MISSING)
post_prob_clear = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=MISSING)
obs_prob_lidar = make_array(nbins_x, nbins_y, nbins_z, /FLOAT, VALUE=0.0)

for i = 0, nbins_x - 1 do begin     ;i point_loop

    x_min_bin = x_min + i*x_bin 
    x_max_bin = x_min + (i+1)*x_bin 

for j = 0, nbins_y - 1 do begin     ;j-point loop

    y_min_bin = y_min + j*y_bin
    y_max_bin = y_min + (j+1)*y_bin

for k = 0, nbins_z - 1 do begin      ;k-point loop

    z_min_bin = z_min + k*z_bin
    z_max_bin = z_min + (k+1)*z_bin


    idx_all = where(x_sfc ge x_min_bin and x_sfc lt x_max_bin and $
                    y_sfc ge y_min_bin and y_sfc lt y_max_bin and $
                    z_sfc ge z_min_bin and z_sfc lt z_max_bin and $
                    phase_true_sfc ge 0, count_all_bin)

    if (count_all_bin lt count_min) then goto, skip_bin

    idx_ice = where(phase_true_sfc[idx_all] eq 2, count_yes_ice_bin)
    idx_water = where(phase_true_sfc[idx_all] eq 1, count_yes_water_bin)
    idx_clear = where(phase_true_sfc[idx_all] eq 0, count_yes_clear_bin)

    ;--- clear
    count_yes_bin = count_yes_clear_bin
    count_yes_all = count_yes_clear_all
    count_no_bin = count_all_bin - count_yes_clear_bin
    count_no_all = count_all - count_yes_clear_all
    prior_yes = clear_fraction[isfc] 
    compute_class_cond, nb_rank, laplace_flag, $
                        count_yes_all, count_no_all, count_yes_bin, count_no_bin, $
                        max_r, prior_yes, cond_yes, cond_no, r, posterior_prob
    cond_ratio_clear[i,j,k] = r
    post_prob_clear[i,j,k] = posterior_prob

    ;--- water
    count_yes_bin = count_yes_water_bin
    count_yes_all = count_yes_water_all
    count_no_bin = count_all_bin - count_yes_water_bin
    count_no_all = count_all - count_yes_water_all
    prior_yes = water_fraction[isfc] 
    compute_class_cond, nb_rank, laplace_flag, $
                        count_yes_all, count_no_all, count_yes_bin, count_no_bin, $
                        max_r, prior_yes, cond_yes, cond_no, r, posterior_prob
    cond_ratio_water[i,j,k] = r
    post_prob_water[i,j,k] = posterior_prob

    ;--- ice
    count_yes_bin = count_yes_ice_bin
    count_yes_all = count_yes_ice_all
    count_no_bin = count_all_bin - count_yes_ice_bin
    count_no_all = count_all - count_yes_ice_all
    prior_yes = ice_fraction[isfc] 
    compute_class_cond, nb_rank, laplace_flag, $
                        count_yes_all, count_no_all, count_yes_bin, count_no_bin, $
                        max_r, prior_yes, cond_yes, cond_no, r, posterior_prob
    cond_ratio_ice[i,j,k] = r
    post_prob_ice[i,j,k] = posterior_prob

    ;---- test
;   print, 'after ice ', r, posterior_prob
;   print, "counts ", count_all_bin, count_yes_clear_bin, count_yes_water_bin, count_yes_ice_bin, $
;                     count_yes_clear_bin + count_yes_water_bin + count_yes_ice_bin
;   print, "probs ", post_prob_clear[i,j,k], post_prob_water[i,j,k], post_prob_ice[i,j,k]
;   print, "prob sum ", post_prob_clear[i,j,k] + post_prob_water[i,j,k] + post_prob_ice[i,j,k]
    ;---- end test

    if (count_all gt count_min) then begin
     obs_prob_lidar[i,j,k] = float(count_all_bin)    
    endif

    skip_bin:

endfor    ;end of k-point loop
endfor    ;end of j-point loop
endfor    ;end of i-point loop

;--- normalize obs_prob
idx = where(obs_prob_lidar gt 0,cc)
if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar)

;-----------------------------------------------------------------------------------
; smooth
;-----------------------------------------------------------------------------------
;goto, skip_smooth
smooth_width = 3
cond_ratio_clear = smooth_nan(cond_ratio_clear,smooth_width,missing)
cond_ratio_water = smooth_nan(cond_ratio_water,smooth_width,missing)
cond_ratio_ice = smooth_nan(cond_ratio_ice,smooth_width,missing)
post_prob_clear = smooth_nan(post_prob_clear,smooth_width,missing)
post_prob_water = smooth_nan(post_prob_water,smooth_width,missing)
post_prob_ice = smooth_nan(post_prob_ice,smooth_width,missing)
;skip_smooth:

;--- normalize obs_prob
idx = where(obs_prob_lidar gt 0,cc)
if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar)

;-------------------
; fill in holes
;-------------------
;goto, skip_fill
obs_prob_interp_value = 1.0e-06
;idx = where(obs_prob_lidar eq 0,cc) & print, 'number with zero = ', cc
for ix = 0, nbins_x - 1 do begin
    ix_start = max([0,ix-1])
    ix_end = min([nbins_x-1,ix+1])
for iy = 0, nbins_y - 1 do begin
    iy_start = max([0,iy-1])
    iy_end = min([nbins_y-1,iy+1])
for iz = 0, nbins_z - 1 do begin
    iz_start = max([0,iz-1])
    iz_end = min([nbins_z-1,iz+1])

    if (obs_prob_lidar[ix,iy,iz] eq 0) then begin

       v = obs_prob_lidar[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
       idx = where(v gt 0,cc)

       if (cc gt 0) then begin

          obs_prob_lidar[ix,iy,iz] = mean(v[idx])

          v = cond_ratio_clear[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
          cond_ratio_clear[ix,iy,iz] = mean(v[idx])

          v = cond_ratio_water[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
          cond_ratio_water[ix,iy,iz] = mean(v[idx])

          v = cond_ratio_ice[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
          cond_ratio_ice[ix,iy,iz] = mean(v[idx])

          v = post_prob_clear[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
          post_prob_clear[ix,iy,iz] = mean(v[idx])

          v = post_prob_water[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
          post_prob_water[ix,iy,iz] = mean(v[idx])

          v = post_prob_ice[ix_start:ix_end,iy_start:iy_end,iz_start:iz_end]
          post_prob_ice[ix,iy,iz] = mean(v[idx])

       endif
    endif

endfor
endfor
endfor

;--- renormalize obs_prob
idx = where(obs_prob_lidar gt 0,cc)
if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar)
skip_fill:

;-----------------------------------------------------------------------------------
; smooth
;-----------------------------------------------------------------------------------
;goto, skip_smooth
smooth_width = 3
cond_ratio_clear = smooth_nan(cond_ratio_clear,smooth_width,missing)
cond_ratio_water = smooth_nan(cond_ratio_water,smooth_width,missing)
cond_ratio_ice = smooth_nan(cond_ratio_ice,smooth_width,missing)
post_prob_clear = smooth_nan(post_prob_clear,smooth_width,missing)
post_prob_water = smooth_nan(post_prob_water,smooth_width,missing)
post_prob_ice = smooth_nan(post_prob_ice,smooth_width,missing)
;skip_smooth:

;--- renormalize obs_prob
idx = where(obs_prob_lidar gt 0,cc)
if (cc gt 0) then obs_prob_lidar[idx] = obs_prob_lidar[idx] / total(obs_prob_lidar)

;------------------------------------------------------------------------------------
; apply table
;------------------------------------------------------------------------------------
n = n_elements(phase_true_sfc)
phase_retrieved = make_array(n, /FLOAT, VALUE = missing)
prob_cloud_retrieved = make_array(n, /FLOAT, VALUE = missing)
for l = 0,  n - 1 do begin

    i = min([nbins_x -1,max([0,fix((x_sfc[l] - x_min) / x_bin)])])
    j = min([nbins_y -1,max([0,fix((y_sfc[l] - y_min) / y_bin)])])
    k = min([nbins_z -1,max([0,fix((z_sfc[l] - z_min) / z_bin)])])

;   max_prob = max([post_prob_clear[i,j,k], post_prob_water[i,j,k] , post_prob_ice[i,j,k]])
;   if (post_prob_clear[i,j,k] eq max_prob) then phase_retrieved[l] = 0
;   if (post_prob_water[i,j,k] eq max_prob) then phase_retrieved[l] = 1
;   if (post_prob_ice[i,j,k] eq max_prob) then phase_retrieved[l] = 2

    phase_retrieved[l] = 0
    prob_cloud_retrieved[l] = 1.0 - post_prob_clear[i,j,k]
    if (post_prob_clear[i,j,k] lt 0.50) then begin
         phase_retrieved[l] = 1
         if (post_prob_ice[i,j,k] gt post_prob_water[i,j,k]) then begin
            phase_retrieved[l] = 2
         endif
    endif

endfor
;make_acc_metrics, phase_true_sfc, prob_cloud_retrieved, acc,bacc,wacc,bwacc,bce, tpr, tnr, fpr, fnr, phase_retrieved
;print, 'make_acc = ', acc,bacc,wacc,bwacc,bce, tpr, tnr, fpr, fnr

idx = where(phase_true_sfc ge 0 and phase_retrieved ge 0, count_all)
idx = where(phase_true_sfc eq 0 and phase_retrieved eq 0, count_agree_clear)
idx = where(phase_true_sfc gt 0 and phase_retrieved gt 0, count_agree_cloud)
idx = where(phase_true_sfc eq 1 and phase_retrieved eq 1, count_agree_water)
idx = where(phase_true_sfc eq 2 and phase_retrieved eq 2, count_agree_ice)
;idx = where(phase_true_sfc eq 0 and phase_retrieved ge 0, count_true_clear)
;idx = where(phase_true_sfc gt 0 and phase_retrieved ge 0, count_true_cloud)

if (count_all gt 0) then begin
  pod_clear = float(count_agree_cloud + count_agree_clear) / count_all
endif else begin
  pod_clear = missing
endelse

if (count_agree_cloud gt 0) then begin
  pod_phase = float(count_agree_water + count_agree_ice) / count_agree_cloud
endif else begin
  pod_phase = missing
endelse

print, "results for Isfc = ", isfc
print, 'pod clear = ', pod_clear
print, 'pod phase = ', pod_phase

;------------------------------------------------------------------------------------
; make images
;------------------------------------------------------------------------------------
;window, 0, xsize = 600, ysize = 600
;load_color_table, 33, 0, 1.0, 249
;erase, color = 0
;
;view, post_prob_clear, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $
;      sds_min=-101, sds_max=1, /mask_missing, title = 'CLEAR'
;
;plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $
;xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $
;yrange = [min(x_y),max(x_y)],xstyle = 1, $
;xtitle = xtitle, ytitle = ytitle
;legend_local, ['CLEAR'], color = 3, /top, /right, charsize = 2, /clear
;legend_local, [string(isfc+1,format="(I02)")], color = 3, /top, /left, charsize = 2, /clear
;saveimage, classifier_name+'_post_prob_clear_sfc'+string(isfc,format="(I01)")+'.png',/png
;
;
;window, 1, xsize = 600, ysize = 600
;load_color_table, 33, 0, 1.0, 249
;erase, color = 0
;
;view, post_prob_water, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $
;      sds_min=-101, sds_max=1, /mask_missing, title = 'WATER'
;
;plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $
;xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $
;yrange = [min(x_y),max(x_y)],xstyle = 1, $
;xtitle = xtitle, ytitle = ytitle
;legend_local, ['WATER'], color = 3, /top, /right, charsize = 2, /clear
;legend_local, [string(isfc+1,format="(I02)")], color = 3, /top, /left, charsize = 2, /clear
;saveimage, classifier_name+'_post_prob_water_sfc'+string(isfc,format="(I01)")+'.png',/png
;
;window, 2, xsize = 600, ysize = 600
;load_color_table, 33, 0, 1.0, 249
;erase, color = 0
;
;view, post_prob_ice, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $
;      sds_min=-101, sds_max=1, /mask_missing, title = 'ICE'
;
;plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $
;xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $
;yrange = [min(x_y),max(x_y)],xstyle = 1, $
;xtitle = xtitle, ytitle = ytitle
;legend_local, ['ICE'], color = 3, /top, /right, charsize = 2, /clear
;legend_local, [string(isfc+1,format="(I02)")], color = 3, /top, /left, charsize = 2, /clear
;saveimage, classifier_name+'_post_prob_ice_sfc'+string(isfc,format="(I01)")+'.png',/png
;
;window, 3, xsize = 600, ysize = 600
;load_color_table, 33, 0, 1.0, 249
;erase, color = 0
;
;view, obs_prob_lidar, coltab=33, /current, position = [0.10,0.06,0.98,0.98], cbar_position=cbar_position, /nocbar, missing_color=1, $
;      sds_min=-101, sds_max=max(obs_prob_lidar), /mask_missing, gam_fac=0.1
;
;plot, [0.0,0.0],[100.0,100.0],color = 3, background = 1 , $
;xrange = [min(x_x),max(x_x)],ystyle = 1, /nodata, /noerase, $
;yrange = [min(x_y),max(x_y)],xstyle = 1, $
;xtitle = xtitle, ytitle = ytitle
;legend_local, ['OBS_PROB'], color = 3, /top, /right, charsize = 2, /clear
;legend_local, [string(isfc+1,format="(I02)")], color = 3, /top, /left, charsize = 2, /clear
;saveimage, classifier_name+'_obs_prob_lidar_sfc'+string(isfc,format="(I01)")+'.png',/png
;
;
;--------------------------------------------
;
;--------------------------------------------
obs_table[*,*,*,isfc] = obs_prob_lidar
ice_table[*,*,*,isfc] = cond_ratio_ice
water_table[*,*,*,isfc] = cond_ratio_water
clear_table[*,*,*,isfc] = cond_ratio_clear

cycle_sfc:

endfor  ;sfc type loop

;-------------------------------------------------------
; constrain
;-------------------------------------------------------
idx = where(obs_table lt 0.0 and obs_table ne missing,cc) & if (cc gt 0) then obs_table[idx] = 0.00
idx = where(ice_table lt 0.0 and ice_table ne missing,cc) & if (cc gt 0) then ice_table[idx] = 0.00
idx = where(water_table lt 0.0 and water_table ne missing,cc) & if (cc gt 0) then water_table[idx] = 0.00
idx = where(clear_table lt 0.0 and clear_table ne missing,cc) & if (cc gt 0) then clear_table[idx] = 0.00

  ;-------------------------------------------------
  ; make a structure to hold data and attributes
  ;-------------------------------------------------

  ;--- assign attributes
  @attributes_define.pro
  attributes.classifier_name = classifier_name
  attributes.length_classifier_name = strlen(classifier_name)
  attributes.description = nc_header
  attributes.sensor = sensor_name
  attributes.training_source = training_source
  attributes.timestamp = timestamp()
  attributes.cod_clear_water_thresh = cod_clear_water_thresh
  attributes.cod_clear_ice_thresh = cod_clear_ice_thresh
  attributes.rank = nb_rank
  attributes.nsfc = nsfc
  attributes.nchan_used = nchan_used
  attributes.x_name = xtitle
  attributes.nbins_x = nbins_x
  attributes.x_min = x_min
  attributes.x_bin = x_bin
  attributes.y_name = ytitle
  attributes.nbins_y = nbins_y
  attributes.y_min = y_min
  attributes.y_bin = y_bin
  attributes.z_name = ztitle
  attributes.nbins_z = nbins_z
  attributes.z_min = z_min
  attributes.z_bin = z_bin
  attributes.lat_min = lat_min
  attributes.lat_max = lat_max
  attributes.zen_min = zen_min
  attributes.zen_max = zen_max
  attributes.solzen_min = solzen_min
  attributes.solzen_max = solzen_max
  attributes.solglintzen_min = solglintzen_min
  attributes.solglintzen_max = solglintzen_max
  attributes.solscatang_min = solscatang_min
  attributes.solscatang_max = solscatang_max
  attributes.lunzen_min = lunzen_min
  attributes.lunzen_max = lunzen_max
  attributes.lunglintzen_min = lunglintzen_min
  attributes.lunglintzen_max = lunglintzen_max
  attributes.lunscatang_min = lunscatang_min
  attributes.lunscatang_max = lunscatang_max
  attributes.tpw_min = slant_tpw_min
  attributes.tpw_max = slant_tpw_max
  attributes.tsfc_min = tsfc_min
  attributes.tsfc_max = tsfc_max
  attributes.zsfc_min = zsfc_min
  attributes.zsfc_max = zsfc_max
  attributes.zsfc_std_min = zsfc_std_min
  attributes.zsfc_std_max = zsfc_std_max
  attributes.snow_class_min = snow_class_min
  attributes.snow_class_max = snow_class_max
  attributes.solglint_mask_min = solglint_mask_min
  attributes.solglint_mask_max = solglint_mask_max
  attributes.lunglint_mask_min = lunglint_mask_min
  attributes.lunglint_mask_max = lunglint_mask_max
  attributes.coast_mask_min = coast_mask_min
  attributes.coast_mask_max = coast_mask_max
  attributes.city_mask_min = city_mask_min
  attributes.city_mask_max = city_mask_max
  attributes.moon_illum_frac_min = moon_illum_frac_min
  attributes.moon_illum_frac_max = moon_illum_frac_max
  attributes.rut_solzen_thresh = rut_solzen_thresh

  ;--- assign tables
  @lut_3d_define.pro
  lut.clear = clear_table
  lut.ice = ice_table
  lut.water = water_table
  lut.obs = obs_table
  lut.observation_count = observation_count
  lut.cloud_fraction = cloud_fraction
  lut.ice_fraction  = ice_fraction
  lut.water_fraction  = water_fraction
  lut.on_flag = sfc_type_on_flag
  channels_used.wvl = channel_wvl_used

;----------------------------------------------
; write to netcdf file
;----------------------------------------------
nc_filename = path_dir + sensor_name + "_" + day_string+"_"+classifier_name+".nc"

sd_id = ncdf_create(nc_filename,/CLOBBER,/NETCDF4_FORMAT)

write_3d_table, sd_id, lut, attributes, channels_used, success_flag, /class_att

ncdf_close, sd_id

print, 'write success = ', success_flag

end