; =======================================================================================
  pro merge_files,yr,hem,version,revision,cv_exists,path_3c,path_3d,            $
                  n_add,rev_add,n_remove,rev_remove,fsum,fcv1,fcv2,             $
                  nhead_3c,head_3c,nhead_3d,head_3d, path_3e,head_3e,nhead_3e,  $ 
                  nstations,station_name,station_indx,lat0,lon0,coinc_type,     $
                  max_dist,lat_min,lat_max,lon_min,lon_max,file_unit_3e,fst1,fst2,nlog, $
                  def,tag_preliminary
; =======================================================================================
; Produces final files, either by merging new and old files or just removing 
; orbits from old files.
; ---------------------------------------------------------------------------------

; =======================
; Merge the 3c files
; =======================

;  --------------------------
;  Restore existing/old file
;  --------------------------
file=path_3c+filename_3c(hem,yr,version,revision,'','sav',tag_preliminary)
retrieve_cips_3c,file,nrev_old,nbin,nthresh,threshold,rev_old,date_old,dfs_old,lat_grid,mode,   $
                 lon_old,sza_old,ltime_old,ut_old,num_cld_old,num_obs_old,alb_old,      $
                 alb_std_old,rad_old,rad_std_old,iwc_old,iwc_std_old,alb_air_old,       $
                 alb_air_std_old,iwc_air_old,iwc_air_std_old, $
                 ndays_old,date_daily_old,dfs_daily_old,num_obs_daily_old,num_cld_daily_old,  $
                 alb_daily_old,iwc_daily_old,rad_daily_old,alb_air_daily_old,iwc_air_daily_old, $
                 iostatus

;  ----------------------------------
;  Restore tmp/new file if it exists
;  ----------------------------------
if(n_add ne 0) then begin
   file=path_3c+filename_3c(hem,yr,version,revision,'tmp','sav',tag_preliminary)
   retrieve_cips_3c,file,nrev_new,nbin,nthresh,threshold,rev_new,date_new,dfs_new,lat_grid,mode,   $
                    lon_new,sza_new,ltime_new,ut_new,num_cld_new,num_obs_new,alb_new,      $
                    alb_std_new,rad_new,rad_std_new,iwc_new,iwc_std_new,alb_air_new,       $
                    alb_air_std_new,iwc_air_new,iwc_air_std_new, $
                    ndays_new,date_daily_new,dfs_daily_new,num_obs_daily_new,num_cld_daily_new,  $
                    alb_daily_new,iwc_daily_new,rad_daily_new,alb_air_daily_new,iwc_air_daily_new, $
                    iostatus
endif

; --------------------------------------------------------------------------------
;  Get final output rev vector by combination of rev_old,
;  rev_add and rev_remove.
; --------------------------------------------------------------------------------

combine_rev_arrays,n_add,n_remove,rev_old,rev_new, $
                   rev_add,rev_remove,nrev_out,rev_out,indx,orig,nlog

;  ---------------------------
;  Generate new output file
;  ---------------------------

; First, handle the per-orbit data:

new=where(orig eq 'new',complement=old)
rev=rev_out
nrev=nrev_out

num_cld=intarr(nthresh,nrev,nbin)
num_obs=intarr(nthresh,nrev,nbin)
date=lonarr(nrev)
dfs=lonarr(nrev)
sza=fltarr(nthresh,nrev,nbin)
ut=fltarr(nthresh,nrev,nbin)
lon=fltarr(nthresh,nrev,nbin)
ltime=fltarr(nthresh,nrev,nbin)
alb=fltarr(nthresh,nrev,nbin)
alb_std=fltarr(nthresh,nrev,nbin)
iwc=fltarr(nthresh,nrev,nbin)
iwc_std=fltarr(nthresh,nrev,nbin)
rad=fltarr(nthresh,nrev,nbin)
rad_std=fltarr(nthresh,nrev,nbin)
alb_air=fltarr(nthresh,nrev,nbin)
alb_air_std=fltarr(nthresh,nrev,nbin)
iwc_air=fltarr(nthresh,nrev,nbin)
iwc_air_std=fltarr(nthresh,nrev,nbin)

