;Title:  CT25k_to_netcdf.pro
;Purpose:   Read the CT25k ceilometer data and make a netCDF file
;
;Inputs:	sarg:  starting date/time string, YYYYMMDD
;		earg:  ending date/time string, YYYYMMDD
;Outputs:	NetCDF file containing the following fields:
;		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
;		backscat: backscatter in units of (srad*km)^-1
;		height: height axis [m].
;Keywords:
;		directory: optionally specify a directory that contains the data.
;
;I/O:    	CT25k_to_netcdf,....
;
;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 ct25k_to_netcdf,sarg,earg,directory=directory

mkDAM=1   ;if set then make DAMOCLES-type files

;Got to the correct directory if specified
if keyword_set(directory) then cd,directory,current=orig_dir

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

;Determine the daily argument list
if sarg ne earg then begin
	sjd=julday(strmid(sarg,4,2),strmid(sarg,6,2),strmid(sarg,0,4))
	ejd=julday(strmid(earg,4,2),strmid(earg,6,2),strmid(earg,0,4))
	numa=ejd-sjd+1
	darg=strarr(numa)
	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)
		darg[i]=yr+mo+da
	endfor
endif else begin
	darg=strmid(sarg,3,5)
	numa=n_elements(darg)
endelse

;Loop over days
for d=0,numa-1 do begin

	;Get mean lat and lon for the day
	ascos_readgps,'200'+darg[d]+'00','200'+darg[d]+'24',gtime,glat,glon,directory='/cruises/ASCOS_2008/ODEN/gps'
	wh=where(glat gt 70 and glat lt 90,nm)
	if nm gt 0 then lat=mean(glat[wh]) else lat=-999.
	wh=where(glon gt -20 and glon lt 50,nm)
	if nm gt 0 then lon=mean(glon[wh]) else lon=-999.
	alt=10.0

	;Initialize variables
	ind=0L
	num=numa*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
	statflag=strarr(num)
	statstring=strarr(num)
	year=intarr(num)
	mon=intarr(num)
	day=intarr(num)
	time=fltarr(num)
	backscat=fltarr(num,256)
	height=indgen(256)*30 + 15 ; 30m range gates

	;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
	arg='*'+darg[d]+'*.DAT'
	file=findfile(arg,count=numf)
	;Loop over all files in the day
	for f=0,numf-1 do begin
		openr,lun,/get_lun,file[f]
		readf,lun,jnk
		readf,lun,jnk
		while not eof(lun) do begin
			readf,lun,str,format="(a10)"
			if strmid(str,0,2) eq '-F' then goto,flag
			readf,lun,jnk
			readf,lun,dt1,dt2,dt3,dt4,dt5,dt6,dt7,dt8,dt9,format="(a2,a1,a5,a1,a5,a1,a5,a1,a8)"
			readf,lun,jnk		;parameter line:  100 N 99 +21 .... etc.
			;Store data
			detstat[ind]=fix(strmid(dt1,0,1))
			statflag[ind]=strmid(dt1,1,1)
			statstring[ind]=dt9
			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,2,2))+float(strmid(str,5,2))/60.+float(strmid(str,8,2))/3600.
			for i=0,15 do begin
				readf,lun,str
				str=strmid(str,3,64)  ;Trim the string down to only the backscatter data
				for j=0,15 do begin
					tmp=execute("backscat[ind,i*16+j]='"+strmid(str,4*j,4)+"'X")   ;execute the conversion between hex to int
				endfor
			endfor
			readf,lun,jnk   ;blank line between records
			ind=ind+1
		endwhile
		flag:
		free_lun,lun
	endfor
	;Deal with data formating
	if ind gt 0 then begin
		;trim the data down to actual size
		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]
		backscat=backscat[0:ind-1,*]
		statflag=statflag[0:ind-1]
		statstring=statstring[0:ind-1]
		;Deal with the 2's complement issue for backscatter
		nbits=16
		wh=where(backscat ge 2.0^(nbits-1),nm)
		if nm gt 0 then backscat[wh]=backscat[wh]-2.0^(nbits)
		;Set units
		backscat=backscat*1e-4  ;Convert units: (sr km)^-1
	endif else begin
		;If there is no data thenb fill dummy data
		time=findgen(1440)/1440*24.
		cbase1=time*0-999
		cbase2=cbase1
		cbase3=cbase1
		vertvis=cbase1
		highsig=cbase1
		detstat=time*0-1
		backscat=backscat[0:1439,*]*0-999.
	endelse
	;compute base_time and time_offset
	sec=dt_to_sec( str_to_dt( strmid(darg[d],1,2)+'/'+strmid(darg[d],3,2)+'/200'+strmid(darg[d],0,1) ) )
	base_time=sec+fix(time[0]*3600)
	time_offset=long((time-time[0])*3600)
	;Write an output file for the day
	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='200'+strmid(darg[d],0,1)+'-'+strmid(darg[d],1,2)+'-'+strmid(darg[d],3,2)+' '+hr+':'+mi+':'+se+' 0:00'

	if mkDAM then outfilename='ascos_vceil25k.200'+darg[d]+'.'+hr+mi+se+'.nc' else outfilename='ascvceil25k.200'+darg[d]+'.'+hr+mi+se+'.cdf'
	nbeam=n_elements(time)
	nheights=n_elements(height)
	print,'WRITING: ',outfilename
	outfid=ncdf_create(outfilename,/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 sinde 1970-1-1 0:00:00 0:00'
	if mkDAM then begin
		toid=ncdf_vardef(outfid,'time',[tdid],/double)
		ncdf_attput,outfid,toid,'long_name','time'
		ncdf_attput,outfid,toid,'standard_name','time'
		ncdf_attput,outfid,toid,'units','seconds since '+basetimestring
		ncdf_attput,outfid,toid,'axis','T'
	endif else begin
		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
	endelse
	if not mkDAM then begin
		tid=ncdf_vardef(outfid,'time',[tdid],/float)
		ncdf_attput,outfid,tid,'long_name','Time'
		ncdf_attput,outfid,tid,'units','Decimal hours'
	endif
	dsid=ncdf_vardef(outfid,'detection_status',[tdid],/long)
	ncdf_attput,outfid,dsid,'long_name','Detection status'
	ncdf_attput,outfid,dsid,'values','0=No significant backscatter; 1-3=Number of cloud bases detected; 4=Full obscuration determined but no cloud; 5=Some obscuration determined but detected to be transparent.'
	ncdf_attput,outfid,dsid,'units','none'
	c1id=ncdf_vardef(outfid,'first_cbh',[tdid],/long)
	ncdf_attput,outfid,c1id,'long_name','Height from ground to lowest cloud base.'
	ncdf_attput,outfid,c1id,'units','meters'
	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','meters'
	ncdf_attput,outfid,vvid,'comment','Field only valid if detection_status = 4.'
	c2id=ncdf_vardef(outfid,'second_cbh',[tdid],/long)
	ncdf_attput,outfid,c2id,'long_name','Height from ground to second cloud base'
	ncdf_attput,outfid,c2id,'units','meters'
	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','meters'
	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','Height from ground to third cloud base"
	ncdf_attput,outfid,c3id,'units','meters'
	ncdf_attput,outfid,c3id,'comment','Field valid if detection_status = 3.'
	sfid=ncdf_vardef(outfid,'status_flag',[tdid],/char)
	ncdf_attput,outfid,sfid,'long_name','Ceilometer status indicator'
	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.
	ncdf_attput,outfid,sfid,'units','none'
	ssid=ncdf_vardef(outfid,'status_string',[tdid,sdid],/char)
	ncdf_attput,outfid,ssid,'long_name','Warning and alarm status bits'
	ncdf_attput,outfid,ssid,'units','unitless'
	ncdf_attput,outfid,ssid,'comment','See Vaisala CT25k manual'
	if mkDAM then begin
		rid=ncdf_vardef(outfid,'height',[hdid],/long)
		ncdf_attput,outfid,rid,'long_name','height'
		ncdf_attput,outfid,rid,'units','meters'
		ncdf_attput,outfid,rid,'comment','Height to center of range bin'
		ncdf_attput,outfid,rid,'axis','Z'
	endif else begin
		rid=ncdf_vardef(outfid,'range',[hdid],/long)
		ncdf_attput,outfid,rid,'long_name','Distance to center of range bin'
		ncdf_attput,outfid,rid,'units','meters'
		ncdf_attput,outfid,rid,'comment','none'
	endelse
	bid=ncdf_vardef(outfid,'backscatter',[tdid,hdid],/float)
	ncdf_attput,outfid,bid,'long_name','Backscatter'
	ncdf_attput,outfid,bid,'units','1/(srad km)'
	ncdf_attput,outfid,bid,'comment','Data is range and sensitivity normalized backscatter.'
	if mkDAM then ltid=ncdf_vardef(outfid,'latitude',/float) else ltid=ncdf_vardef(outfid,'lat',/float)
	ncdf_attput,outfid,ltid,'long_name','latitude'
	ncdf_attput,outfid,ltid,'standard_name','latitude'
	ncdf_attput,outfid,ltid,'units','degree_north'
	ncdf_attput,outfid,ltid,'valid_min','-90.f'
	ncdf_attput,outfid,ltid,'valid_max','90.f'
	if mkDAM then lnid=ncdf_vardef(outfid,'longitude',/float) else lnid=ncdf_vardef(outfid,'lon',/float)
	ncdf_attput,outfid,lnid,'long_name','longitude'
	ncdf_attput,outfid,lnid,'standard_name','longitude'
	ncdf_attput,outfid,lnid,'units','degree_east'
	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,'title','Ceilometer backscatter and cloud base height measurements.'
	ncdf_attput,outfid,/global,'abstract','This data set includes measurements from a ground-based, vertically-pointing, ceilometer, which includes backscatter profiles, cloud base heights for up to three layers, and the vertical visibility.'
	ncdf_attput,outfid,/global,'topiccategory','ClimatologyMeteorologyAtmosphere'
	ncdf_attput,outfid,/global,'keywords','Atmospheric Ceiling Cloud Observation Measurement Polar Reflectivity Tropospheric'
	ncdf_attput,outfid,/global,'gcmd_keywords',' '
	ncdf_attput,outfid,/global,'activity_type','Cruise, Ice station'
	ncdf_attput,outfid,/global,'Conventions','CF-1.0'
	ncdf_attput,outfid,/global,'product_name','vceil25kC1.b1'
	ncdf_attput,outfid,/global,'history','2009-08-20 creation'
	ncdf_attput,outfid,/global,'area','Arctic Ocean'
	ncdf_attput,outfid,/global,'southernmost_latitude','56.48'
	ncdf_attput,outfid,/global,'northernmost_latitude','87.49'
	ncdf_attput,outfid,/global,'westernmost_longitude','-11.08'
	ncdf_attput,outfid,/global,'easternmost_longitude','20.78'
	ncdf_attput,outfid,/global,'start_date','2008-08-13 20:11:00 UTC'
	ncdf_attput,outfid,/global,'stop_date','2008-09-16 06:00:00 UTC'
	ncdf_attput,outfid,/global,'institution','Stockholm University'
	ncdf_attput,outfid,/global,'PI_name','Matthew Shupe'
	ncdf_attput,outfid,/global,'contact','matthew.shupe@noaa.gov'
	ncdf_attput,outfid,/global,'distribution_statement','Free'
	ncdf_attput,outfid,/global,'project_name','Arctic Summer Cloud Ocean Study'

	ncdf_attput,outfid,/global,'site_id','asc'
	ncdf_attput,outfid,/global,'facility_id','Icebreaker ODEN during ASCOS cruise, Arctic Ocean'
	ncdf_attput,outfid,/global,'instrument','NOAA/ESRL/PSD3 CT25k ceilometer'
	ncdf_attput,outfid,/global,'PI_affiliation','University of Colorado/CIRES and NOAA ESRL'
	ncdf_attput,outfid,/global,'support','This material is based upon work supported by the National Science Foundation under Grant No. ARC0732925.'

	;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
	if not mkDAM then ncdf_varput,outfid,tid,time
	ncdf_varput,outfid,dsid,detstat
	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,sfid,transpose(byte(statflag))
	ncdf_varput,outfid,ssid,transpose(byte(statstring))
	ncdf_varput,outfid,rid,height
	ncdf_varput,outfid,bid,backscat
	ncdf_varput,outfid,ltid,lat
	ncdf_varput,outfid,lnid,lon
	ncdf_varput,outfid,alid,alt
	ncdf_close,outfid

endfor ;days
if keyword_set(directory) then cd,orig_dir
end