;Title:  acse_ingest_ceilo.pro
;Purpose:   Read the ceilometer data from ACSE and create netCDF files
;
;Inputs:	sarg:  date/time string, YYYYMMDDHH
;			earg:  data/time string, YYYYMMDDHH
;Outputs:	time:  decimal hour of day relative to the first day.
;			detstat:  ceilometer detection status: 0 - clear, 1-3 - number of cloud bases, 4 - use vertical visibility
;			cbase1:  1st cloud base height [m].  Available for detstat=[1,2,3]
;			cbase2:  2nd cloud base height [m].  Available for detstat=[1,2]
;			cbase3:  3rd cloud base height [m].  Available for detstat=3
;			vertvis: vertical visibility [m]. Available for detstat=4
;			highsig: altitude of highest signal [m].  Available for detstat=4
;Keywords:
;
;Notes:  	Ceilometer raw file naming convention:  CYMMDDHH.DAT
;			where C is a specified character, Y is the year
;			files are typically daily, unless the software has been restarted
;------------------------------------------------------------------------------
pro acse_ingest_ceilo

in_dir='/psd3data/cruises/ACSE_2014/ODEN/ceilometer/raw'
out_dir='/psd3data/cruises/ACSE_2014/ODEN/ceilometer/ingested'

sarg='20140704'
earg='20141003'

;Set up time reference
base_date=julday(1,1,1970,0,0,0)
base=jul_to_dt(base_date)
!dt_base=base


;Determine the files to open
sjd=julday(strmid(sarg,4,2),strmid(sarg,6,2),strmid(sarg,0,4))
if strmid(sarg,6,2) ne strmid(earg,6,2) then begin
	ejd=julday(strmid(earg,4,2),strmid(earg,6,2),strmid(earg,0,4))
	numa=ejd-sjd+1
	darg=strarr(numa)
	farg=darg
	for i=0,numa-1 do begin
		caldat,sjd+i,mo,da,yr
		if mo lt 10 then mo='0'+strtrim(mo,2) else mo=strtrim(mo,2)
		if da lt 10 then da='0'+strtrim(da,2) else da=strtrim(da,2)
		yr=strmid(strtrim(yr,2),3,1)
		farg[i]=yr+mo+da
		darg[i]='201'+farg[i]
	endfor
endif else begin
	darg=strmid(sarg,0,8)
	farg=strmid(sarg,3,5)
	numa=n_elements(darg)
endelse
sjd=sjd-julday(1,1,2014)+1