date[old]=date_old[indx[old]]
dfs[old]=dfs_old[indx[old]]
num_cld[*,old,*]=num_cld_old[*,indx[old],*]
num_obs[*,old,*]=num_obs_old[*,indx[old],*]
sza[*,old,*]=sza_old[*,indx[old],*]
ut[*,old,*]=ut_old[*,indx[old],*]
lon[*,old,*]=lon_old[*,indx[old],*]
ltime[*,old,*]=ltime_old[*,indx[old],*]
alb[*,old,*]=alb_old[*,indx[old],*]
alb_std[*,old,*]=alb_std_old[*,indx[old],*]
iwc[*,old,*]=iwc_old[*,indx[old],*]
iwc_std[*,old,*]=iwc_std_old[*,indx[old],*]
rad[*,old,*]=rad_old[*,indx[old],*]
rad_std[*,old,*]=rad_std_old[*,indx[old],*]
alb_air[*,old,*]=alb_air_old[*,indx[old],*]
alb_air_std[*,old,*]=alb_air_std_old[*,indx[old],*]
iwc_air[*,old,*]=iwc_air_old[*,indx[old],*]
iwc_air_std[*,old,*]=iwc_air_std_old[*,indx[old],*]

date[new]=date_new[indx[new]]
dfs[new]=dfs_new[indx[new]]
num_cld[*,new,*]=num_cld_new[*,indx[new],*]
num_obs[*,new,*]=num_obs_new[*,indx[new],*]
sza[*,new,*]=sza_new[*,indx[new],*]
ut[*,new,*]=ut_new[*,indx[new],*]
lon[*,new,*]=lon_new[*,indx[new],*]
ltime[*,new,*]=ltime_new[*,indx[new],*]
alb[*,new,*]=alb_new[*,indx[new],*]
alb_std[*,new,*]=alb_std_new[*,indx[new],*]
iwc[*,new,*]=iwc_new[*,indx[new],*]
iwc_std[*,new,*]=iwc_std_new[*,indx[new],*]
rad[*,new,*]=rad_new[*,indx[new],*]
rad_std[*,new,*]=rad_std_new[*,indx[new],*]
alb_air[*,new,*]=alb_air_new[*,indx[new],*]
alb_air_std[*,new,*]=alb_air_std_new[*,indx[new],*]
iwc_air[*,new,*]=iwc_air_new[*,indx[new],*]
iwc_air_std[*,new,*]=iwc_air_std_new[*,indx[new],*]

