;  Directly download granules
;https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/

;  Search for granules in a date time range
;https://search.earthdata.nasa.gov/search
;  In the search box put MYD06_L2 and select the dataset.
;  Enter time using calendar button.  Enter box using spacial rectangle SW=-76,58  NE=-49,152
;  Click download all
;  direct download
;  click download data
;  view/download data links
;  download links file


pro read_mod06_clouds_sb_test

;  imac or chpc
path_prefix='/uufs/chpc.utah.edu/common/home/'
;path_prefix='/Volumes/'

;  Output path
;output_path=''
;output_path=path_prefix+'mace-group4/modis/hysplit/modis_histograms_v3/'
output_path=path_prefix+'mace-group4/modis/hysplit/modis_histograms_sm/'

;  Modis file directory
fdir06=path_prefix+'mace-group4/modis/MYD06_L2/'
fdir03=path_prefix+'mace-group4/modis/MYD03/'

;  Loop through a range of day of the years
for n=60,90 do begin
;for n=4,4 do begin
  
  print,'n',n,'*************'
  nstr=string(n,format='(I03)')
  ;mod06_files=file_search(fdir06+'MYD06_L2.A2017*hdf',count=num_mod06)
  ;mod06_files=file_search(fdir06+'MYD06_L2.A2017'+nstr+'*hdf',count=num_mod06)
  ;mod06_files=file_search(fdir06+'MYD06_L2.A2017354*hdf',count=num_mod06)
  ;mod06_files=file_search(fdir06+'MYD06_L2.A2017345.0545.061.2018011134700*.hdf',count=num_mod06)
  mod06_files=file_search(fdir06+'MYD06_L2.A2018'+nstr+'*hdf',count=num_mod06)

  ;  Loop through the modis granules for that day of the year
  for f=0,num_mod06-1 do begin
    
    mod06_file=mod06_files[f]
    print, mod06_file

    ;  Pulls out the string  'MYD06_L2.A2018006.0625.061.2018006194906'
    output_file_string=strmid(file_basename(mod06_file),0,40)

    ;  Pick up the matching myd03 file
    parts=strsplit(file_basename(mod06_file),'.',/extract)
    syear=strmid(parts[1],1,4)
    sdoy=strmid(parts[1],5,3)
    mod03_file=file_search(fdir03+'MYD03.'+parts[1]+'.'+parts[2]+'.'+parts[3]+'*hdf',count=numfiles)
    mod03_file=mod03_file[0]    

    ;  If found a matching myd03 file    
    if numfiles eq 1 then begin
      print,'found ',mod03_file  

      ;  Output directory
      out_dir=output_path+syear+'/'+sdoy+'/'
      file_mkdir,out_dir
      
      ;  Read myd03 and myd06 hdf files
      file_id=hdf_sd_start(mod03_file)
      x_id=hdf_sd_nametoindex(file_id,'Latitude')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,lat_1km
      hdf_sd_endaccess,x_id2

      x_id=hdf_sd_nametoindex(file_id,'Longitude')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,lon_1km
      hdf_sd_endaccess,x_id2
      hdf_sd_end,file_id

      file_id=hdf_sd_start(mod06_file)
      x_id=hdf_sd_nametoindex(file_id,'Latitude')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,lat_5km
      hdf_sd_endaccess,x_id2

      x_id=hdf_sd_nametoindex(file_id,'Longitude')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,lon_5km
      hdf_sd_endaccess,x_id2

      x_id=hdf_sd_nametoindex(file_id,'Scan_Start_Time')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,scan_start_time
      hdf_sd_endaccess,x_id2

      day=julday(1,1,1993,0,0,0) & julian_day=day+(scan_start_time/86400.d)
      ;caldat, julian_day[0,0], mm, dd, yy, hh, mi, ss

      x_id=hdf_sd_nametoindex(file_id,'Solar_Zenith')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,solar_zenith_angle
      hdf_sd_endaccess,x_id2
      solar_zenith_angle=float(solar_zenith_angle)*0.01

      x_id=hdf_sd_nametoindex(file_id,'Sensor_Zenith')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,sensor_zenith_angle
      hdf_sd_endaccess,x_id2
      sensor_zenith_angle=float(sensor_zenith_angle)*0.01

      ;  scale_factor=0.01  offset=0  fillvalue=-9999 1KM
      x_id=hdf_sd_nametoindex(file_id,'Cloud_Effective_Radius')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,cloud_effective_radius2
      hdf_sd_endaccess,x_id2
      scaleid=hdf_sd_attrfind(x_id2,'scale_factor')
      hdf_sd_attrinfo,x_id2,scaleid,data=scale_factor  ;0.01
      scale_factor=scale_factor[0]
      cloud_effective_radius=float(cloud_effective_radius2)*scale_factor
      r=where(cloud_effective_radius2 eq -9999,c)
      if c gt 0 then cloud_effective_radius[r]=-9999

      ;  scale_factor=0.01  offset=-15000  fillvalue=-32768 5KM
      x_id=hdf_sd_nametoindex(file_id,'Cloud_Top_Temperature')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,cloud_top_temperature2
      hdf_sd_endaccess,x_id2
      scaleid=hdf_sd_attrfind(x_id2,'scale_factor')
      hdf_sd_attrinfo,x_id2,scaleid,data=scale_factor  ;0.01
      scale_factor=scale_factor[0]
      offsetid=hdf_sd_attrfind(x_id2,'add_offset')
      hdf_sd_attrinfo,x_id2,offsetid,data=offset  ;-15000.000 
      offset=offset[0]
      cloud_top_temperature=((float(cloud_top_temperature2)-offset)*scale_factor)  ;modis
      ;cloud_top_temperature=((float(cloud_top_temperature2))*scale_factor)+150.0  ;jay
      cloud_top_temperature=cloud_top_temperature-273.  ;convert to celcius
      r=where(cloud_top_temperature2 eq -32768,c)
      if c gt 0 then cloud_top_temperature[r]=-9999

      ;  scale_factor=0.01  offset=-15000  fillvalue=-32768
      x_id=hdf_sd_nametoindex(file_id,'Surface_Temperature')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,surface_temperature2
      hdf_sd_endaccess,x_id2
      scaleid=hdf_sd_attrfind(x_id2,'scale_factor')
      hdf_sd_attrinfo,x_id2,scaleid,data=scale_factor  ;0.01
      scale_factor=scale_factor[0]
      offsetid=hdf_sd_attrfind(x_id2,'add_offset')
      hdf_sd_attrinfo,x_id2,offsetid,data=offset  ;-15000.000
      offset=offset[0]
      surface_temperature=((float(surface_temperature2)-offset)*scale_factor)  ;modis
      ;surface_temperature=((float(surface_temperature2))*scale_factor)+150.  ;jay
      surface_temperature=surface_temperature-273.  ;convert to celcius

      ;  scale_factor=0.1  offset=0  fillvalue=-32768 5KM
      x_id=hdf_sd_nametoindex(file_id,'Cloud_Top_Pressure')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,cloud_top_pressure2
      hdf_sd_endaccess,x_id2
      scaleid=hdf_sd_attrfind(x_id2,'scale_factor')
      hdf_sd_attrinfo,x_id2,scaleid,data=scale_factor  ;0.1
      scale_factor=scale_factor[0]
      cloud_top_pressure=(float(cloud_top_pressure2)*scale_factor)*100. ; pa
      r=where(cloud_top_pressure2 eq -32768,c)
      if c gt 0 then cloud_top_pressure[r]=-9999

      ;  scale_factor=0.01  offset=0  fillvalue=-9999 1KM
      x_id=hdf_sd_nametoindex(file_id,'Cloud_Optical_Thickness')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,cloud_optical_thickness2
      hdf_sd_endaccess,x_id2
      scaleid=hdf_sd_attrfind(x_id2,'scale_factor')
      hdf_sd_attrinfo,x_id2,scaleid,data=scale_factor  ;0.01
      scale_factor=scale_factor[0]
      cloud_optical_thickness=float(cloud_optical_thickness2)*scale_factor
      r=where(cloud_optical_thickness2 eq -9999,c)
      if c gt 0 then cloud_optical_thickness[r]=-9999

      ;  scale_factor=1.0  offset=0  fillvalue=0  1KM
      x_id=hdf_sd_nametoindex(file_id,'Cloud_Phase_Optical_Properties')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2,cloud_phase
      hdf_sd_endaccess,x_id2

      ;  scale_factor=1  offset=0  fillvalue=-9999
      x_id=hdf_sd_nametoindex(file_id,'Cloud_Water_Path')
      x_id2=hdf_sd_select(file_id,x_id)
      hdf_sd_getdata,x_id2, cloud_water_path
      hdf_sd_endaccess,x_id2
      cloud_water_path=float(cloud_water_path)
      hdf_sd_end,file_id
      print,'finished reading data'

      ;  The size of the data arrays
      s=size(lat_1km,/dimensions)
      numx_1km=s[0]  ;lat_1km[*,0]
      numy_1km=s[1]  ;lat_1km[0,*]
      s=size(lat_5km,/dimensions)
      numx_5km=s[0]  ;lon_low_res[*,0]
      numy_5km=s[1]  ;lat_low_res[0,*]

      ;  Calculate colongitude
      colon_1km=lon_1km
      result=where(colon_1km lt 0,count)
      if count gt 0 then colon_1km[result]=360.0+colon_1km[result]
      colon_5km=lon_5km
      result=where(colon_5km lt 0,count)
      if count gt 0 then colon_5km[result]=360.0+colon_5km[result]

      ;****************************************************
      ;  Create an Nd array on the MODIS grid
      ;****************************************************

      print,'start nd array'

      ;  Increase resolution from 5km to 1km
      cloud_top_temp_1km=congrid(cloud_top_temperature,numx_1km,numy_1km)
      sensor_zenith_angle_1km=congrid(sensor_zenith_angle,numx_1km,numy_1km)
      solar_zenith_angle_1km=congrid(solar_zenith_angle,numx_1km,numy_1km)

      ;p0=contour(cloud_top_temp_1km,lon_1km,lat_1km,irregular=0)
      ;p1=contour(cloud_top_temperature,lon_5km,lat_5km,color='red',irregular=0)

      ;  These arrrays are calculated
      nd_array=make_array(/float,numx_1km,numy_1km,value=-9999)
      cw=make_array(/float,numx_1km,numy_1km,value=-9999)

      r=where(cloud_effective_radius gt 0 and cloud_water_path gt 0 and cloud_phase eq 2,count)
      if count gt 0 then begin
        cw[r]=0.6+((0.9/20.)*((cloud_top_temp_1km[r])+20.)) ; parameterized fro Wood's figure 1
        r2=where(cw lt 0.4 and cw ne -9999,c2)
        if c2 gt 0 then begin
          cw[r2]=0.4
        endif
        nd_array[r]=((sqrt(5.)/(2.*!pi*0.8))*sqrt((0.8*cw[r]*(1.e-6)*$
          cloud_optical_thickness[r])/(2.*1000.*(((1.e-6)*(cloud_effective_radius[r]))^5))))*1.e-6 ; Grosvenor et al 2018, eq 11 converted to 1/cm3
      endif
      print,'end nd array'