;Loop over all daily arguments
for d=0,numa-1 do begin

  ;Initialize variables
  ind=0L
  num=10000
  detstat=intarr(num)*0-1
  cbase1=intarr(num)*0
  cbase2=intarr(num)*0
  cbase3=intarr(num)*0
  vertvis=intarr(num)*0
  highsig=intarr(num)*0
  year=intarr(num)
  mon=intarr(num)
  day=intarr(num)
  time=fltarr(num)
  
  ;Open files, read, concatenate data
  jnk=' '
  str=' '
  dt1=' ' & dt2=' ' & dt3=' ' & dt4=' ' & dt5=' ' & dt6=' ' & dt7=' ' & dt8=' ' & dt9=' '

	;Get the list of files for the daily argument
	cd,in_dir
	arg='*'+farg[d]+'*.DAT'
	file=findfile(arg,count=numf)
	;Loop over all files in the day
	for f=0,numf-1 do begin
	  print,'Reading ',file[f]
		openr,lun,/get_lun,file[f]
		readf,lun,jnk
		readf,lun,jnk
		while not eof(lun) do begin
			readf,lun,str,format="(a20)"
			if strmid(str,0,2) eq '-F' then goto,flag
			if strmid(str,0,2) ne '-2' then begin
			   print,'off'
			   stop
			endif
			readf,lun,jnk
			readf,lun,dt1,dt2,dt3,dt4,dt5,dt6,dt7,dt8,dt9,format="(a2,a1,a5,a1,a5,a1,a5,a1,a12)"
			
			for i=0,4 do readf,lun,jnk
			;Store data
			detstat[ind]=fix(strmid(dt1,0,1))
			case detstat[ind] of
				1: cbase1[ind]=fix(dt3)
				2: begin
					cbase1[ind]=fix(dt3)
					cbase2[ind]=fix(dt5)
				end
				3: begin
					cbase1[ind]=fix(dt3)
					cbase2[ind]=fix(dt5)
					cbase3[ind]=fix(dt7)
				end
				4: begin
					vertvis[ind]=fix(dt3)
					highsig[ind]=fix(dt5)
				end
				else:
			endcase
			time[ind]=float(strmid(str,12,2))+float(strmid(str,15,2))/60.+float(strmid(str,18,2))/3600.
			ind=ind+1
		endwhile
		flag:
		free_lun,lun
	endfor

  if ind gt 0 then begin
  	detstat=detstat[0:ind-1]
  	cbase1=float(cbase1[0:ind-1])
  	cbase2=float(cbase2[0:ind-1])
  	cbase3=float(cbase3[0:ind-1])
  	highsig=float(highsig[0:ind-1])
  	vertvis=float(vertvis[0:ind-1])
  	time=time[0:ind-1]
    ;compute base_time and time_offset
    sec=dt_to_sec( str_to_dt( strmid(darg[d],4,2)+'/'+strmid(darg[d],6,2)+'/'+strmid(darg[d],0,4) ) )
    base_time=sec+fix(time[0]*3600)
    time_offset=long((time-time[0])*3600)
    
    ;Write output
    cd,out_dir
    hr=fix(time[0])
    if hr lt 10 then hr='0'+strtrim(hr,2) else hr=strtrim(hr,2)
    mi=fix((time[0] mod 1)*60.)
    if mi lt 10 then mi='0'+strtrim(mi,2) else mi=strtrim(mi,2)
    se=fix((((time[0] mod 1)*60.) mod 1)*60.)
    if se lt 10 then se='0'+strtrim(se,2) else se=strtrim(se,2)
    basetimestring=strmid(darg[d],0,4)+'-'+strmid(darg[d],4,2)+'-'+strmid(darg[d],6,2)+' '+hr+':'+mi+':'+se+' 0:00'

    outfilename='acsvceilcbC1.b1.'+darg[d]+'.'+hr+mi+se
    nbeam=n_elements(time)
    nheights=n_elements(height)
    outfile=outfilename+'.cdf'
    print,'WRITING: ',outfile
    outfid=ncdf_create(outfile,/clobber)
    ;Define Dimensions
    ;----------------------
    tdid=ncdf_dimdef(outfid,'time',nbeam)
    hdid=ncdf_dimdef(outfid,'range',nheights)
    sdid=ncdf_dimdef(outfid,'string_len',8)
    ;Define Variables
    ;-----------------------
    btid=ncdf_vardef(outfid,'base_time',/long)
    ncdf_attput,outfid,btid,'string',basetimestring
    ncdf_attput,outfid,btid,'long_name','Base time in Epoch'
    ncdf_attput,outfid,btid,'units','seconds since 1970-1-1 0:00:00 0:00'

    toid=ncdf_vardef(outfid,'time_offset',[tdid],/double)
    ncdf_attput,outfid,toid,'long_name','Time offset from base_time'
    ncdf_attput,outfid,toid,'units','seconds since '+basetimestring

    tid=ncdf_vardef(outfid,'time',[tdid],/double)   ;**** This is different from ARM... they have seconds since start of day
    ncdf_attput,outfid,tid,'long_name','Time'
    ncdf_attput,outfid,tid,'units','Decimal hours'

;    rid=ncdf_vardef(outfid,'range',[hdid],/float)
;    ncdf_attput,outfid,rid,'long_name','Distance to center of range bin'
;    ncdf_attput,outfid,rid,'units','m, AGL'
;    ncdf_attput,outfid,rid,'comment','none'

    dsid=ncdf_vardef(outfid,'detection_status',[tdid],/long)
    ncdf_attput,outfid,dsid,'long_name','Detection status'
    ncdf_attput,outfid,dsid,'units','unitless'
    ncdf_attput,outfid,dsid,'values','0=No significant backscatter; 1-3=Number of cloud bases detected; 4=Full obscuration determined but no cloud base detected; 5=Some obscuration detected but determined to be transparent.'

;    sfid=ncdf_vardef(outfid,'status_flag',[tdid],/char)
;    ncdf_attput,outfid,sfid,'long_name','Ceilometer status indicator'
;    ncdf_attput,outfid,sfid,'units','unitless'
;    ncdf_attput,outfid,sfid,'values','Warning and alarm information as follows: 0=self-check OK; W = At least one Warning active, no Alarms; A = At least one Alarm active.

    c1id=ncdf_vardef(outfid,'first_cbh',[tdid],/long)
    ncdf_attput,outfid,c1id,'long_name','Lowest cloud base height detected.'
    ncdf_attput,outfid,c1id,'units','m'
    ncdf_attput,outfid,c1id,'missing_value','-9999.f'
    ncdf_attput,outfid,c1id,'comment','Field valid if detection_status = 1,2 or 3.'

    vvid=ncdf_vardef(outfid,'vertical_visibility',[tdid],/long)
    ncdf_attput,outfid,vvid,'long_name','Vertical visibility'
    ncdf_attput,outfid,vvid,'units','m'
    ncdf_attput,outfid,vvid,'missing_value','-9999.f'
    ncdf_attput,outfid,vvid,'values','Field only valid if detection_status = 4.'
    ncdf_attput,outfid,vvid,'comment','This is calculated using a proprietary contrast threshold.'

    c2id=ncdf_vardef(outfid,'second_cbh',[tdid],/long)
    ncdf_attput,outfid,c2id,'long_name','Second lowest cloud base height'
    ncdf_attput,outfid,c2id,'units','m'
    ncdf_attput,outfid,c2id,'missing_value','-9999.f'
    ncdf_attput,outfid,c2id,'comment','Field valid if detection_status = 2 or 3'

    hsid=ncdf_vardef(outfid,'alt_highest_signal',[tdid],/long)
    ncdf_attput,outfid,hsid,'long_name','Altitude of highest signal'
    ncdf_attput,outfid,hsid,'units','m'
    ncdf_attput,outfid,hsid,'missing_value','-9999.f'
    ncdf_attput,outfid,hsid,'comment','Field only valid if detection_status = 4.'

    c3id=ncdf_vardef(outfid,'third_cbh',[tdid],/long)
    ncdf_attput,outfid,c3id,'long_name','Third lowest cloud base height'
    ncdf_attput,outfid,c3id,'units','m'
    ncdf_attput,outfid,c3id,'missing_value','-9999.f'
    ncdf_attput,outfid,c3id,'comment','Field valid if detection_status = 3.'

