; ================================================================================
  pro output_3e,unn,lat0,lon0,max_dist,latmin,latmax,lonmin,lonmax,            $
                latmin_loose,latmax_loose,lonmin_loose,lonmax_loose,           $
                coinc_type,cross_prime,max_dist_loose,hem,sza_min,sza_max,rev0,$
                date,rad_min,qf_in,ut_in,ltime_in,lat_in,lon_in,sza_in,alb_in, $
                rad_in,iwc_in,alb_air_in,iwc_air_in,cp_in,def,f1,f2
; ================================================================================
; This routine finds pixels coincident with the given station and 
; outputs results to the appropriate 3e file for one orbit.
; --------------------------------------------------------------------------------
; Modifications/updates:
;
;    May 2016  - Made changes to avoid I/O overflow for anomalous ALB/IWC values.
;  April 2018  - Modified to output the AIR albedo and IWC values.
; --------------------------------------------------------------------------------
min_coinc=100 ;key parameter - need more pixels than this for either reg or loose coinc.

; -------------------------------------------
; Copy input data into temporary arrays.
; -------------------------------------------
lat_station=lat0
lon_station=lon0
lat=lat_in
lon=lon_in
qf=qf_in
ut=ut_in
ltime=ltime_in
sza=sza_in
alb=alb_in
rad=rad_in
iwc=iwc_in
alb_air=alb_air_in
iwc_air=iwc_air_in
cp=cp_in

; If working in SH, Convert input L2 latitudes back to negative values.
if(hem eq 'south') then lat=-lat

; ------------------------------------------------------------------------
; If this station has been flagged as straddling the prime meridian, we
; need to convert CIPS longitude grid back to [-180,180]. Otherwise convert
; station lontitude arrays to [0,360] if necessary.
; ------------------------------------------------------------------------
if(cross_prime) then begin
  west=where(lon gt 180.,nwest)
  if(nwest ne 0) then lon(west)-=360.
endif else begin
  if(lon_station lt 0) then lon_station+=360.
  if(lonmin lt 0) then lonmin+=360.
  if(lonmax lt 0) then lonmax+=360.
  if(lonmin_loose lt 0) then lonmin_loose+=360.
  if(lonmax_loose lt 0) then lonmax_loose+=360.
endelse

; --------------------------------------------------------------------------------------
; Calculate average over loose coincidence range defined by max_dist_loose.
; First narrow things down a bit so we don't have to calculate distance for entire orbit.
; --------------------------------------------------------------------------------------
loose=where((sza ge sza_min) and (sza le sza_max) and           $
            (lat ge latmin_loose) and (lat le latmax_loose) and $
            (lon ge lonmin_loose) and (lon le lonmax_loose) and $
            (distance(lat_station,lon_station,lat,lon) le max_dist_loose),nloose)

alb_loose=0.
rad_loose=0.
iwc_loose=0.
alb_air_loose=0.
iwc_air_loose=0.
frac_loose=0.
min_cld=50
if(nloose eq 0) then goto,noloose
cld=where(cp[loose] gt 0,ncld)
if(ncld ge min_cld) then begin
  alb_loose=median(alb[loose[cld]])
  alb_air_loose=median(alb_air[loose[cld]])
  iwc_air_loose=median(iwc_air[loose[cld]])
  frac_loose=100.*float(ncld)/float(nloose)
endif
gd=where(cp[loose] gt 0 and rad[loose] ge rad_min,ngd)
if(ngd ge min_cld) then begin
  rad_loose=median(rad[loose[gd]])
  iwc_loose=median(iwc[loose[gd]])
endif else begin
  if(frac_loose gt 0.) then begin
    rad_loose=def
    iwc_loose=def
  endif
endelse
noloose:

; ------------------------------------------------------------------------
; Now find desired coincident points. If there are none then get out.
; ------------------------------------------------------------------------
ncoinc=0
if(coinc_type eq 0) then begin
; Work with a subset of "loose" criteria points here - 1st do sanity check. 
   if(max_dist ge max_dist_loose) then message,'distance parameters are screwed up!'
   if(nloose ne 0) then begin
      dist=distance(lat_station,lon_station,lat[loose],lon[loose])
      x=where(dist le max_dist,ncoinc)
      if(ncoinc ne 0) then coinc=loose[x]
   endif