;****************************************************************
;  Carve out the microphysics for histogram writing in lat-lon regions.  
;  Only liquid phase cloud data with LWP < 250 g/m2 is used.  
;  Require view zenith to be less than 30 and solar zenith to be less than 60 
;****************************************************************

      ;  Check to see if there is any data that meets our conditions
      ;  0<re< 50, phase=2(liq), sensor zenith < 30 and solar zenith < 60
      good_data=make_array(/float,numx_1km,numy_1km,value=0)
      result=where(cloud_effective_radius gt 0. and cloud_effective_radius lt 50. and $
        cloud_phase eq 2 and sensor_zenith_angle_1km lt 30. and solar_zenith_angle_1km lt 60.,count_good)
      ;print,count_good,'count_good'

      ;  Continue on if there is good data
      if count_good gt 0 then begin
        good_data[result]=1  ;good data
        r=where(sensor_zenith_angle_1km gt 30.0,c)
        if c gt 0 then good_data[r]=2  ;mark these excluded values
        r=where(solar_zenith_angle_1km gt 60.0,c)
        if c gt 0 then good_data[r]=3        
        r=where(solar_zenith_angle_1km gt 60.0 and sensor_zenith_angle_1km gt 30,c)
        if c gt 0 then good_data[r]=4
        
        ;  Bins for histograms
        max_lwp=300. & dlwp=20.
        lwp_bins=0. & while max(lwp_bins) lt max_lwp do lwp_bins=[lwp_bins, max(lwp_bins)+dlwp] & histo_lwp=fltarr(n_elements(lwp_bins))
        max_re=30. & dre=2.
        re_bins=0. & while max(re_bins) lt max_re do re_bins=[re_bins, max(re_bins)+dre] & histo_re=fltarr(n_elements(re_bins))
        max_tau=50. & dtau=2.5
        tau_bins=0. & while max(tau_bins) lt max_tau do tau_bins=[tau_bins, max(tau_bins)+dtau] & histo_tau=fltarr(n_elements(tau_bins))
        max_solar_zenith=90. & dsolar_zenith=2.5
        solar_zenith_bins=30. & while max(solar_zenith_bins) lt max_solar_zenith do solar_zenith_bins=[solar_zenith_bins, max(solar_zenith_bins)+dsolar_zenith] & histo_szen=fltarr(n_elements(solar_zenith_bins))
        max_view_zenith=45. & dview_zenith=2.5
        view_zenith_bins=0. & while max(view_zenith_bins) lt max_view_zenith do view_zenith_bins=[view_zenith_bins, max(view_zenith_bins)+dview_zenith] & histo_vzen=fltarr(n_elements(view_zenith_bins))
        max_nd=300. & dnd=10
        nd_bins=0. & while max(nd_bins) lt max_nd do nd_bins=[nd_bins, max(nd_bins)+dnd] & histo_nd=fltarr(n_elements(nd_bins))
        max_temp=20.0 & dtemp=2.0
        temp_bins=-65.0 & while max(temp_bins) lt max_temp do temp_bins=[temp_bins, max(temp_bins)+dtemp] & histo_temp=fltarr(n_elements(temp_bins))
        max_phase=5 & dphase=1.0
        phase_bins=0 & while max(phase_bins) lt max_phase do phase_bins=[phase_bins, max(phase_bins)+dphase] & histo_phase=fltarr(n_elements(phase_bins))

        ;  Half bin width of the lat and lon grid
        ;dlon=2.5 & dlat=1.5  ;large grid
        dlon=1.0 & dlat=0.5  ;small grid
        
        ;  spaced lon array - turned this to a constant grid
        colon_vector_histo=(findgen((360.0/2*dlon))*2*dlon)+dlon

        ;  spaced lat array - turned this to a constant grid
        lat_vector_histo=((findgen(180.0/(2*dlat))*2*dlat)+dlat)-90.0        
        
        ;  Find the bounds of the granule and only loop through those 
        min_lat=floor(min(lat_1km)) & max_lat=ceil(max(lat_1km))
        min_colon=floor(min(colon_1km)) & max_colon=ceil(max(colon_1km))
        rlat=where(lat_vector_histo ge min_lat and lat_vector_histo le max_lat,clat)
        rcolon=where(colon_vector_histo ge min_colon and colon_vector_histo le max_colon,ccolon)
        ;print, min_lat,max_lat,min_colon,max_colon,'bounds of granule'
        ;  Step across the lat and lon 
        for jj=0,ccolon-1 do begin
          for kk=0,clat-1 do begin
            ship_latitude=lat_vector_histo[rlat[kk]]
            ship_longitude=colon_vector_histo[rcolon[jj]]
            ;print,ship_latitude,ship_longitude,'lat,lon,center box *********'
        
            ;  Find the number of points in the 1KM box
            r1km=where(abs(lat_1km-ship_latitude) lt dlat and abs(colon_1km-ship_longitude) lt dlon and $
              good_data eq 1,count_1km) 
            ;  Find the number of points in the 5KM box
            r5km=where(abs(lat_5km-ship_latitude) lt dlat and abs(colon_5km-ship_longitude) lt dlon,count_5km)

            ;  Calculate these mean values and do a temperature check
            if count_5km gt 0 then begin
              solar_zenith_int=solar_zenith_angle[r5km]
              sensor_zenith_int=sensor_zenith_angle[r5km]
              top_temp_int=cloud_top_temperature[r5km]
              rcold=where(top_temp_int lt -20 and top_temp_int ne -9999,ccold)
              if float(ccold)/float(count_5km) gt 0.10 then cold_flag=1 else cold_flag=0
            endif else begin
              solar_zenith_int=1e6
              sensor_zenith_int=1e6
              cold_flag=1e6
            endelse
           
            ;  Now see if there are enough points in the grid box  
            if count_5km gt 100 and count_1km gt 10 and $
              cold_flag eq 0 and $
              mean(sensor_zenith_int) lt 30. and mean(solar_zenith_int) lt 60. then begin
          
              ;  Subset the high res modis  
              ;print,ship_latitude,ship_longitude,count_1km,'  1KM lat and colon in box must be greater than 10'
              lat_vector_int=lat_1km[r1km]
              colon_vector_int=colon_1km[r1km]
              re_vector_int=cloud_effective_radius[r1km]
              lwp_vector_int=cloud_water_path[r1km]
              tau_vector_int=cloud_optical_thickness[r1km]
              nd_vector_int=nd_array[r1km]
              cloud_phase_int=cloud_phase[r1km]
              ;print,max(colon_vector_int),' ',min(colon_vector_int),' using these lons 1km'
              ;print,max(lat_vector_int),' ',min(lat_vector_int),' using these lats 1km'
      
              ;  subset the low res modis  
              solar_zenith_int=solar_zenith_angle[r5km]
              sensor_zenith_int=sensor_zenith_angle[r5km]
              colon_low_int=colon_5km[r5km]
              lat_low_int=lat_5km[r5km]
              julian_day_int=julian_day[r5km]
              cloud_top_temp_int=cloud_top_temperature[r5km]
              ;print,max(colon_low_int),' ',min(colon_low_int),' using these lons 5km'
              ;print,max(lat_low_int),' ',min(lat_low_int),' using these lats 5km'
              ;print,mean(sensor_zenith_int),mean(solar_zenith_int),' sensorZ<30,solarZ<60 5km'     

              ;  Reinitialize the histograms to zero
              histo_lwp=make_array(/float,n_elements(lwp_bins),value=0)
              histo_tau=make_array(/float,n_elements(tau_bins),value=0)
              histo_re=make_array(/float,n_elements(re_bins),value=0)
              histo_nd=make_array(/float,n_elements(nd_bins),value=0)
              histo_szen=make_array(/float,n_elements(solar_zenith_bins),value=0)
              histo_vzen=make_array(/float,n_elements(view_zenith_bins),value=0)
              histo_temp=make_array(/float,n_elements(temp_bins),value=0)
              histo_phase=make_array(/float,n_elements(phase_bins),value=0)  
                          
              ;  Create the histograms
              for j=0,n_elements(lwp_bins)-1 do begin
                d=where(abs(lwp_vector_int-lwp_bins[j]) le dlwp/2., count)
                histo_lwp[j]=count
              endfor
              for j=0,n_elements(tau_bins)-1 do begin
                d=where(abs(tau_vector_int-tau_bins[j]) le dtau/2., count)
                histo_tau[j]=count
              endfor
              for j=0,n_elements(re_bins)-1 do begin
                d=where(abs(re_vector_int-re_bins[j]) le dre/2., count)
                histo_re[j]=count
              endfor
              for j=0,n_elements(nd_bins)-1 do begin
                d=where(abs(nd_vector_int-nd_bins[j]) le dnd/2., count)
                histo_nd[j]=count
              endfor
              for j=0,n_elements(solar_zenith_bins)-1 do begin
                d=where(abs(solar_zenith_int-solar_zenith_bins[j]) le dsolar_zenith/2., count)
                histo_szen[j]=count
              endfor
              for j=0,n_elements(view_zenith_bins)-1 do begin
                d=where(abs(sensor_zenith_int-view_zenith_bins[j]) le dview_zenith/2., count)
                histo_vzen[j]=count
              endfor
              for j=0,n_elements(temp_bins)-1 do begin
                d=where(abs(cloud_top_temp_int-temp_bins[j]) le dtemp/2., count)
                histo_temp[j]=count
              endfor
              for j=0,n_elements(phase_bins)-1 do begin
                d=where(abs(cloud_phase_int-phase_bins[j]) le dphase/2., count)
                histo_phase[j]=count
              endfor

              distfreq = Histogram(cloud_phase_int, MIN=Min(cloud_phase_int))

              ; Find the maximum of the frequency distribution.
              maxfreq = Max(distfreq)

              ; Find the mode.
              cloud_phase_mode = Where(distfreq EQ maxfreq, count) + Min(cloud_phase_int)

              ; Find the 10th percentile value
              r=where(cloud_top_temp_int ne -9999)
              cloud_top_temp_int=cloud_top_temp_int[r]
              sidx=sort(cloud_top_temp_int)
              sort_cloud_top_temp_int=cloud_top_temp_int[sidx]
              p10=round(.10*n_elements(sort_cloud_top_temp_int))
              tenth_per_temp=sort_cloud_top_temp_int[p10]
              ;print,p10,n_elements(sort_cloud_top_temp_int),sort_cloud_top_temp_int[p10]
              
              r=where(cloud_top_temp_int ne -9999,c)
              if c gt 0 then begin
                temp_median=median(cloud_top_temp_int[where(cloud_top_temp_int ne -9999.0000)])              
              endif else begin
                temp_median=-9999
              endelse
              
              
              ; write out the histogram file
              lat_string=strtrim(string(fix(mean(lat_vector_int)), format='(i3)'),2)
              ;lon_string=strtrim(string(fix(mean(lon_vector_int)), format='(i4)'),2)
              lon_string=strtrim(string(fix(mean(colon_vector_int)), format='(i4)'),2)
              histogram_file_string=output_file_string+'_lat_'+lat_string+'_lon_'+lon_string+'_histo.cdf
              print,histogram_file_string
              cdfid=ncdf_create(out_dir+histogram_file_string, /clobber)
              nd_bins_did=ncdf_dimdef(cdfid,'nd',n_elements(nd_bins))
              re_bins_did=ncdf_dimdef(cdfid,'re',n_elements(re_bins))
              lwp_bins_did=ncdf_dimdef(cdfid,'lwp',n_elements(lwp_bins))
              szen_bins_did=ncdf_dimdef(cdfid,'solar_zenith',n_elements(solar_zenith_bins))
              vzen_bins_did=ncdf_dimdef(cdfid,'view_zenith',n_elements(view_zenith_bins))
              temp_bins_did=ncdf_dimdef(cdfid,'cloud_top_temp',n_elements(temp_bins))
              phase_bins_did=ncdf_dimdef(cdfid,'cloud_phase',n_elements(phase_bins))
                           
              julian_day_id=ncdf_vardef(cdfid,'julian_day',/double)

              center_lat_id=ncdf_vardef(cdfid,'center_latitude',/float)
              center_lon_id=ncdf_vardef(cdfid,'center_longitude',/float)
              
              count_5km_id=ncdf_vardef(cdfid,'count_5km',/float)
              count_1km_id=ncdf_vardef(cdfid,'count_1km',/float)
              
              mean_lat_id=ncdf_vardef(cdfid, 'mean_latitude', /float)
              mean_lon_id=ncdf_vardef(cdfid, 'mean_longitude', /float)
              mean_nd_id=ncdf_vardef(cdfid, 'mean_nd', /float)
              mean_re_id=ncdf_vardef(cdfid, 'mean_re', /float)
              mean_lwp_id=ncdf_vardef(cdfid, 'mean_lwp', /float)
              mean_temp_id=ncdf_vardef(cdfid, 'mean_cloud_top_temp', /float)
              median_temp_id=ncdf_vardef(cdfid, 'median_cloud_top_temp', /float)              
              tenth_temp_id=ncdf_vardef(cdfid, 'tenth_percentile_cloud_top_temp', /float)  
              mode_phase_id=ncdf_vardef(cdfid, 'mode_cloud_phase',/short)

              nd_bins_id=ncdf_vardef(cdfid, 'nd_bins', [nd_bins_did], /float)
              nd_histo_id=ncdf_vardef(cdfid, 'nd_histo', [nd_bins_did], /float)
              re_bins_id=ncdf_vardef(cdfid, 're_bins', [re_bins_did], /float)
              re_histo_id=ncdf_vardef(cdfid, 're_histo', [re_bins_did], /float)
              lwp_bins_id=ncdf_vardef(cdfid, 'lwp_bins', [lwp_bins_did], /float)
              lwp_histo_id=ncdf_vardef(cdfid, 'lwp_histo', [lwp_bins_did], /float)
              szen_bins_id=ncdf_vardef(cdfid, 'solar_zenith_bins', [szen_bins_did], /float)
              szen_histo_id=ncdf_vardef(cdfid, 'solar_zenith_histo', [szen_bins_did], /float)
              vzen_bins_id=ncdf_vardef(cdfid, 'view_zenith_bins', [vzen_bins_did], /float)
              vzen_histo_id=ncdf_vardef(cdfid, 'view_zenith_histo', [vzen_bins_did], /float)
              temp_bins_id=ncdf_vardef(cdfid, 'cloud_top_temp_bins', [temp_bins_did], /float)
              temp_histo_id=ncdf_vardef(cdfid, 'cloud_top_temp_histo', [temp_bins_did], /float)
              phase_bins_id=ncdf_vardef(cdfid, 'cloud_phase_bins', [phase_bins_did], /short)
              phase_histo_id=ncdf_vardef(cdfid, 'cloud_phase_histo', [phase_bins_did], /short)
  
              cold_flag_id=ncdf_vardef(cdfid,'cold_flag',/short)

              ncdf_control, cdfid, /endef

              ncdf_varput, cdfid, julian_day_id, mean(julian_day_int)
              ncdf_varput, cdfid, mean_lat_id, mean(lat_vector_int)
              ncdf_varput, cdfid, mean_lon_id, mean(colon_vector_int)
              ncdf_varput, cdfid, mean_nd_id, mean(nd_vector_int[where(nd_vector_int gt 0. and nd_vector_int lt 300. and re_vector_int gt 0. and lwp_vector_int lt 250.)])
              ncdf_varput, cdfid, mean_re_id, mean(re_vector_int[where(nd_vector_int gt 0. and nd_vector_int lt 300. and re_vector_int gt 0. and lwp_vector_int lt 250.)])
              ncdf_varput, cdfid, mean_lwp_id, mean(lwp_vector_int[where(nd_vector_int gt 0. and nd_vector_int lt 300. and re_vector_int gt 0. and lwp_vector_int lt 250.)])
              ncdf_varput, cdfid, mean_temp_id, mean(cloud_top_temp_int[where(cloud_top_temp_int ne -9999.0000)])
              ncdf_varput, cdfid, median_temp_id,temp_median
              ncdf_varput, cdfid, tenth_temp_id,tenth_per_temp

              ncdf_varput, cdfid, center_lat_id,ship_latitude
              ncdf_varput, cdfid, center_lon_id,ship_longitude
              ncdf_varput,cdfid,count_5km_id,count_5km
              ncdf_varput,cdfid,count_1km_id,count_1km
              ncdf_varput, cdfid, nd_bins_id, nd_bins
              ncdf_varput, cdfid, nd_histo_id, float(histo_nd)
              ncdf_varput, cdfid, re_bins_id, re_bins
              ncdf_varput, cdfid, re_histo_id, float(histo_re)
              ncdf_varput, cdfid, lwp_bins_id, lwp_bins
              ncdf_varput, cdfid, lwp_histo_id, float(histo_lwp)
              ncdf_varput, cdfid, szen_bins_id, solar_zenith_bins
              ncdf_varput, cdfid, szen_histo_id, float(histo_szen)
              ncdf_varput, cdfid, vzen_bins_id, view_zenith_bins
              ncdf_varput, cdfid, vzen_histo_id, float(histo_vzen)

              ncdf_varput, cdfid, temp_bins_id, temp_bins
              ncdf_varput, cdfid, temp_histo_id, float(histo_temp)
              ncdf_varput, cdfid, phase_bins_id, phase_bins
              ncdf_varput, cdfid, phase_histo_id, float(histo_phase)
              
              ncdf_varput,cdfid,cold_flag_id,cold_flag

              ncdf_close, cdfid
            endif ; count_5km gt 100 and count_1km gt 10 and sensor_zenith lt 30. and solar_zenith
          endfor  ; for k=0,n_elements(lat_vector_histo)-1 do begin
        endfor  ;for j=0,n_elements(lon_vector_histo)-1 do begin

        ;  Plot the histograms
        files=file_search(out_dir+output_file_string+'*_histo.cdf',count=num_hist_files)
        if 1 eq 0 then begin
        files=file_search(out_dir+output_file_string+'*_histo.cdf',count=num_hist_files)
        for i=0,num_hist_files-1 do begin
          print,files[i]
          fid=ncdf_open(files[i])
          vid=ncdf_varid(fid,'nd_bins') & ncdf_varget,fid,vid,nd_bins
          vid=ncdf_varid(fid,'nd_histo') & ncdf_varget,fid,vid,nd_histo
          vid=ncdf_varid(fid,'re_bins') & ncdf_varget,fid,vid,re_bins
          vid=ncdf_varid(fid,'re_histo') & ncdf_varget,fid,vid,re_histo
          vid=ncdf_varid(fid,'lwp_bins') & ncdf_varget,fid,vid,lwp_bins    
          vid=ncdf_varid(fid,'lwp_histo') & ncdf_varget,fid,vid,lwp_histo 
          vid=ncdf_varid(fid,'cloud_top_temp_bins') & ncdf_varget,fid,vid,temp_bins
          vid=ncdf_varid(fid,'cloud_top_temp_histo') & ncdf_varget,fid,vid,temp_histo          
          vid=ncdf_varid(fid,'solar_zenith_bins') & ncdf_varget,fid,vid,solar_zenith_bins
          vid=ncdf_varid(fid,'solar_zenith_histo') & ncdf_varget,fid,vid,solar_zenith_histo
          vid=ncdf_varid(fid,'view_zenith_bins') & ncdf_varget,fid,vid,view_zenith_bins
          vid=ncdf_varid(fid,'view_zenith_histo') & ncdf_varget,fid,vid,view_zenith_histo 
          vid=ncdf_varid(fid,'cloud_phase_bins') & ncdf_varget,fid,vid,phase_bins
          vid=ncdf_varid(fid,'cloud_phase_histo') & ncdf_varget,fid,vid,phase_histo           
          vid=ncdf_varid(fid,'mean_latitude') & ncdf_varget,fid,vid,mean_latitude  
          vid=ncdf_varid(fid,'mean_longitude') & ncdf_varget,fid,vid,mean_longitude     
          vid=ncdf_varid(fid,'mean_nd') & ncdf_varget,fid,vid,mean_nd     
          vid=ncdf_varid(fid,'mean_re') & ncdf_varget,fid,vid,mean_re    
          vid=ncdf_varid(fid,'mean_lwp') & ncdf_varget,fid,vid,mean_lwp   
          vid=ncdf_varid(fid,'mean_cloud_top_temp') & ncdf_varget,fid,vid,mean_temp   
                    
          ncdf_close,fid          
          pxdim=900
          pydim=800
          xl=0.08 & xr=0.90
          yb=0.07 & yt=0.95
          sx=0.10
          sy=0.08
          numplots_x=2
          numplots_y=3
          position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
          dx=0.01
          dy=0.01
          pnum=0
          p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],axis_style=4,/nodata)
          p1=plot(lwp_bins,lwp_histo/total(lwp_histo),/hist,position=pos[2,*],/current,/xstyle,$
            xtitle='LWP (g/m2)',ytitle='Frequency',yrange=[0,1],font_size=12)
          p1=plot(nd_bins,nd_histo/total(nd_histo),/hist,position=pos[0,*],/current,/xstyle,$
            xtitle='Nd (cm-3)',ytitle='Frequency',yrange=[0,1],font_size=12)  
          p1=plot(re_bins,re_histo/total(re_histo),/hist,position=pos[1,*],/current,/xstyle,$
            xtitle='Re (microns)',ytitle='Frequency',yrange=[0,1],font_size=12)
          p1=plot(solar_zenith_bins,solar_zenith_histo/total(solar_zenith_histo),/hist,$
            position=pos[3,*],/current,/xstyle,$
            xtitle='Solar Zenith',ytitle='Frequency',yrange=[0,1],font_size=12) 
          p1=plot(view_zenith_bins,view_zenith_histo/total(view_zenith_histo),/hist,$
            position=pos[4,*],/current,/xstyle,$
            xtitle='View Zenith',ytitle='Frequency',yrange=[0,1],font_size=12)  
          t1=text(pos[5,0],pos[5,3]-dy,string(mean_latitude),font_size=12)  
          t1=text(pos[5,0]+10*dx,pos[5,3]-dy,string(mean_longitude),font_size=12)  
          t1=text(pos[5,0]+3*dx,pos[5,3]-3*dy,'mean_nd',font_size=12)
          t1=text(pos[5,0]+13*dx,pos[5,3]-3*dy,'mean_re',font_size=12)
          t1=text(pos[5,0]+23*dx,pos[5,3]-3*dy,'mean_lwp',font_size=12)
          t1=text(pos[5,0],pos[5,3]-6*dy,string(mean_nd),font_size=12) 
          t1=text(pos[5,0]+10*dx,pos[5,3]-6*dy,string(mean_re),font_size=12)     
          t1=text(pos[5,0]+20*dx,pos[5,3]-6*dy,string(mean_lwp),font_size=12)               
          parts=strsplit(file_basename(files[i]),'cdf',/extract)  
          p0.save,out_dir+parts[0]+'png',height=pydim     
        endfor    ;  end of loop through histogram plot
        endif
    
        ;**********************************************
        ;  Regrid and plot
        ;**********************************************

        ;  This does the regridding and plotting
        if 1 eq 0 and num_hist_files gt 0 then begin
    
        ;*********************************************
        ;  Create a regular lat-lon grid at 1km
        ;*********************************************

        ;radius_earth=6357.d  ;km
        ;res=1. ;km
        ;dlat_1km=res*(180./(!dpi*radius_earth))
        ;dlon_1km=res*((180./(radius_earth*cos(mean(abs(lat_1km))*!dpi/180.d))))  
        ;lat_vector2=min(lat_1km) & while max(lat_vector2) lt max(lat_1km) do lat_vector2=[lat_vector2,max(lat_vector2)+dlat_1km]
        ;colon_vector2=min(colon_1km) & while max(colon_vector2) lt max(colon_1km) do colon_vector2=[colon_vector2,max(colon_vector2)+dlon_1km]

        dlat_1km=0.01
        dlon_1km=0.1
        lat_vector2=floor(min(lat_1km)) & while max(lat_vector2) lt max(lat_1km) do lat_vector2=[lat_vector2,max(lat_vector2)+dlat_1km]
        colon_vector2=floor(min(colon_1km)) & while max(colon_vector2) lt max(colon_1km) do colon_vector2=[colon_vector2,max(colon_vector2)+dlon_1km]

        lon_vector2=colon_vector2
        result=where(colon_vector2 gt 180,count)
        if count gt 0 then lon_vector2[result]=colon_vector2[result]-360.0

        plot_mod06_data,lat_5km,lon_5km,colon_5km,lat_1km,lon_1km,colon_1km,out_dir,$
          lat_vector2,lon_vector2,colon_vector2,lat_vector_histo,colon_vector_histo,$  
          solar_zenith_angle_1km,sensor_zenith_angle_1km,good_data,$
          cloud_top_temp_1km,cloud_top_temperature,cloud_optical_thickness,$
          cloud_effective_radius,cloud_water_path,cloud_phase,nd_array,cw,output_file_string,$
          plot_effective_radius,plot_nd_array,plot_liquid_water_path,plot_optical_depth,$  ;output regridded arrays
          plot_cloud_top_temp,plot_cloud_phase,plot_good_data,plot_sensor_zenith,plot_solar_zenith
    
        if 1 eq 1 then begin
        cdfid=ncdf_create(out_dir+output_file_string+'.cdf', /clobber)

        lat_did=ncdf_dimdef(cdfid, 'lat_did', n_elements(lat_vector2))
        lon_did=ncdf_dimdef(cdfid, 'lon_did', n_elements(lon_vector2))

        lat_vector_id=ncdf_vardef(cdfid, 'latitude_vector', [lat_did], /float)
        lon_vector_id=ncdf_vardef(cdfid, 'longitude_vector', [lon_did], /float)

        re_id=ncdf_vardef(cdfid, 're_interp_um', [lon_did, lat_did], /float)
        nd_id=ncdf_vardef(cdfid, 'nd_interp_per_cubic_cm', [lon_did, lat_did], /float)
        lwp_id=ncdf_vardef(cdfid, 'lwp_interp_g_per_sq_m', [lon_did, lat_did], /float)
        tau_id=ncdf_vardef(cdfid, 'tau_interp_unitless', [lon_did, lat_did], /float)
        ctt_id=ncdf_vardef(cdfid, 'cloud_top_temp_1km', [lon_did, lat_did], /float)
        cloud_phase_id=ncdf_vardef(cdfid,'cloud_phase_1km', [lon_did, lat_did], /float)
        good_data_id=ncdf_vardef(cdfid,'good_data_flag_1km', [lon_did, lat_did], /float)
        sensor_zenith_id=ncdf_vardef(cdfid,'sensor_zenith_1km', [lon_did, lat_did], /float)
        solar_zenith_id=ncdf_vardef(cdfid,'solar_zenith_1km', [lon_did, lat_did], /float)
  
        ncdf_control, cdfid, /endef
        ncdf_varput,cdfid,re_id,plot_effective_radius
        ncdf_varput,cdfid,nd_id,plot_nd_array
        ncdf_varput,cdfid,lwp_id,plot_liquid_water_path
        ncdf_varput,cdfid,tau_id,plot_optical_depth
        ncdf_varput,cdfid,lat_vector_id, lat_vector2
        ncdf_varput,cdfid,lon_vector_id, lon_vector2
        ncdf_varput,cdfid,ctt_id,plot_cloud_top_temp
        ncdf_varput,cdfid,cloud_phase_id,plot_cloud_phase
        ncdf_varput,cdfid,good_data_id,plot_good_data
        ncdf_varput,cdfid,sensor_zenith_id,plot_sensor_zenith
        ncdf_varput,cdfid,solar_zenith_id,plot_solar_zenith

        ncdf_close, cdfid
        endif  ;write regrid file
        endif  ;do plotting if there are histogram files
      endif else begin  ;found good data
        spawn,'rm '+mod06_file
        spawn,'rm '+mod03_file
        print,'no good data, files removed'
      endelse  ;remove files with no good data
    endif  ;found mod03
  endfor  ;end of loop through mod06 files