;    lpid=ncdf_vardef(outfid,'laser_pulse_energy',[tdid],/float)
;    ncdf_attput,outfid,lpid,'long_name','Laser pulse energy'
;    ncdf_attput,outfid,lpid,'units','%'
;    ncdf_attput,outfid,lpid,'missing_value','-9999.f'
;    ncdf_attput,outfid,lpid,'values','0 to 999 %'
;    ncdf_attput,outfid,lpid,'comment','Percentage of nominal factory setting.'
;
;    lsid=ncdf_vardef(outfid,'laser_temperature',[tdid],/float)
;    ncdf_attput,outfid,lsid,'long_name','Laser temperature'
;    ncdf_attput,outfid,lsid,'units','C'
;    ncdf_attput,outfid,lsid,'missing_value','-9999.f'
;    ncdf_attput,outfid,lsid,'values','-50 to +99'
;
;    rcid=ncdf_vardef(outfid,'receiver_sensitivity',[tdid],/float)
;    ncdf_attput,outfid,rcid,'long_name','Receiver sensitivity'
;    ncdf_attput,outfid,rcid,'units','%'
;    ncdf_attput,outfid,rcid,'missing_value','-9999.f'
;    ncdf_attput,outfid,rcid,'values','0 to 999 %'
;    ncdf_attput,outfid,rcid,'comment','Percentage of nominal factory setting.'
;
;    wcid=ncdf_vardef(outfid,'window_contamination',[tdid],/float)
;    ncdf_attput,outfid,wcid,'long_name','Window contamination'
;    ncdf_attput,outfid,wcid,'units','mV'
;    ncdf_attput,outfid,wcid,'missing_value','-9999.f'
;    ncdf_attput,outfid,wcid,'values','0 to 2500 mV'
;    ncdf_attput,outfid,wcid,'comment','Millivolts measured at the internal ADC input'
;
;    taid=ncdf_vardef(outfid,'tilt_angle',[tdid],/float)
;    ncdf_attput,outfid,taid,'long_name','Tilt angle'
;    ncdf_attput,outfid,taid,'units','deg'
;    ncdf_attput,outfid,taid,'missing_value','-9999.f'
;    ncdf_attput,outfid,taid,'values','-15 to +90 degrees from vertical'
;    ncdf_attput,outfid,taid,'comment','If Angle Correction is Off detection values are distances not altitudes.'
;
;    blid=ncdf_vardef(outfid,'background_light',[tdid],/float)
;    ncdf_attput,outfid,blid,'long_name','Background light'
;    ncdf_attput,outfid,blid,'units','mV'
;    ncdf_attput,outfid,blid,'missing_value','-9999.f'
;    ncdf_attput,outfid,blid,'values','0 to 2500 mV'
;    ncdf_attput,outfid,blid,'comment','Millivolts measured at the internal ADC input.'
;
;    sbid=ncdf_vardef(outfid,'sum_backscatter',[tdid],/float)
;    ncdf_attput,outfid,sbid,'long_name','SUM of detected an dnormalized backscatter'
;    ncdf_attput,outfid,sbid,'units','1/srad'
;    ncdf_attput,outfid,sbid,'missing_value','-9999.f'
;    ncdf_attput,outfid,sbid,'values','Multipliled by scaling factors times 10^4.'
;
;    bid=ncdf_vardef(outfid,'backscatter',[hdid,tdid],/float)
;    ncdf_attput,outfid,bid,'long_name','Backscatter'
;    ncdf_attput,outfid,bid,'units','1/(srad km)'
;    ncdf_attput,outfid,bid,'missing_value','-9999.f'
;    ncdf_attput,outfid,bid,'comment','Data is range and sensitivity normalized backscatter.'
;
;    mpid=ncdf_vardef(outfid,'measurement_parameters',[sdid,tdid],/char)
;    ncdf_attput,outfid,mpid,'long_name','6 character string describing instrument measurement parameters'
;    ncdf_attput,outfid,mpid,'units','unitless'
;    ncdf_attput,outfid,mpid,'missing_value','-9999'
;    ncdf_attput,outfid,mpid,'values','(L)ong/(S)hort - (F)req - pulse qty r^(7)+1 - gain (H)igh/(L)ow - bandwidth (N)arrow/(W)ide - samping (5) / (1)0 / (2)0 Mhz'
;
;    ssid=ncdf_vardef(outfid,'status_string',[sdid,tdid],/char)
;    ncdf_attput,outfid,ssid,'long_name','Warning and alarm status bits'
;    ncdf_attput,outfid,ssid,'units','unitless'
;    ncdf_attput,outfid,ssid,'missing_value','-9999'
;    ncdf_attput,outfid,ssid,'comment','See Vaisala CT25k manual'   ;ARM format has more information included here
;
;    ltid=ncdf_vardef(outfid,'lat',/float)
;    ncdf_attput,outfid,ltid,'long_name','north latitude'
;    ncdf_attput,outfid,ltid,'units','degrees'
;    ncdf_attput,outfid,ltid,'valid_min','-90.f'
;    ncdf_attput,outfid,ltid,'valid_max','90.f'
;
;    lnid=ncdf_vardef(outfid,'lon',/float)
;    ncdf_attput,outfid,lnid,'long_name','east longitude'
;    ncdf_attput,outfid,lnid,'units','degrees'
;    ncdf_attput,outfid,lnid,'valid_min','-180.f'
;    ncdf_attput,outfid,lnid,'valid_max','180.f'
;
;    alid=ncdf_vardef(outfid,'alt',/float)
;    ncdf_attput,outfid,alid,'long_name','altitude'
;    ncdf_attput,outfid,alid,'units','meters above Mean Sea Level'

    ;write global attributes
    ncdf_attput,outfid,/global,'date_created',systime(0)
    ncdf_attput,outfid,/global,'product_name','vceilcbX1.b1'
    ncdf_attput,outfid,/global,'proc_level','b1'
    ncdf_attput,outfid,/global,'input_data_stream_01','hourly raw Vaisala files (i.e.,'+file[0]+')'
    ncdf_attput,outfid,/global,'description','Vaisala CL31 ceilometer raw measurement data set converted from raw ASCII files (cloud boundaries only)'
    ncdf_attput,outfid,/global,'PI_names','Matthew Shupe'
    ncdf_attput,outfid,/global,'contact','matthew.shupe@noaa.gov'
    ncdf_attput,outfid,/global,'project_name','ACSE: Arctic Clouds in Summer Experiment (SWERUS-C3)'
    ncdf_attput,outfid,/global,'site_id','acs'
    ncdf_attput,outfid,/global,'facility_id','Oden icebreaker in Arctic Ocean'
    ncdf_attput,outfid,/global,'instrument','NOAA ESRL CL31 ceilometer'

    ;Write data
    ncdf_control,outfid,/fill
    ncdf_control,outfid,/endef
    ;Assign data to variables
    ;-------------------------
    ncdf_varput,outfid,btid,base_time
    ncdf_varput,outfid,toid,time_offset
    ncdf_varput,outfid,tid,time
;    ncdf_varput,outfid,rid,height
    ncdf_varput,outfid,dsid,detstat
;    ncdf_varput,outfid,sfid,string(transpose(byte(statflag)))
    ncdf_varput,outfid,c1id,cbase1
    ncdf_varput,outfid,vvid,vertvis
    ncdf_varput,outfid,c2id,cbase2
    ncdf_varput,outfid,hsid,highsig
    ncdf_varput,outfid,c3id,cbase3
;    ncdf_varput,outfid,lpid,lasener
;    ncdf_varput,outfid,lsid,lastemp
;    ncdf_varput,outfid,rcid,recsens
;    ncdf_varput,outfid,wcid,wincont
;    ncdf_varput,outfid,taid,tiltang
;    ncdf_varput,outfid,blid,bckglgt
;    ncdf_varput,outfid,sbid,sumback
;    ncdf_varput,outfid,bid,transpose(backscat)
;    ncdf_varput,outfid,mpid,string(byte(measparam))
;    ncdf_varput,outfid,ssid,string(byte(statstring))
;    ncdf_varput,outfid,ltid,lat
;    ncdf_varput,outfid,lnid,lon
;    ncdf_varput,outfid,alid,alt
    ncdf_close,outfid
    
    
  endif

endfor
end