pro readskyrad,arg,ktime,swd,lwd,KJTIME=kjtime,NO_NAN=no_nan,FILE=file,PLATFORM=platform,TIRK=tirk,SWDDIFF=swddiff,LWDSHAD=lwdshad
;----------------------------------------------------------------------
;Title:		readskyrad.pro
;
;Purpose:	Read the NSA sky radiation (skyrad) netCDF files. 
;
;Input:		arg: 	a string filename argument for the desired file to read.
;
;Output:	gtime: time [decimal hours]
;		swd:  downwelling SW hemispheric irradiance, unshaded [W/m2]
;		lwd:  downwelling LW hemispheric irradiance, unshaded [W/m2]
;Keywords:	kjtime: optional variable for outputed decimal Julian day time
;		nonan:  if set then instead of filling bad data with NaN, the data are left as is.
;		file:  set to an exact filename if only one file is to be read.
;		platform: string variable containing the platform (i.e. nsagndrad60sC1.b1)
;		tirk: Sky/cloud IR temperature [C]
;		swddiff: Downwelling SW diffuse hemispheric irradiance [W/m2]
;		lwdshad: Downwelling LW hemispheric irradiance, shaded [W/m2]
;		
;Note:		RULES:  if swd < -20 W/m2 then NaN
;			if lwd < -20 W/m2 then NaN
;			if swddiff < -20 W/m2 then NaN
;			if lwdshad < -20 W/m2 then NaN
;			if tirk < -100 C then NaN
;
;I/O format:	readskyrad,farg,ktime,swd,lwd,KJTIME=gjtime,TIRK=tirk
;
;Author:	Matthew Shupe
;Date:		5/29/00
;Modified:	5/29/00
;----------------------------------------------------------------------

base_date=julday(1,1,1970,0,0,0)
base=jul_to_dt(base_date)
!dt_base=base

if not(keyword_set(FILE)) then begin
	farg='nsaskyrad*b1.'+arg+'*.cdf'
	file=findfile(farg,count=nfil)
	if nfil eq 0 then begin
		farg='nsaskyrad*a1.'+arg+'*.cdf'
		file=findfile(farg,count=nfil)
	endif
endif else nfil=1

if nfil eq 0 then begin
	ktime=-1
	platform=''
	return
endif

platform=strmid(file[0],0,strpos(file[0],'.2'))
for i=0,nfil-1 do begin
	fid=ncdf_open(file[i])
	ncdf_varget,fid,ncdf_varid(fid,'base_time'),base_time
	ncdf_varget,fid,ncdf_varid(fid,'time_offset'),time_offset
	if ncdf_varid(fid,'down_short_hemisp') ne -1 then begin
		ncdf_varget,fid,ncdf_varid(fid,'down_short_hemisp'),sud
		if ncdf_varid(fid,'down_long_hemisp_shaded1') ne -1 then begin
			ncdf_varget,fid,ncdf_varid(fid,'down_long_hemisp_shaded1'),lud
			ncdf_varget,fid,ncdf_varid(fid,'down_long_hemisp_shaded2'),lsd
		endif else begin
			ncdf_varget,fid,ncdf_varid(fid,'down_long_hemisp_unshaded'),lud
			ncdf_varget,fid,ncdf_varid(fid,'down_long_hemisp_shaded'),lsd
		endelse
		ncdf_varget,fid,ncdf_varid(fid,'down_short_diffuse_hemisp'),ssd
		ncdf_varget,fid,ncdf_varid(fid,'sky_ir_temp'),irt
	endif else begin
		ncdf_varget,fid,ncdf_varid(fid,'psp1_mean'),sud
		ncdf_varget,fid,ncdf_varid(fid,'psps_mean'),ssd
		ncdf_varget,fid,ncdf_varid(fid,'pir1_mean'),lud
		ncdf_varget,fid,ncdf_varid(fid,'pirs_mean'),lsd
		ncdf_varget,fid,ncdf_varid(fid,'irt1_mean'),irt
	endelse
	ncdf_close,fid
	;Create the time axis
	starttime=sec_to_dt(base_time)
	kt=float(starttime.hour+(starttime.minute + (starttime.second+time_offset)/60.)/60.)
	kjt=float(floor(starttime.julian-julday(1,1,starttime.year,0,0,0)+1))+kt/24.
	if n_elements(where(kt gt 24.0)) gt n_elements(kt)/2 then kt=kt-24.
	if kjt[0] ge 365 and fix(strmid(arg,4,2)) eq 1 then kjt=kjt-fix(kjt[0])+1.0
	if i eq 0 then begin
		swd=sud
		swddiff=ssd
		lwd=lud
		lwdshad=lsd
		tirk=irt
		ktime=kt
		kjtime=kjt
	endif else begin
		swd=[swd,sud]
		swddiff=[swddiff,ssd]
		lwd=[lwd,lud]
		lwdshad=[lwdshad,lsd]
		tirk=[tirk,irt]
		ktime=[ktime,kt]
		kjtime=[kjtime,kjt]
	endelse
endfor

;Set all "Missing" or obviously bad data to NaN  (swd<0,lwd<0,swddiff<0,lwdshad<0,tirk<-100C)
if not(keyword_set(NO_NAN)) then begin
	iwh=where(swd lt -20)
	if iwh[0] ne -1 then swd[iwh]=!values.f_nan
	iwh=where(swd lt 0 and swd ge -20)
	if iwh[0] ne -1 then swd[iwh]=0
	iwh=where(swddiff lt -20)
	if iwh[0] ne -1 then swddiff[iwh]=!values.f_nan
	iwh=where(swddiff lt 0 and swddiff ge -20)
	if iwh[0] ne -1 then swddiff[iwh]=0
	iwh=where(lwd lt -20)
	if iwh[0] ne -1 then lwd[iwh]=!values.f_nan
	iwh=where(lwd lt 0 and lwd ge -20)
	if iwh[0] ne -1 then lwd[iwh]=0
	iwh=where(lwdshad lt -20)
	if iwh[0] ne -1 then lwdshad[iwh]=!values.f_nan
	iwh=where(lwdshad lt 0 and lwdshad ge -20)
	if iwh[0] ne -1 then lwdshad[iwh]=0
	iwh=where(tirk lt -100)
	if iwh[0] ne -1 then tirk[iwh]=!values.f_nan
endif
end