endfor  ;end of loop through doy
stop
end

pro plot_mod06_data,lat_5km,lon_5km,colon_5km,lat_1km,lon_1km,colon_1km,out_dir,$
    lat_vector2,lon_vector2,colon_vector2,lat_vector_histo,colon_vector_histo,$  
    solar_zenith_angle_1km,sensor_zenith_angle_1km,good_data,$
    cloud_top_temp_1km,cloud_top_temperature,cloud_optical_thickness,$
    cloud_effective_radius,cloud_water_path,cloud_phase,nd_array,cw,output_file_string,$
    plot_effective_radius,plot_nd_array,plot_liquid_water_path,plot_optical_depth,$  ;output regridded arrays
    plot_cloud_top_temp,plot_cloud_phase,plot_good_data,plot_sensor_zenith,plot_solar_zenith

print,'start plotting'

;  Top is the last color to scale 256 colors, 0-255
top_color=252
;  Colortable  0-252  253=white
;mytable=colortable(39,ncolors=254)
mytable=colortable(5,ncolors=254)
;254=hot pink               ;gray=255
mytable=[mytable,transpose([238,18,137]),transpose([180,180,180])]
mycbtable=mytable[0:top_color,*]

;***********************
;  Set up the plot
;***********************

pxdim=1000
pydim=1000