; Now, since we have merged orbits data, we need to recalculate the
; daily-averaged cloud properties (it's kind of kludgy to do this here,
; rather than in create_level_3cde, but without significantly
; restructuring the code it's the most logical place:

calc_3c_daily_avg,nthresh,nrev,nbin,def,date,   $
                  num_obs,num_cld,alb,iwc,rad,alb_air,iwc_air,   $
                  ndays,date_daily,  $
                  num_obs_daily,num_cld_daily,frequency_daily,  $
                  alb_daily,iwc_daily,rad_daily,  $
                  alb_air_daily,iwc_air_daily

date_to_yymmdd_doy, date_daily, yy, mm, dd, day_of_year
dfs_from_doy,hem,fix(yr),day_of_year,udfs
dfs_daily = long(udfs)
;date_to_yymmdd_doy, date_daily, uniq_yy, uniq_mm, uniq_dd, uniq_doy
;   dfs_daily = date_daily*0L
;   for didx = 0, (ndays - 1) do begin
;      dfs_from_doy,hem,uniq_yy[didx],uniq_doy[didx],udfs
;      dfs_daily[didx] = udfs
;   end

final_3c_file,path_3c,version,revision,hem,fix(yr),threshold,nthresh,0,       $
   nrev,nbin,rev,date,dfs,lat_grid,num_cld,num_obs,lon,sza,ltime,ut,alb,rad,iwc,  $
   alb_std,rad_std,iwc_std,alb_air,iwc_air,alb_air_std,iwc_air_std,           $
   ndays,date_daily,dfs_daily,num_obs_daily,num_cld_daily,                 $
   alb_daily,iwc_daily,rad_daily,alb_air_daily,iwc_air_daily,tag_preliminary

; =======================
; Merge the 3d files
; =======================
if(cv_exists) then begin

;  --------------------
;  Read existing file
;  --------------------
   cvfile=filename_3d(hem,yr,version,revision,tag_preliminary)
   read_cips_3d,path_3d+cvfile,npix_old,nrev_old,rev_old,date_old,ut_old,ltime_old,ncv_old, $
                ncld_old,cld_presence_old,cld_frac_old,lat_old,lon_old,sza_old,rad_old,alb_old, $
                iwc_old,alb_air_old,iwc_air_old,dist_old,cld_map_old,qf_old,iostatus

;  ----------------------------
;  Read tmp file if it exists
;  ----------------------------
   if(n_add ne 0) then begin
     p0=strpos(cvfile,'.txt')
     tmpfile=strmid(cvfile,0,p0)+'_tmp.txt'
     read_cips_3d,path_3d+tmpfile,npix_new,nrev_new,rev_new,date_new,ut_new,ltime_new,ncv_new,  $
                ncld_new,cld_presence_new,cld_frac_new,lat_new,lon_new,sza_new,rad_new,    $
                alb_new,iwc_new,alb_air_new,iwc_air_new,dist_new,cld_map_new,qf_new,iostatus
   endif

;  ----------------------
;  Write output file
;  ----------------------
   newfile=path_3d+cvfile
   get_lun,unit
   openw,unit,newfile
   for ihead=0,nhead_3d-1 do printf,unit,head_3d[ihead]

   for irev=0,nrev_out-1 do begin   ; loop over orbits
      rev0=rev_out[irev]
      id=indx[irev]

      if(orig[irev] eq 'old') then begin
         if(ncv_old[id] eq 0) then begin
           printf,unit,format=fcv1,rev0,date_old[id],0.,0.,ncv_old[id],0,0,0.
         endif else begin
           printf,unit,format=fcv1,rev0,date_old[id],ut_old[id],ltime_old[id],  $
                              ncv_old[id],cld_presence_old[id],ncld_old[id],cld_frac_old[id]
           for m=0,ncv_old[id]-1 do begin
              printf,unit,format=fcv2,lat_old[id,m],lon_old[id,m],sza_old[id,m],dist_old[id,m],  $
                          rad_old[id,m],alb_old[id,m],iwc_old[id,m],alb_air_old[id,m],           $
                          iwc_air_old[id,m],qf_old[id,m],cld_map_old[id,m]
           endfor
         endelse
      endif else begin
         if(ncv_new[id] eq 0) then begin
           printf,unit,format=fcv1,rev0,date_new[id],0.,0.,ncv_new[id],0,0,0.
         endif else begin
           printf,unit,format=fcv1,rev0,date_new[id],ut_new[id],ltime_new[id],  $
                              ncv_new[id],cld_presence_new[id],ncld_new[id],cld_frac_new[id]
           for m=0,ncv_new[id]-1 do begin
              printf,unit,format=fcv2,lat_new[id,m],lon_new[id,m],sza_new[id,m],dist_new[id,m],  $
                          rad_new[id,m],alb_new[id,m],iwc_new[id,m],alb_air_new[id,m],           $
                          iwc_air_new[id,m],qf_new[id,m],cld_map_new[id,m]
           endfor
         endelse
      endelse
   endfor   ; loop over orbits
   close,unit
   free_lun,unit
endif

; =======================
; Merge the 3e files
; =======================
for ist=0,nstations-1 do close,file_unit_3e[ist]

for ist=0,nstations-1 do begin
   
;  --------------------
;  Read existing file
;  --------------------
   file_3e=filename_3e(hem,yr,version,revision,station_indx[ist],tag_preliminary)
   read_cips_3e,path_3e+file_3e,nmax_old,nrev_old,rev_old,date_old,ut_old,ltime_old, $
                npix_old,ncld_old,cld_presence_old,cld_frac_old,alb_loose_old,       $
                rad_loose_old,iwc_loose_old,alb_air_loose_old,iwc_air_loose_old,     $
                frac_loose_old,lat_old,lon_old,sza_old,                              $
                radius_old,albedo_old,iwc_old,albedo_air_old,iwc_air_old,dist_old,   $
                cld_map_old,qf_old,station_old,lat0_old,lon0_old,iostatus

;  ----------------------------
;  Read tmp file if it exists
;  ----------------------------
   if(n_add ne 0) then begin
     p0=strpos(file_3e,'.txt')
     tmpfile=strmid(file_3e,0,p0)+'_tmp.txt'
     read_cips_3e,path_3e+tmpfile,nmax_new,nrev_new,rev_new,date_new,ut_new,ltime_new, $
                  npix_new,ncld_new,cld_presence_new,cld_frac_new,alb_loose_new,       $
                  rad_loose_new,iwc_loose_new,alb_air_loose_new,iwc_air_loose_new,     $
                  frac_loose_new,lat_new,lon_new,sza_new,                              $
                  radius_new,albedo_new,iwc_new,albedo_air_new,iwc_air_new,dist_new,   $
                  cld_map_new,qf_new,station_new,lat0_new,lon0_new,iostatus
   endif

; ---------------------------------------------------------------------------
;  Need to re-define rev_add array to do the reconciliation. It's different 
;  for 3e files than 3c/3d
; ---------------------------------------------------------------------------
   combine_rev_arrays,nrev_new,n_remove,rev_old,rev_new,rev_new,rev_remove, $
                      nrev_out,rev_out,indx,orig,nlog

;  ----------------------
;  Write output file
;  ----------------------
   newfile=path_3e+file_3e
   unit=file_unit_3e[ist]
   openw,unit,newfile
   printf,unit,'#CIPS Coincidence Data for Station # '+strtrim(station_indx[ist],2)
   printf,unit,'#'+station_name[ist]
   printf,unit,'#'+strtrim(lat0[ist],2)+'  '+strtrim(lon0[ist],2)
   if(coinc_type[ist] eq 0) then begin
     printf,unit,'#Coincidence criteria type 0:   Max distance(km) = '+strtrim(fix(max_dist[ist]),2)
   endif else begin
     printf,unit,'#Coincidence criteria type 1: '  $
         +'min/max lat = '+ strtrim(lat_min[ist],2)+' '+strtrim(lat_max[ist],2) $
         +' ;  min/max lon = '+strtrim(lon_min[ist],2)+' '+strtrim(lon_max[ist],2)
   endelse
   for ihead=0,nhead_3e-1 do printf,unit,head_3e[ihead]

   for irev=0,nrev_out-1 do begin
      rev0=rev_out[irev]
      id=indx[irev]
      if(orig[irev] eq 'old') then begin
         printf,unit,format=fst1,rev0,date_old[id],ut_old[id],ltime_old[id],npix_old[id],  $
                     cld_presence_old[id],ncld_old[id],cld_frac_old[id],alb_loose_old[id], $
                     rad_loose_old[id],iwc_loose_old[id],alb_air_loose_old[id],            $
                     iwc_air_loose_old[id],frac_loose_old[id]
         for m=0,npix_old[id]-1 do begin
            printf,unit,format=fst2,lat_old[id,m],lon_old[id,m],sza_old[id,m],dist_old[id,m], $
                        radius_old[id,m],albedo_old[id,m],iwc_old[id,m],albedo_air_old[id,m], $
                        iwc_air_old[id,m],qf_old[id,m],cld_map_old[id,m]
         endfor
      endif else begin
         printf,unit,format=fst1,rev0,date_new[id],ut_new[id],ltime_new[id],npix_new[id],  $
                     cld_presence_new[id],ncld_new[id],cld_frac_new[id],alb_loose_new[id], $
                     rad_loose_new[id],iwc_loose_new[id],alb_air_loose_new[id],            $
                     iwc_air_loose_new[id],frac_loose_new[id]
         for m=0,npix_new[id]-1 do begin
            printf,unit,format=fst2,lat_new[id,m],lon_new[id,m],sza_new[id,m],dist_new[id,m], $
                        radius_new[id,m],albedo_new[id,m],iwc_new[id,m],albedo_air_new[id,m], $
                        iwc_air_new[id,m],qf_new[id,m],cld_map_new[id,m]
         endfor
      endelse

   endfor   ; loop over orbits
   close,unit
   free_lun,unit

endfor  ; loop over stations

return
end