endif else begin
  coinc=where(sza ge sza_min and sza le sza_max and $
        lat ge latmin and lat le latmax and lon ge lonmin and lon le lonmax,ncoinc)
endelse
if(ncoinc lt min_coinc) then return

; ---------------------------------------------
; Cut down arrays to just coincident points
; ---------------------------------------------
qf=qf[coinc]
ut=ut[coinc]
ltime=ltime[coinc]
lat=lat[coinc]
lon=lon[coinc]
sza=sza[coinc]
alb=alb[coinc]
rad=rad[coinc]
iwc=iwc[coinc]
alb_air=alb_air[coinc]
iwc_air=iwc_air[coinc]
cp=cp[coinc]
dist=distance(lat_station,lon_station,lat,lon)

; ----------------------------------------------------------------
; Convert longitude back to [0,360] if it was changed at top
; ----------------------------------------------------------------
if(cross_prime) then begin
  west=where(lon lt 0,nwest)
  if(nwest gt 0) then lon(west)+=360.
endif

; ----------------------------------------------------------------
; Calculate average UT and local time and cloud fraction.
; ----------------------------------------------------------------
meanut=mean(ut,/nan)
ltime=mean(ltime,/nan)

; --------------------------------
; Flag bad radius/IWC retrievals. 
; --------------------------------
bad=where(cp eq 1 and rad lt rad_min,nbad)
if(nbad gt 0) then begin
  rad[bad]=def
  iwc[bad]=def
endif

; ---------------------------------------------------------------------
; If there are NaN's in the RAD/IWC values, replace with def
; ---------------------------------------------------------------------
xxx=where(~finite(rad))
if(xxx[0] ne -1) then rad[xxx]=def
xxx=where(~finite(iwc))
if(xxx[0] ne -1) then iwc[xxx]=def

; ---------------------------------------------------------------------
; Next set all cloud parameters to 0 for non-cloud pixels. 
; ---------------------------------------------------------------------
xxx=where(cp eq 0)
if(xxx[0] ne -1) then begin
  alb[xxx] = 0.
  rad[xxx] = 0.
  iwc[xxx] = 0.
  alb_air[xxx] = 0.
  iwc_air[xxx] = 0.
endif

; ---------------------------------------------------------------------
; Search for high ALB/IWC values that will exceed format and replace
; ---------------------------------------------------------------------
maxval=99999.99
if(alb_loose gt maxval) then alb_loose=maxval
if(iwc_loose gt maxval) then iwc_loose=maxval
if(alb_air_loose gt maxval) then alb_air_loose=maxval
if(iwc_air_loose gt maxval) then iwc_air_loose=maxval
bad=where(alb gt maxval,nbad)
if(nbad ne 0) then alb[bad]=maxval
bad=where(iwc gt maxval,nbad)
if(nbad ne 0) then iwc[bad]=maxval
bad=where(alb_air gt maxval,nbad)
if(nbad ne 0) then alb_air[bad]=maxval
bad=where(iwc_air gt maxval,nbad)
if(nbad ne 0) then iwc_air[bad]=maxval

; --------------------------------
; Calculate cloud fraction.
; --------------------------------
cld=where(cp eq 1,ncld)
cld_frac=100.*float(ncld)/float(ncoinc)
if(ncld gt 0) then cld_presence=1 else cld_presence=0

; ----------------------
; Write output to file
; ----------------------

printf,unn,format=f1,rev0,date,meanut,ltime,ncoinc,cld_presence,ncld,cld_frac, $
                  alb_loose,rad_loose,iwc_loose,alb_air_loose,iwc_air_loose,frac_loose
for m=0,ncoinc-1 do begin
   printf,unn,format=f2,lat[m],lon[m],sza[m],dist[m],rad[m],alb[m],iwc[m],alb_air[m],iwc_air[m],qf[m],cp[m]
endfor

heap_gc
return
end