;  Set up the positions
xl=0.06 & xr=0.90
yb=0.05 & yt=0.95
sx=0.13
sy=0.05
numplots_x=2
numplots_y=2;num_files/numplots_x
position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
;  Colorbar position
cbpos=pos
cbpos[*,0]=pos[*,2]+0.07
cbpos[*,2]=cbpos[*,0]+0.007
cbpos[*,1]=cbpos[*,1]+0.02
cbpos[*,3]=cbpos[*,3]-0.02
dx=0.01
dy=0.01
;  Date format for plotting
dummy=label_date(date_format=['%Z/%M/%D %H:%I'])
pnum=0
p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],axis_style=4,/nodata)
d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device)
isx=(d[0,1,0]-d[0,0,0])
isy=(d[1,1,0]-d[1,0,0])

ulat=max(lat_1km)  ;upper lat of box
llat=min(lat_1km)  ;lower lat of box
llon=min(colon_1km)  ;left lon of box
rlon=max(colon_1km)  ;right lon of box

bulat=-50
bllat=-75
bllon=60
brlon=150
if bulat gt ulat then bulat=ulat  ;upper lat of box
if bllat lt llat then bllat=llat  ;lower lat of box
if bllon lt llon then bllon=llon  ;left lon of box
if brlon gt rlon then brlon=rlon  ;right lon of box

if (max(lon_1km) gt 179 and min(lon_1km) lt -179) or $
  (max(lat_1km) lt -60 ) or $
  (min(lat_1km) lt -75 ) then begin
  map_projection='Polar Stereographic'
  label_pos=1
  print,map_projection
endif else begin
  map_projection='Mercator'
  label_pos=0
endelse

;if mean(sensor_zenith_int) lt 30. and mean(solar_zenith_int) lt 60. then begin

;  Calculate regular grid array of the plotting box
;delta=0.01
;ynum=fix((ulat-llat)/delta)
;grid_lat=findgen(ynum)*delta+llat
;xnum=fix((rlon-llon)/delta)
;grid_lon=findgen(xnum)*delta+llon
grid_lat=lat_vector2
grid_lon=colon_vector2

;**************************
;  *** Start new plot ***
;**************************
pnum=0
p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],axis_style=4,/nodata)

;  CLOUD_WATER_PATH
if 1 eq 1 then begin
  print,'cloud water path'
  pnum=0
  data_var=cloud_water_path
  data_lat=lat_1km
  data_lon=colon_1km
  r=where(data_var gt 0,c)
  data_var[r]=alog10(data_var[r])
  r=where(data_var ne -9999,c)
  dmax_wp=max(data_var[r])
  dmin_wp=min(data_var[r])
  good_lat=data_lat[r]
  good_lon=data_lon[r]
  good_var=data_var[r]
  grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
  lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
  ;  Triangulate the data
  ;qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
  qhull,lon,lat,triangles,/delaunay,sphere=s
  ;triangulate,good_lon,good_lat,triangles,sphere=s
  ;grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
  grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
  ;;   Bytscal the data
  ;data_image=bytscl(grid_var,top=top_color,min=dmin_wp,max=dmax_wp)
  ;result=where(grid_var eq -9999,count)
  ;if count gt 0 then data_image[result]=255  ;gray
  ;p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
  ;  limit=[llat,llon,ulat,rlon],$  ;map limit vector
  ;  /current,$  ;Put it in the current open plot window p0
  ;  grid_units='degrees',$
  ;  center_longitude=((rlon-llon)/2.0)+llon,$
  ;  /box_axes,$
  ;  label_position=label_pos,$
  ;  font_size=12,$  ;bigger font
  ;  linestyle=2,$  ;grid will be dashed line
  ;  rgb_table=mytable,$  ;use color table 5
  ;  map_projection=map_projection,$  ;map projection
  ;  position=pos[pnum,*])
  ;p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
  ;p3=plot([bllon,brlon,brlon,bllon,bllon],[bllat,bllat,bulat,bulat,bllat],/overplot,/data,color='green',thick=3)
  ;;  Put on a colorbar
  ;c0=colorbar(range=[dmin_wp,dmax_wp],$  ;the colorbar is to match the image labeled p1
  ;  orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
  ;  tickdir=1,border_on=1,title='LOG Cloud Water Path 1km')
  ;t1=text(pos[pnum,0],pos[pnum,3]+dy,'Cloud Water Path',font_size=12)
  plot_liquid_water_path=grid_var
endif

;  CW 1km
if 1 eq 0 then begin
  print,'cw'
  data_var=cw
  data_lat=lat_1km
  data_lon=colon_1km
  r=where(data_var ne -9999,c)
  dmax_cw=max(data_var[r])
  dmin_cw=min(data_var[r])
  good_lat=data_lat[r]
  good_lon=data_lon[r]
  good_var=data_var[r]
  grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
  lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
  ;  Triangulate the data
  qhull,lon,lat,triangles,/delaunay,sphere=s
  grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
  ;   Bytscal the data
  data_image=bytscl(grid_var,top=top_color,min=dmin_cw,max=dmax_cw)
  result=where(grid_var eq -9999,count)
  if count gt 0 then data_image[result]=255  ;gray
  pnum=6
  p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
    limit=[llat,llon,ulat,rlon],$  ;map limit vector
    /current,$  ;Put it in the current open plot window p0
    grid_units='degrees',$
    /box_axes,$
    label_position=label_pos,$
    font_size=12,$  ;bigger font
    linestyle=2,$  ;grid will be dashed line
    rgb_table=mytable,$  ;use color table 5
    map_projection=map_projection,$  ;map projection
    position=pos[pnum,*])
  p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
  ;  Put on a colorbar
  c0=colorbar(range=[dmin_cw,dmax_cw],$  ;the colorbar is to match the image labeled p1
    orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
    tickdir=1,border_on=1,title='CW 1km')
endif

;  CLoud optical thickness 1km
if 1 eq 1 then begin
  print,'optical thickness'
  pnum=1
  data_var=cloud_optical_thickness
  data_lat=lat_1km
  data_lon=colon_1km
  r=where(data_var ne -9999,c)
  dmax_tau=50;max(data_var[r])
  dmin_tau=0;min(data_var[r])
  good_lat=data_lat[r]
  good_lon=data_lon[r]
  good_var=data_var[r]
  grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
  lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
  ;  Triangulate the data
  qhull,lon,lat,triangles,/delaunay,sphere=s
  grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
  ;   Bytscal the data
  data_image=bytscl(grid_var,top=top_color,min=dmin_tau,max=dmax_tau)
  result=where(grid_var eq -9999,count)
  if count gt 0 then data_image[result]=255  ;gray
  p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
    limit=[llat,llon,ulat,rlon],$  ;map limit vector
    /current,$  ;Put it in the current open plot window p0
    grid_units='degrees',$
    center_longitude=((rlon-llon)/2.0)+llon,$
    /box_axes,$
    label_position=label_pos,$
    font_size=12,$  ;bigger font
    linestyle=2,$  ;grid will be dashed line
    rgb_table=mytable,$  ;use color table 5
    map_projection=map_projection,$  ;map projection
    position=pos[pnum,*])
  p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
  ;  Put on a colorbar
  c0=colorbar(range=[dmin_tau,dmax_tau],$  ;the colorbar is to match the image labeled p1
    orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
    tickdir=1,border_on=1,title='Cloud Optical Thickness 1km')
  t1=text(pos[pnum,0],pos[pnum,3]+dy,'Cloud Optical Thickness 1km',font_size=12)
  plot_optical_depth=grid_var
endif

;  Solar Zenith Angle 1km
if 1 eq 1 then begin
  print,'solar zenith'
  pnum=2
  data_var=solar_zenith_angle_1km
  data_lat=lat_1km
  data_lon=colon_1km
  r=where(data_var ne -9999,c)
  dmax_solz=max(data_var[r])
  dmin_solz=min(data_var[r])
  good_lat=data_lat[r]
  good_lon=data_lon[r]
  good_var=data_var[r]
  grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
  lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
  ;  Triangulate the data
  qhull,lon,lat,triangles,/delaunay,sphere=s
  grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    ; min_points fewer, set to missing
    ;/kriging,min_points=8,sectors=8,empty_sectors=3,missing=-9999) ;original
    ;/kriging,min_points=4,sectors=8,empty_sectors=6,missing=-9999)  ;better
    /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)  ;much better 12,8,3
  ;   Bytscal the data
  data_image=bytscl(grid_var,top=top_color,min=dmin_solz,max=dmax_solz)
  result=where(grid_var eq -9999,count)
  if count gt 0 then data_image[result]=255  ;gray
  p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
    limit=[llat,llon,ulat,rlon],$  ;map limit vector
    /current,$  ;Put it in the current open plot window p0
    grid_units='degrees',$
    center_longitude=((rlon-llon)/2.0)+llon,$
    /box_axes,$
    label_position=label_pos,$
    font_size=12,$  ;bigger font
    linestyle=2,$  ;grid will be dashed line
    rgb_table=mytable,$  ;use color table 5
    map_projection=map_projection,$  ;map projection
    position=pos[pnum,*])
  p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
  ;  Put on a colorbar
  c0=colorbar(range=[dmin_solz,dmax_solz],$  ;the colorbar is to match the image labeled p1
    orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
    tickdir=1,border_on=1,title='Solar Zenith Angle 1km')
  t1=text(pos[pnum,0],pos[pnum,3]+dy,'Solar Zenith lt 60',font_size=12)
  plot_solar_zenith=grid_var
endif

;  Good Data 1km
if 1 eq 1 then begin
  print,'good data'
  pnum=3
  data_var=good_data
  data_lat=lat_1km
  data_lon=colon_1km
  r=where(data_var ne -9999,c)
  dmax_good=max(data_var[r])
  dmin_good=min(data_var[r])
  good_lat=data_lat[r]
  good_lon=data_lon[r]
  good_var=data_var[r]
  grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='First',epsilon=0.01
  lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
  ;  Triangulate the data
  qhull,lon,lat,triangles,/delaunay,sphere=s
  ;grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  ;  /grid,/degrees,/sphere,triangles=triangles,$
  ;  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
  grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    /nearest_neighbor,missing=-9999)
  r=where(plot_solar_zenith eq -9999,c)
  if c gt 0 then grid_var[r]=-9999
  ;   Bytscal the data
  data_image=bytscl(grid_var,top=top_color,min=dmin_good,max=dmax_good)
  result=where(grid_var eq -9999,count)
  if count gt 0 then data_image[result]=255  ;gray
  p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
    limit=[llat,llon,ulat,rlon],$  ;map limit vector
    /current,$  ;Put it in the current open plot window p0
    grid_units='degrees',$
    center_longitude=((rlon-llon)/2.0)+llon,$
    /box_axes,$
    label_position=label_pos,$
    font_size=12,$  ;bigger font
    linestyle=2,$  ;grid will be dashed line
    rgb_table=mytable,$  ;use color table 5
    map_projection=map_projection,$  ;map projection
    position=pos[pnum,*])
  p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
  ;  Put on a colorbar
  c0=colorbar(range=[dmin_good,dmax_good],$  ;the colorbar is to match the image labeled p1
    orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
    tickdir=1,border_on=1,title='good data 1km')
  t1=text(pos[pnum,0],pos[pnum,3]+dy,'Good data=solzen<60 senzen<30 re0-50 liq phase',font_size=12)
  plot_good_data=grid_var
endif

;  Sensor Zenith Angle 1km
if 1 eq 1 then begin
  print,'sensor zenith'
  pnum=0
  data_var=sensor_zenith_angle_1km
  data_lat=lat_1km
  data_lon=colon_1km
  r=where(data_var ne -9999,c)
  dmax=max(data_var[r])
  dmin=min(data_var[r])
  good_lat=data_lat[r]
  good_lon=data_lon[r]
  good_var=data_var[r]
  grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
  lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
  ;  Triangulate the data
  qhull,lon,lat,triangles,/delaunay,sphere=s
  grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
    /grid,/degrees,/sphere,triangles=triangles,$
    /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
  ;r=where(grid_var ge 0.5,c)
  ;if c gt 0 then grid_var[r]=1
  ;r=where(grid_var lt 0.5,c)
  ;if c gt 0 then grid_var[r]=0
  ;   Bytscal the data
  data_image=bytscl(grid_var,top=top_color,min=dmin,max=dmax)
  result=where(grid_var eq -9999,count)
  if count gt 0 then data_image[result]=255  ;gray
  p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
    limit=[llat,llon,ulat,rlon],$  ;map limit vector
    /current,$  ;Put it in the current open plot window p0
    grid_units='degrees',$
    center_longitude=((rlon-llon)/2.0)+llon,$
    /box_axes,$
    label_position=label_pos,$
    font_size=12,$  ;bigger font
    linestyle=2,$  ;grid will be dashed line
    rgb_table=mytable,$  ;use color table 5
    map_projection=map_projection,$  ;map projection
    position=pos[pnum,*])
  p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
  ;  Put on a colorbar
  c0=colorbar(range=[dmin,dmax],$  ;the colorbar is to match the image labeled p1
    orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
    tickdir=1,border_on=1,title='Sensor Zenith 1km')
  t1=text(pos[pnum,0],pos[pnum,3]+dy,'Sensor Zenith lt 30',font_size=12)
  plot_sensor_zenith=grid_var
endif


p0.save,out_dir+output_file_string+'.lwp_tau.png',height=pydim



;****************************
;  Start next image
;****************************
pnum=0
p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],axis_style=4,/nodata)

;  CLOUD_TOP_TEMP_1km 
print,'cloud top temp 1km'
pnum=0
data_var=cloud_top_temp_1km
data_var=data_var;[0:*:10,0:*:10]
data_lat=lat_1km;[0:*:10,0:*:10]
data_lon=colon_1km;[0:*:10,0:*:10]
r=where(data_var ne -9999,c)
dmax_temp=max(data_var[r])
dmin_temp=min(data_var[r])
good_lat=data_lat[r]
good_lon=data_lon[r]
good_var=data_var[r]
grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degrees,/sphere,duplicates='Avg',epsilon=0.01
lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
;  Triangulate the data
;qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
qhull,lon,lat,triangles,/delaunay,sphere=s
;grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
;  Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin_temp,max=dmax_temp)
result=where(grid_var eq -9999,count)
if count gt 0 then data_image[result]=255  ;gray
p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
  limit=[llat,llon,ulat,rlon],$  ;map limit vector
  /current,$  ;Put it in the current open plot window p0
  grid_units='degrees',$
  center_longitude=((rlon-llon)/2.0)+llon,$
  ;max_value=dmax_temp,min_value=dmin_temp,$
  /box_axes,$
  label_position=label_pos,$
  font_size=12,$  ;bigger font
  linestyle=2,$  ;grid will be dashed line
  rgb_table=mytable,$  ;use color table 5
  map_projection=map_projection,$  ;map projection
  position=pos[pnum,*])
p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
p3=plot([bllon,brlon,brlon,bllon,bllon],[bllat,bllat,bulat,bulat,bllat],/overplot,/data,color='green',thick=3)
;  Put on a colorbar
c0=colorbar(range=[dmin_temp,dmax_temp],$  ;the colorbar is to match the image labeled p1
  orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
  tickdir=1,border_on=1,title='Cloud Top Temp 1km')
t1=text(pos[pnum,0],pos[pnum,3]+dy,'Cloud Top Temp 1km',font_size=12)  
plot_cloud_top_temp=grid_var   

;  CLOUD_TOP_TEMPERATURE
if 1 eq 0 then begin
print,'cloud top temp 5km'
data_var=cloud_top_temperature
data_var=data_var;[0:*:2,0:*:2]
data_lat=lat_5km;[0:*:2,0:*:2]
data_lon=colon_5km;[0:*:2,0:*:2]
r=where(data_var ne -9999,c)
dmax_temp=max(data_var[r])
dmin_temp=min(data_var[r])
good_lat=data_lat[r]
good_lon=data_lon[r]
good_var=data_var[r]
;  Triangulate the data
qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
;  Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin_temp,max=dmax_temp)
result=where(grid_var eq -9999,count)
if count gt 0 then data_image[result]=255  ;gray
pnum=1
p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
  limit=[llat,llon,ulat,rlon],$  ;map limit vector
  /current,$  ;Put it in the current open plot window p0
  grid_units='degrees',$
  center_longitude=((rlon-llon)/2.0)+llon,$
  /box_axes,$
  label_position=label_pos,$
  font_size=12,$  ;bigger font
  linestyle=2,$  ;grid will be dashed line
  rgb_table=mytable,$  ;use color table 5
  map_projection=map_projection,$  ;map projection
  position=pos[pnum,*])
p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
;  Put on a colorbar
c0=colorbar(range=[dmin_temp,dmax_temp],$  ;the colorbar is to match the image labeled p1
  orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
  tickdir=1,border_on=1,title='Cloud Top Temp 5km')
endif  
  
  ;  CLOUD_EFFECTIVE_RADIUS
print,'cloud effective radius'
if 1 eq 1 then begin
pnum=1
data_var=cloud_effective_radius
data_var=data_var
data_lat=lat_1km
data_lon=colon_1km
r=where(data_var ne -9999,c)
dmax_re=30;max(data_var[r])
dmin_re=0;min(data_var[r])
good_lat=data_lat[r]
good_lon=data_lon[r]
good_var=data_var[r]
grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
;  Triangulate the data
;qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
qhull,lon,lat,triangles,/delaunay,sphere=s
;triangulate,good_lon,good_lat,triangles,sphere=s
;grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
;   Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin_re,max=dmax_re)
result=where(grid_var eq -9999,count)
if count gt 0 then data_image[result]=255  ;gray
p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
  limit=[llat,llon,ulat,rlon],$  ;map limit vector
  /current,$  ;Put it in the current open plot window p0
  grid_units='degrees',$
  center_longitude=((rlon-llon)/2.0)+llon,$
  /box_axes,$
  label_position=label_pos,$
  font_size=12,$  ;bigger font
  linestyle=2,$  ;grid will be dashed line
  rgb_table=mytable,$  ;use color table 5
  map_projection=map_projection,$  ;map projection
  position=pos[pnum,*])
p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
p3=plot([bllon,brlon,brlon,bllon,bllon],[bllat,bllat,bulat,bulat,bllat],/overplot,/data,color='green',thick=3)
;  Put on a colorbar
c0=colorbar(range=[dmin_re,dmax_re],$  ;the colorbar is to match the image labeled p1
  orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
  tickdir=1,border_on=1,title='Cloud Effective Radius 1km')
t1=text(pos[pnum,0],pos[pnum,3]+dy,'Cloud Effective Radius 1km',font_size=12)     
plot_effective_radius=grid_var  ;output regridded array
endif

;  ND_ARRAY
print,'nd array'
if 1 eq 1 then begin
pnum=2
data_var=nd_array
;r=where(data_var gt 0,c)
;data_var[r]=alog10(data_var[r])
data_lat=lat_1km
data_lon=colon_1km
r=where(data_var ne -9999,c)
dmax_nd=200;max(data_var[r])
dmin_nd=0;min(data_var[r])
good_lat=data_lat[r]
good_lon=data_lon[r]
good_var=data_var[r]
grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
;  Triangulate the data
;qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
qhull,lon,lat,triangles,/delaunay,sphere=s
;triangulate,good_lon,good_lat,triangles,sphere=s
;grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
;   Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin_nd,max=dmax_nd)
result=where(grid_var eq -9999,count)
if count gt 0 then data_image[result]=255  ;gray
p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
  limit=[llat,llon,ulat,rlon],$  ;map limit vector
  /current,$  ;Put it in the current open plot window p0
  grid_units='degrees',$
  center_longitude=((rlon-llon)/2.0)+llon,$
  /box_axes,$
  label_position=label_pos,$
  font_size=12,$  ;bigger font
  linestyle=2,$  ;grid will be dashed line
  rgb_table=mytable,$  ;use color table 5
  map_projection=map_projection,$  ;map projection
  position=pos[pnum,*])
p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
p3=plot([bllon,brlon,brlon,bllon,bllon],[bllat,bllat,bulat,bulat,bllat],/overplot,/data,color='green',thick=3)
;  Put on a colorbar
c0=colorbar(range=[dmin_nd,dmax_nd],$  ;the colorbar is to match the image labeled p1
  orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=mycbtable,$
  tickdir=1,border_on=1,title='Nd array 1km')
t1=text(pos[pnum,0],pos[pnum,3]+dy,'Nd Array 1km',font_size=12)     
plot_nd_array=grid_var
endif

;  CLOUD_PHASE 0=missing, 1=clear, 2=liquid, 3=ice, 4=undetermined
print,'cloud phase'
if 1 eq 1 then begin
pnum=3
data_var=cloud_phase
data_lat=lat_1km
data_lon=colon_1km
r=where(data_var ne -9999,c)
dmax_phase=max(data_var[r])
dmin_phase=min(data_var[r])
good_lat=data_lat[r]
good_lon=data_lon[r]
good_var=data_var[r]
grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='First',epsilon=0.01
lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
;  Triangulate the data
;qhull,good_lon,good_lat,triangles,/delaunay,sphere=s
qhull,lon,lat,triangles,/delaunay,sphere=s
;triangulate,good_lon,good_lat,triangles,sphere=s
;grid_var=griddata(good_lon,good_lat,good_var,xout=grid_lon,yout=grid_lat,$
;grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
;  /grid,/degrees,/sphere,triangles=triangles,$
;  /kriging,min_points=15,sectors=8,empty_sectors=3,missing=-9999)
grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /nearest_neighbor,missing=-9999)  
r=where(plot_solar_zenith eq -9999,c)
if c gt 0 then grid_var[r]=-9999
;   Bytscal the data
phase_table=['light gray','blue','green','gold','red']
data_image=byte(grid_var)
r=where(grid_var lt 0,c)
if c gt 0 then data_image[r]=0
p1=image(data_image,grid_lon,grid_lat,$  ;  I have called this one p1
  limit=[llat,llon,ulat,rlon],$  ;map limit vector
  /current,$  ;Put it in the current open plot window p0
  grid_units='degrees',$
  center_longitude=((rlon-llon)/2.0)+llon,$
  /box_axes,$
  label_position=label_pos,$
  font_size=12,$  ;bigger font
  linestyle=2,$  ;grid will be dashed line
  rgb_table=phase_table,$  ;use color table 5
  map_projection=map_projection,$  ;map projection
  position=pos[pnum,*])
p2=mapcontinents(color='white',/hires,thick=2)  ;add the continent lines
p3=plot([bllon,brlon,brlon,bllon,bllon],[bllat,bllat,bulat,bulat,bllat],/overplot,/data,color='black',thick=3)
;  Put on a colorbar
c0=colorbar(range=[dmin_phase,dmax_phase],$  ;the colorbar is to match the image labeled p1
  orientation=1,position=reform(cbpos[pnum,*]),font_size=12,rgb_table=phase_table,$
  tickdir=1,border_on=1,title='Cloud Phase 1km')
t1=text(pos[pnum,0],pos[pnum,3]+dy,'0=miss 1=clear 2=liq 3=ice 4=undet',font_size=14) 
plot_cloud_phase=grid_var
endif
p0.save,out_dir+output_file_string+'.temp_re_nd_phase.png',height=pydim


end