pro plot_physprop_input, date1
;date1='20001128'

common five_min_PhysProp, hrfrac,jday,height,reflectivity,reflectp,reflectc,cbase,ctop, $
	velocity, spectrum, temp, pressure,fluxclear,solarratio,lwp,vceil,vap, $
	IWC_best_estimate, re_ice_best_estimate, ice_micro_error_estimate, iwc_source, $
	ice_effective_radius_source, ice_tau_solar, ice_omega_solar, ice_g_solar, ice_tau_IR,$
    ice_omega_IR, ice_g_IR, ice_RadiationParameterization_source_solar, $
	ice_RadiationParameterization_source_ir, lwc_best_estimate, re_liquid_best_estimate, $
    liquid_micro_error_estimate, lwc_source, liquid_effective_radius_source, liquid_tau_solar, $
	liquid_omega_solar, liquid_g_solar, liquid_tau_IR, liquid_omega_IR, liquid_g_IR, $
	liquid_RadiationParameterization_source_solar, liquid_RadiationParameterization_source_ir, $
	ice_fraction_CCM3, solar_flux_down, solar_flux_up, IR_flux_up, IR_flux_down, $
	downwelling_solar_diffuse, downwelling_solar_direct, TOA_solar_up, TOA_solar_down, TOA_IR_up, $
	solar_flux_down_clear, solar_flux_up_clear, IR_flux_up_clear, IR_flux_down_clear, $
	downwelling_solar_diffuse_clear, downwelling_solar_direct_clear, $
	TOA_solar_up_clear, TOA_solar_down_clear, TOA_IR_up_clear

common five_min_PhysProp2, down_long,down_short,up_long,up_short

common sgpgoes8minnisX1_data, hrfrac_M, jday_M, num_times_M, num_longs_M, num_lats_M, num_levels, num_views, $
	latitude_M, longitude_M, level, view, Cloud_Amount, Visible_Optical_Depth, IR_Optical_Depth, Emissivity, $
	Cloud_Center_Height, Cloud_Top_Height, Cloud_Temperature, Cloud_Thickness, Reflectance, Albedo, $
	Cloud_Center_Temperature, Cloud_Top_Temperature, Visible_Optical_Depth_SD, Cloud_Center_Temperature_SD, $
	Broadband_LW_Flux, Broadband_SW_Albedo, Narrowband_VIS_Albedo, Clear_Temperature, Clear_Temperature_SD, $
	Narrowband_VIS_Albedo_SD, Clear_VIS_Reflectance, Average_Total_Temperature, Solar_Zenith_Angle, $
	Viewing_Zenith_Angle, Relative_Azimuth_Angle

;[sfc_temp_vector],[sfc_mxrat_vector], [hrfrac], [vap_vector]


 mlr_coeff_up=[-1379.1984,6.1345426,-0.78215353,1.2752958, -2.9191544]


;IDL> print, mlr_coeff_up
;       6.1345426
;     -0.78215353
;       1.2752958
;      -2.9191544
;IDL> print, mlr_constu
;      -1379.1984



mlr_coeff_down=[-284.89975,1.7988579,3.9060878,0.13083693,26.946193]
;IDL> print, mlr_constd
;      -284.89975
;IDL> print, mlr_coeff_down
;       1.7988579
;       3.9060878
;      0.13083693
;       26.946193



;date1='20000328'
date2='000000'
site='sgp'
;fname='d:\sgp.average.5min\sgp.average.5min.20000313.000000.temp.cdf '
;fname='c:\mace\20010411\'+site+'.average.5min.'+date1+'.'+date2+'.temp.cdf '
;fname_g_prefix='c:\mace\20010411\'+site+'.average.5min.'+date1+'.'+date2+'.'

;fname='e:\'+site+'.average.5min\FullRet\'+site+'.average.5min.'+date1+'.'+date2+'.temp.cdf '
;fname_g_prefix='c:\mace\'+site+'.average.5min\FullRet\'+site+'.average.5min.'+date1+'.'+date2+'.'

path='/data/mace3/arm_data/sgp/sgp.average.5min/'
minnis_path='/home/mace/sgp/sgpgoes8minnis-acfX1.c1/'

fname=path+site+'.average.5min.'+date1+'.'+date2+'.temp.cdf '
print, 'fname is ',strtrim(fname,2)
fname=strtrim(fname,2)
spawn, 'gunzip '+fname
fname_clear=path+site+'.average.5min.'+date1+'.'+date2+'.temp.cdf'
fname_g_prefix=path+site+'.average.5min.'+date1+'.'+date2+'.'


fname_g_suffix='.png'

title_string='ARM Integrated Cloud Product from '+site+' on '+date1

close, /all
OPENR, 1, strtrim(fname,2), ERROR = err
print, 'err is ',err
if err eq 0 then begin

close, 1

spawn, 'gunzip -f '+strtrim(fname,2)+'.gz'

get_5min_avg_PhysProp, strtrim(fname,2), fname_clear, sfc_mxrat, sfc_temp, solar_zenith_anglexx

;spawn, 'gzip -f '+strtrim(fname,2)

flag_it, data_flag

print, 'data flag is ',data_flag
if data_flag ne 0 then begin

fname_minnis=minnis_path+'sgpgoes8minnis-acfX1.c1.'+date1+'.000000.cdf'
;fname_minnis=' '

spawn, 'gunzip -f '+fname_minnis+'.gz'
if fname_minnis ne ' ' then get_sgpgoes8minnisx1, fname_minnis
;spawn, 'gzip -f '+fname_minnis

print, jday[0:20],fname_minnis
print, hrfrac_M, num_times_M

; need the ability to plot just part of the arrays:

start_time=0. & end_time=23.9
start_height=0.0 & end_height=14000.0
loadct, 15
tvlct, red, green, blue, /get

start_time_index=0
while hrfrac(start_time_index) lt start_time do start_time_index=start_time_index+1 & end_time_index=start_time_index
while hrfrac(end_time_index) lt end_time do end_time_index=end_time_index+1
start_height_index=0
while height(start_height_index) lt start_height do start_height_index=start_height_index+1 & end_height_index=start_height_index
while height(end_height_index) lt end_height do end_height_index=end_height_index+1

image=fltarr(end_time_index-start_time_index+1,end_height_index-start_height_index+1)
plot_xaxis=hrfrac[start_time_index:end_time_index]
plot_yaxis=height[start_height_index:end_height_index]

one_left_pos=[0.05,0.825,0.45,0.925] & one_left_pos2=[0.45+0.045,0.825,0.45+0.05,0.925]; upper left
one_right_pos=[0.535,0.825,0.935,0.925]	& one_right_pos2=[0.935+0.045,0.825,0.935+0.05,0.925]; upper right
two_left_pos=[0.05,0.675,0.45,0.775] & two_left_pos2=[0.05+0.445,0.675,0.45+0.05,0.775]	; 2nd left
two_right_pos=[0.535,0.675,0.935,0.775] & two_right_pos2=[0.935+0.045,0.675,0.935+0.05,0.775]	; 2nd right
three_left_pos=[0.05,0.525,0.45,0.625] & three_left_pos2=[0.45+0.045,0.525,0.45+0.05,0.625]; 3rd left
three_right_pos=[0.535,0.525,0.935,0.625] & three_right_pos2=[0.935+0.045,0.525,0.935+0.05,0.625]	; 3rd right
four_left_pos=[0.05,0.375,0.45,0.475] & four_left_pos2=[0.45+0.045,0.375,0.45+0.05,0.475]	; 4th left
four_right_pos=[0.535,0.375,0.935,0.475] & four_right_pos2=[0.935+0.045,0.375,0.935+0.05,0.475]	; 4th right
five_left_pos=[0.05,0.225,0.45,0.325] & five_left_pos2=[0.45+0.045,0.225,0.45+0.05,0.325]	; 5th left
five_right_pos=[0.535,0.225,0.935,0.325] & five_right_pos2=[0.935+0.045,0.225,0.935+0.05,0.325]	;5th right
six_left_pos=[0.05,0.075,0.45,0.175] & six_left_pos2=[0.05+0.445,0.075,0.45+0.05,0.175]	; 6th left
six_right_pos=[0.535,0.075,0.935,0.175] & six_right_pos2=[0.935+0.045,0.075,0.935+0.05,0.175]	; 6th right

;print, n_elements(solar_flux_down), 'is the fucking n_elements(solar_flux_down)'
if n_elements(solar_flux_down) gt 1 then begin
window_1='yes'
window_2='yes'
window_3='yes'
window_4='no'
window_5='yes'
window_6='yes'
endif else begin
window_1='no'
window_2='no'
window_3='no'
window_4='no'
window_5='no'
window_6='no'
print, window_1, window_2, window_3, window_4, window_5, window_6, 'is the fucking window variables'
endelse


;if window_1 ne 'no' then begin


if window_1 ne 'no' then WINDOW, 1, XSIZE=1000, YSIZE=1000, /pixmap,TITLE='ARM Integrated Column Product Plots (ice)'



!p.multi = [0,6,6]
;height=height/1000.
!p.charsize=2.0

if window_1 ne 'no' then begin

;solar flux down

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_down[kk,jj] gt 0. and solar_flux_down[kk,jj] lt 1500. then image[j,k]=solar_flux_down[kk,jj] $
			else if solar_flux_down[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

; radar reflectivity

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if reflectivity[kk,jj] ne -9999. then begin
				image[j,k]=reflectivity[kk,jj]/100.
			endif else begin
				if reflectivity[kk,jj] eq -9999. then image[j,k]=-9999.
			endelse
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.825 & xr=0.45 & yr=0.925	; upper left


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Radar Reflectivity', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'DbZ',10., -50.

xyouts, 0.5,0.97, title_string, /normal, alignment=0.5, color=0

; Doppler Velocity

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if velocity[kk,jj] ne -9999. then image[j,k]=velocity[kk,jj]/1000. $
			else if velocity[kk,jj] eq -9999. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.825 & xr=0.935 & yr=0.925	; upper right


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Doppler Velocity', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'm/s',1.5, -0.5

; Temperature

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if temp[kk,jj] gt -9999. then image[j,k]=temp[kk,jj]-273. $
			else if temp[kk,jj] le -9999. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.675 & xr=0.45 & yr=0.775 ; middle left


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Temperature', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'Celcius',max(image), -80.

; Pressure

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if pressure[kk,jj] ne -9999. then image[j,k]=pressure[kk,jj]/100. $
			else if pressure[kk,jj] eq -9999. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.675 & xr=0.935 & yr=0.775 ; middle right


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Pressure', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'mb',max(image), min(image)

; IWC

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IWC_best_estimate[kk,jj] ne -8000 then image[j,k]=(float(IWC_best_estimate[kk,jj])/1000.) $
			else if IWC_best_estimate[kk,jj] lt -8045 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.525 & xr=0.45 & yr=0.625 ; middle left 2


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Ice Water', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'log(g/m3)',0.1, -5.

; ice effective radius

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if re_ice_best_estimate[kk,jj] ne -8000 then image[j,k]=(1.e6)*(10.^(float(re_ice_best_estimate[kk,jj])/1000.)) $
			else if re_ice_best_estimate[kk,jj] lt -8045 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.525 & xr=0.935 & yr=0.625 ; middle right 2


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'ice effective radius', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'microns',150., 0.

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if iwc_source[kk,jj] gt 0 then image[j,k]=float(iwc_source[kk,jj]) $
			else if iwc_source[kk,jj] le 0 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.375 & xr=0.45 & yr=0.475 ; middle left 2


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'IWC source', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'IWC Source',max(image), 0.


;ice_effective_radius_source
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if ice_effective_radius_source[kk,jj] gt 0 then image[j,k]=float(ice_effective_radius_source[kk,jj]) $
			else if ice_effective_radius_source[kk,jj] le 0 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.375 & xr=0.935 & yr=0.475 ; middle right 2


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'ice re source', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 're Source',max(image), 0.

endif	; window 1

; visible extinction

jj=start_time_index
vector2=fltarr(n_elements(plot_xaxis))
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if ice_effective_radius_source[kk,jj] gt 0 then begin
				image[j,k]=((1.e3)/90.)*(10.^(float(ice_tau_solar[kk,jj])/1000.))
				vector2[j]=vector2[j]+(image[j,k]*.090)
			endif else if ice_effective_radius_source[kk,jj] le 0 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.225 & xr=0.45 & yr=0.325 ; middle left 2


if window_1 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'ice visible extinction', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], '1/km',max(image), -0.25



; IR extinction

jj=start_time_index
vector1=fltarr(n_elements(plot_xaxis))
image2=image
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if ice_effective_radius_source[kk,jj] gt 0 then begin
				image[j,k]=((1.e3)/90.)*(10.^(float(ice_tau_IR[kk,jj])/1000.))
				ssa=10.^((float(ice_omega_IR[kk,jj]))/1000.)
				;print, ssa,image[j,k],10.^((float(ice_g_IR[kk,jj]))/1000.)
				image2[j,k]=image[j,k]*(1.-ssa)
			endif else begin
				if ice_effective_radius_source[kk,jj] le 0 then begin
						image[j,k]=-9999. & image2[j,k]=-9999.
				endif
			endelse
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.225 & xr=0.935 & yr=0.325 ; middle right 2


if window_1 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'ice IR extinction', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], '1/km',max(image), -0.25

xl=0.05 & yl=0.075 & xr=0.45 & yr=0.175 ; middle left 2


if window_1 ne 'no' then plot, plot_xaxis, vector2, psym=-2, pos=[xl,yl,xr,yr], ytitle='Tau', /ylog, $
	title='ice Visible Optical Depth',background=!d.n_colors-1, color=0, yrange=[0.01, max(vector2)]

ice_tau_vector_solar=vector2

;plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
;	 'DownWelling Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
;	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',1400., 0.


jj=start_time_index
vector1=fltarr(n_elements(plot_xaxis))
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if ice_effective_radius_source[kk,jj] gt 0 then begin
				vector1[j]=vector1[j]+(0.090*image2[j,k])
			endif else begin
				if ice_effective_radius_source[kk,jj] le 0 then begin
						image[j,k]=-9999. & image2[j,k]=-9999.
				endif
			endelse
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.075 & xr=0.935 & yr=0.175 ; middle right 2

if window_1 ne 'no' then plot, plot_xaxis, (1.-exp(-vector1)), psym=-2, pos=[xl,yl,xr,yr], ytitle='ice IR Emissivity', $
	title='ice IR Emissivity and IR Optical Depth',background=!d.n_colors-1, color=0, yrange=[1.e-2, 10.], /ylog

if window_1 ne 'no' then oplot, plot_xaxis, vector1, psym=-1, color=0

ice_tau_vector_ir=vector1

;plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
;	 'DownWelling Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
;	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',1400., 0.



;solar_flux_down, solar_flux_up, IR_flux_up, IR_flux_down, $
;	downwelling_solar_diffuse, downwelling_solar_direct, TOA_solar_up, TOA_solar_down, TOA_IR_up


if window_1 ne 'no' then write_png, fname_g_prefix+'plot1'+fname_g_suffix, tvrd(), red, green, blue
if window_1 ne 'no' then spawn, 'convert '+fname_g_prefix+'plot1'+fname_g_suffix+' '+fname_g_prefix+'plot1.gif'

;endif ; if window_1 ne 'no' then begin


;if window_2 ne 'no' then begin

if window_2 ne 'no' then WINDOW, 2, XSIZE=1000, YSIZE=1000, TITLE='ARM Integrated Column Product Plots (liquid)', /pixmap

!p.multi = [0,6,6]
;height=height/1000.
!p.charsize=2.0

if window_2 ne 'no' then begin

; radar reflectivity	; 1st left plot

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if reflectivity[kk,jj] ne -9999. then image[j,k]=reflectivity[kk,jj]/100. $
			else if reflectivity[kk,jj] eq -9999. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Radar Reflectivity', ' ', 'Height (km)', one_left_pos, $
	one_left_pos2, 'DbZ',10., -40.

xyouts, 0.5,0.97, title_string, /normal, alignment=0.5, color=0

; Doppler Velocity	; first right plot

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if lwc_source[kk,jj] gt 0 or iwc_source[kk,jj] gt 0 then begin
					image[j,k]=10.^((float(ice_fraction_CCM3[kk,jj]))/1000.)
					;print, image[j,k]
			endif else begin
				 image[j,k]=-9999.
			endelse
			kk=kk+1
	endfor
jj=jj+1
endfor


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'CCM3 liquid-Ice Fraction', ' ', 'Height (km)', one_right_pos, $
	one_right_pos2, 'fraction',1.5, 0.

; LWC	; 2nd right plot

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if LWC_best_estimate[kk,jj] ne -8000 then image[j,k]=(float(LWC_best_estimate[kk,jj])/1000.) $
			else if LWC_best_estimate[kk,jj] lt -8045 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Liquid Water', ' ', 'Height (km)', two_left_pos, $
	two_left_pos2, 'log(g/m3)',2., -2.

; liquid effective radius

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if re_liquid_best_estimate[kk,jj] ne -8000 then image[j,k]=(1.e6)*(10.^(float(re_liquid_best_estimate[kk,jj])/1000.)) $
			else if re_liquid_best_estimate[kk,jj] lt -8045 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'liquid effective radius', ' ', 'Height (km)', two_right_pos, $
	two_right_pos2, 'microns',50., 0.

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if lwc_source[kk,jj] gt 0 then image[j,k]=float(lwc_source[kk,jj]) $
			else if lwc_source[kk,jj] le 0 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'LWC source', ' ', 'Height (km)', three_left_pos, $
	three_left_pos2, 'LWC Source',max(image)+2, 0.


;ice_effective_radius_source
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if liquid_effective_radius_source[kk,jj] gt 0 then image[j,k]=float(liquid_effective_radius_source[kk,jj]) $
			else if liquid_effective_radius_source[kk,jj] le 0 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.375 & xr=0.935 & yr=0.475 ; middle right 2


plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'liquid re source', ' ', 'Height (km)', three_right_pos, $
	three_right_pos2, 're Source',max(image)+2., 0.

endif
; visible extinction

jj=start_time_index
vector2=fltarr(n_elements(plot_xaxis))
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if liquid_effective_radius_source[kk,jj] gt 0 then begin
				image[j,k]=((1.e3)/90.)*(10.^(float(liquid_tau_solar[kk,jj])/1000.))
				vector2[j]=vector2[j]+(image[j,k]*.090)
			endif else if liquid_effective_radius_source[kk,jj] le 0 then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.225 & xr=0.45 & yr=0.325 ; middle left 2


if window_2 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'liquid visible extinction', ' ', 'Height (km)', four_left_pos, $
	four_left_pos2, '1/km',100., 0.

;endif

; IR extinction

jj=start_time_index
vector1=fltarr(n_elements(plot_xaxis))
image2=image
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if liquid_effective_radius_source[kk,jj] gt 0 then begin
				image[j,k]=((1.e3)/90.)*(10.^(float(liquid_tau_IR[kk,jj])/1000.))
				ssa=10.^((float(liquid_omega_IR[kk,jj]))/1000.)
				;print, ssa,image[j,k],10.^((float(ice_g_IR[kk,jj]))/1000.)
				image2[j,k]=image[j,k]*(1.-ssa)
			endif else begin
				if liquid_effective_radius_source[kk,jj] le 0 then begin
						image[j,k]=-9999. & image2[j,k]=-9999.
				endif
			endelse
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.225 & xr=0.935 & yr=0.325 ; middle right 2


if window_2 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'liquid IR extinction', ' ', 'Height (km)', four_right_pos, $
	four_right_pos2, '1/km',100., 0.

xl=0.05 & yl=0.075 & xr=0.45 & yr=0.175 ; middle left 2


if window_2 ne 'no' then plot, plot_xaxis, vector2, psym=-3, pos=five_left_pos, ytitle='Tau', /ylog, $
	title='liquid and Total (dashed) Visible Optical Depth',background=!d.n_colors-1, $
	color=0, yrange=[1.0, 100.]

total_visible_optical_depth=vector2

if n_elements(ice_tau_vector_solar) eq n_elements(plot_xaxis) then begin

if window_2 ne 'no' then oplot, plot_xaxis, ice_tau_vector_solar+vector2, psym=-3, linestyle=2, color=0

total_visible_optical_depth=vector2+ice_tau_vector_solar

endif

;plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
;	 'DownWelling Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
;	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',1400., 0.


jj=start_time_index
vector1=fltarr(n_elements(plot_xaxis))
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if liquid_effective_radius_source[kk,jj] gt 0 then begin
				vector1[j]=vector1[j]+(0.090*image2[j,k])
			endif else begin
				if liquid_effective_radius_source[kk,jj] le 0 then begin
						image[j,k]=-9999. & image2[j,k]=-9999.
				endif
			endelse
			kk=kk+1
	endfor
jj=jj+1
endfor


if window_2 ne 'no' then plot, plot_xaxis, (1.-exp(-vector1)), psym=-3, pos=five_right_pos, ytitle='liquid IR Emissivity/ O.D.', $
	title='liquid IR Emissivity, IR (+) and Total (dashed) Optical Depth',$
	background=!d.n_colors-1, color=0, yrange=[1.e-1, 100.], /ylog

if window_2 ne 'no' then oplot, plot_xaxis, vector1, psym=-1, color=0

total_ir_optical_depth=vector1

if n_elements(ice_tau_vector_ir) eq n_elements(plot_xaxis) then begin

if window_2 ne 'no' then oplot, plot_xaxis, ice_tau_vector_ir+vector1, psym=-3, linestyle=2, color=0

total_ir_optical_depth=ice_tau_vector_ir+vector1

endif



; liquid and vapor

jj=start_time_index
vector3=fltarr(n_elements(plot_xaxis))
vector4=vector3
for j=0,n_elements(plot_xaxis)-1 do begin

			if vap[jj] gt 0 then begin
				vector3[j]=vap[jj]
			endif
			if lwp[jj] gt 0 then begin
				vector4[j]=lwp[jj]
			endif

jj=jj+1
endfor

if window_2 ne 'no' then plot, plot_xaxis, vector3, psym=-3, pos=six_left_pos, ytitle='Water Paths', $
	title='MWR Vapor (cm) and MWR liquid (mm - dashed)',$
	background=!d.n_colors-1, color=0, yrange=[1.e-2, 100.], /ylog

if window_2 ne 'no' then oplot, plot_xaxis, vector4, psym=-3, linestyle=3, color=0

; cloud base

jj=start_time_index
vector3=fltarr(n_elements(plot_xaxis))
for j=0,n_elements(plot_xaxis)-1 do begin

			if cbase[jj] gt 0 then begin
				vector3[j]=cbase[jj]/1000.
			endif

jj=jj+1
endfor

if window_2 ne 'no' then plot, plot_xaxis, vector3, psym=-3, pos=six_right_pos, ytitle='Water Paths', $
	title='Cloud Base (km)',$
	background=!d.n_colors-1, color=0, yrange=[-1., 15.]

;plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
;	 'DownWelling Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
;	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',1400., 0.



;solar_flux_down, solar_flux_up, IR_flux_up, IR_flux_down, $
;	downwelling_solar_diffuse, downwelling_solar_direct, TOA_solar_up, TOA_solar_down, TOA_IR_up

if window_2 ne 'no' then write_png, fname_g_prefix+'plot2'+fname_g_suffix, tvrd(), red, green, blue
if window_2 ne 'no' then spawn, 'convert '+fname_g_prefix+'plot2'+fname_g_suffix+' '+fname_g_prefix+'plot2.gif'


;endif ; if window_2 ne 'no' then begin
if window_3 ne 'no' then begin

if window_3 ne 'no' then  WINDOW, 3, XSIZE=1000, YSIZE=1000, TITLE='ARM Integrated Column Product Plots', /pixmap

;print, window_3, ' is window fucking 3'

!p.multi = [0,6,6]
;height=height/1000.
!p.charsize=2.0



;solar flux down

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_down[kk,jj] gt 0. and solar_flux_down[kk,jj] lt 1500. then image[j,k]=solar_flux_down[kk,jj] $
			else if solar_flux_down[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor


xl=0.05 & yl=0.825 & xr=0.45 & yr=0.925	; upper left


if window_3 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'DownWelling Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',200.+max(image[where(image lt 1500.)]), min(image[where(image gt 0.0)])

if window_3 ne 'no' then xyouts, 0.5,0.97, title_string, /normal, alignment=0.5, color=0

; Downwelling IR flux

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IR_flux_down[kk,jj] gt 0. and IR_flux_down[kk,jj] lt 1500. then image[j,k]=IR_flux_down[kk,jj] $
			else if IR_flux_down[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.825 & xr=0.935 & yr=0.925	; upper right


if window_3 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'DownWelling IR Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image[where(image lt 1500.)])+200., min(image[where(image gt 0.0)])

; Upwelling Solar

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_up[kk,jj] gt 0. and solar_flux_up[kk,jj] lt 1500. then image[j,k]=solar_flux_up[kk,jj] $
			else if solar_flux_up[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.675 & xr=0.45 & yr=0.775 ; middle left


if window_3 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'UpWelling Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image[where(image lt 1500.)])+200., min(image[where(image gt 0.0)])

; Upwelling IR

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IR_flux_up[kk,jj] gt 0. and IR_flux_up[kk,jj] lt 1500. then image[j,k]=IR_flux_up[kk,jj] $
			else if IR_flux_up[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.675 & xr=0.935 & yr=0.775 ; middle right


if window_3 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'UpWelling IR Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image[where(image lt 1500.)])+200., min(image[where(image gt 0.0)])

; Net Solar down

image_solar_net=image
density=image_solar_net
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_up[kk,jj] gt 0. $
			and solar_flux_up[kk,jj] lt 1500. $
			and solar_flux_down[kk,jj] gt 0. and $
			solar_flux_down[kk,jj] lt 1500. then begin
				image_solar_net[j,k]=(solar_flux_down[kk,jj]-solar_flux_up[kk,jj])
				density[j,k]=pressure[kk,jj]/(temp[kk,jj]*287.04)
			endif else begin
				image_solar_net[j,k]=-9999.
				density[j,k]=pressure[kk,jj]/(temp[kk,jj]*287.04)
			endelse

			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.525 & xr=0.45 & yr=0.625 ; middle left 2

if n_elements((where(image_solar_net lt 1500. and image_solar_net gt -1500.))) eq 1  then window_3='no'
print, max(image_solar_net), min(image_solar_net), n_elements((where(image_solar_net lt 1500. and image_solar_net gt -1500.)))

if window_3 ne 'no' then plot_2d_image_999_maxminin, image_solar_net, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Net Solar Down', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image_solar_net)+200., min(image_solar_net)

; Net IR down

image_ir_net=image
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IR_flux_up[kk,jj] gt 0. $
			and IR_flux_up[kk,jj] lt 1500. $
			and IR_flux_down[kk,jj] gt 0. and $
			IR_flux_down[kk,jj] lt 1500. then begin
				image_ir_net[j,k]=(IR_flux_down[kk,jj]-IR_flux_up[kk,jj])
			endif else begin
				image_ir_net[j,k]=-9999.
			endelse

			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.525 & xr=0.935 & yr=0.625 ; middle right 2


if window_3 ne 'no' and n_elements(where(image_ir_net gt -1500.)) gt 1 then plot_2d_image_999_maxminin, image_ir_net, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Net IR Down', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image_ir_net[where(image_ir_net lt 1500.)])+200., min(image_ir_net[where(image_ir_net gt -1500.)])

; Solar Heating

solar_heating=image
solar_heating[*,*]=0.

for j=0,n_elements(plot_xaxis)-1 do begin

	for k=1,n_elements(plot_yaxis)-2 do begin

		if image_solar_net[j,k-1] ne -9999. and image_solar_net[j,k+1] ne -9999. and density[j,k] lt 2.0 and $
		density[j,k] gt 0.0 then begin


			solar_heating[j,k]=((86400.)/(density[j,k]*1004.))*((image_solar_net[j,k+1]-image_solar_net[j,k-1])/180.)


		endif


	endfor

endfor

xl=0.05 & yl=0.375 & xr=0.45 & yr=0.475 ; middle left 2


if window_3 ne 'no' then plot_2d_image_999_maxminin, solar_heating, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Solar Heating (k/day)', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'K/day',10.+max(solar_heating[where(solar_heating lt 100.)]), min(solar_heating[where(solar_heating gt -100.)])


; IR Heating

IR_heating=image
IR_heating[*,*]=-9999.
net_heating=IR_heating

for j=0,n_elements(plot_xaxis)-1 do begin

	for k=1,n_elements(plot_yaxis)-2 do begin

		if image_ir_net[j,k-1] ne -9999. and image_ir_net[j,k+1] ne -9999. and density[j,k] lt 2.0 and $
		density[j,k] gt 0.0 then begin


			IR_heating[j,k]=((86400.)/(density[j,k]*1004.))*((image_ir_net[j,k+1]-image_ir_net[j,k-1])/180.)
			net_heating[j,k]=((86400.)/(density[j,k]*1004.))*(((image_ir_net[j,k+1]-image_ir_net[j,k-1])/180.)+((image_solar_net[j,k+1]-image_solar_net[j,k-1])/180.))

		endif


	endfor

endfor

xl=0.535 & yl=0.375 & xr=0.935 & yr=0.475 ; middle right 2


if window_3 ne 'no' and n_elements(where(ir_heating gt -50.)) gt 1 then plot_2d_image_999_maxminin, IR_heating, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'IR Heating (k/day)', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'K/day',max(ir_heating[where(ir_heating lt 100.)]), min(ir_heating[where(ir_heating gt -50.)])


xl=0.05 & yl=0.225 & xr=0.45 & yr=0.325 ; middle left 2


if window_3 ne 'no' and n_elements(where(net_heating gt -100.)) gt 1 then plot_2d_image_999_maxminin, net_heating, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Net Heating (k/day)', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'K/day',max(net_heating[where(net_heating lt 100.)]), min(net_heating[where(net_heating gt -100.)])


xl=0.535 & yl=0.225 & xr=0.935 & yr=0.325 ; middle right 2

; broadbnd LW flux - Minnis vs Retrieval comparison

if window_3 ne 'no' then plot, plot_xaxis, total_visible_optical_depth, psym=-3, pos=five_right_pos, ytitle='Tau', $
	title='Visible Optical Depth (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[0.0, max(total_visible_optical_depth)]

if n_elements(Visible_Optical_Depth) gt 0 then begin

; find the coordinates of the site

site_lat=36.6 & lat_index=0 & lat_index=where(abs(site_lat-latitude_M) eq min(abs(site_lat-latitude_M)) )
site_lon=-97.5 & lon_index=0 & lon_index=where(abs(site_lon-longitude_M) eq min(abs(site_lon-longitude_M)) )

if window_3 ne 'no' then oplot, hrfrac_M, Visible_Optical_Depth[lon_index[0],lat_index[0],3,*], psym=2, linestyle=2, color=0

endif


xl=0.05 & yl=0.075 & xr=0.45 & yr=0.175 ; middle left 2


vector_olr=fltarr(n_elements(plot_xaxis)) & vector_albedo=vector_olr
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

	vector_olr[j]=TOA_IR_up[jj]
	vector_albedo[j]=(TOA_solar_up[jj]/TOA_solar_down[jj])

jj=jj+1
endfor

if window_3 ne 'no' then plot, plot_xaxis, vector_olr, psym=-3, pos=six_left_pos, ytitle='OLR', $
	title='Outgoing Longwave Flux (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[100.0, 400.]

if n_elements(Visible_Optical_Depth) gt 0 then begin

if window_3 ne 'no' then  oplot, hrfrac_M, Broadband_LW_Flux[lon_index[0],lat_index[0],1,*], psym=2, linestyle=2, color=0

	;print, hrfrac_M[where(hrfrac_M gt 17 and hrfrac_M lt 18)],Broadband_LW_Flux[lon_index[0],lat_index[0],1,where(hrfrac_M gt 17 and hrfrac_M lt 18)]

endif

if window_3 ne 'no' then plot, plot_xaxis, vector_albedo*100., psym=-3, pos=six_right_pos, ytitle='OLR', $
	title='Broadband Albedo (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[0.0, 100.]

if n_elements(Visible_Optical_Depth) gt 0 then begin

if window_3 ne 'no' then oplot, hrfrac_M, Broadband_SW_Albedo[lon_index[0],lat_index[0],1,*], psym=2, linestyle=2, color=0

endif

if window_3 ne 'no' then write_png, fname_g_prefix+'plot3'+fname_g_suffix, tvrd(), red, green, blue
if window_3 ne 'no' then spawn, 'convert '+fname_g_prefix+'plot3'+fname_g_suffix+' '+fname_g_prefix+'plot3.gif'


endif ; if window_3 ne 'no' then begin

;;;;;;;;;;;;;;;;;;;;; window 4 - clear stuff

if window_4 ne 'no' then begin

if window_4 ne 'no' then WINDOW, 4, XSIZE=1000, YSIZE=1000, TITLE='ARM Integrated Column Product Plots', /pixmap

!p.multi = [0,6,6]
;height=height/1000.
!p.charsize=2.0



;solar flux down

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_down_clear[kk,jj] gt 0. and solar_flux_down_clear[kk,jj] lt 1500. then image[j,k]=solar_flux_down_clear[kk,jj] $
			else if solar_flux_down_clear[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor


xl=0.05 & yl=0.825 & xr=0.45 & yr=0.925	; upper left


if window_4 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'DownWelling Clear Sky Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',200.+max(image[where(image lt 1500.)]), min(image[where(image gt 0.0)])

if window_4 ne 'no' then xyouts, 0.5,0.97, title_string, /normal, alignment=0.5, color=0

; Downwelling IR flux

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IR_flux_down_clear[kk,jj] gt 0. and IR_flux_down_clear[kk,jj] lt 1500. then image[j,k]=IR_flux_down_clear[kk,jj] $
			else if IR_flux_down_clear[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.825 & xr=0.935 & yr=0.925	; upper right


if window_4 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'DownWelling Clear Sky IR Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image[where(image lt 1500.)])+200., min(image[where(image gt 0.0)])

; Upwelling Solar

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_up_clear[kk,jj] gt 0. and solar_flux_up_clear[kk,jj] lt 1500. then image[j,k]=solar_flux_up_clear[kk,jj] $
			else if solar_flux_up_clear[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.675 & xr=0.45 & yr=0.775 ; middle left


if window_4 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'UpWelling Clear Sky Solar Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image[where(image lt 1500.)])+200., min(image[where(image gt 0.0)])

; Upwelling IR

jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IR_flux_up_clear[kk,jj] gt 0. and IR_flux_up_clear[kk,jj] lt 1500. then image[j,k]=IR_flux_up_clear[kk,jj] $
			else if IR_flux_up_clear[kk,jj] le 0. then image[j,k]=-9999.
			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.675 & xr=0.935 & yr=0.775 ; middle right


if window_4 ne 'no' then plot_2d_image_999_maxminin, image, plot_xaxis, plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'UpWelling Clear Sky IR Flux', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',max(image[where(image lt 1500.)])+200., min(image[where(image gt 0.0)])

; Net Solar down

image_solar_net_clear=image
density=image_solar_net_clear
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if solar_flux_up_clear[kk,jj] gt 0. $
			and solar_flux_up_clear[kk,jj] lt 1500. $
			and solar_flux_down_clear[kk,jj] gt 0. and $
			solar_flux_down_clear[kk,jj] lt 1500. then begin
				image_solar_net_clear[j,k]=(solar_flux_down_clear[kk,jj]-solar_flux_up_clear[kk,jj])
				density[j,k]=pressure[kk,jj]/(temp[kk,jj]*287.04)
			endif else begin
				image_solar_net_clear[j,k]=-9999.
				density[j,k]=pressure[kk,jj]/(temp[kk,jj]*287.04)
			endelse

			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.05 & yl=0.525 & xr=0.45 & yr=0.625 ; middle left 2


if window_4 ne 'no' then plot_2d_image_999_maxminin, image_solar_net_clear, plot_xaxis, $
plot_yaxis/1000., n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Net Clear Sky Solar Down', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',$
	max(image_solar_net_clear[where(image_solar_net_clear lt 1500.)])+200., $
	min(image_solar_net_clear[where(image_solar_net_clear gt -1500.)])

; Net IR down

image_ir_net_clear=image
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin
kk=start_height_index
	for k=0,n_elements(plot_yaxis)-1 do begin
			if IR_flux_up_clear[kk,jj] gt 0. $
			and IR_flux_up_clear[kk,jj] lt 1500. $
			and IR_flux_down_clear[kk,jj] gt 0. and $
			IR_flux_down_clear[kk,jj] lt 1500. then begin
				image_ir_net_clear[j,k]=(IR_flux_down_clear[kk,jj]-IR_flux_up_clear[kk,jj])
			endif else begin
				image_ir_net_clear[j,k]=-9999.
			endelse

			kk=kk+1
	endfor
jj=jj+1
endfor

xl=0.535 & yl=0.525 & xr=0.935 & yr=0.625 ; middle right 2


if window_4 ne 'no' then plot_2d_image_999_maxminin, image_ir_net_clear, plot_xaxis, plot_yaxis/1000., $
n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Net Clear Sky IR Down', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'W/m^2',$
	max(image_ir_net_clear[where(image_ir_net_clear lt 1500.)])+200., $
	min(image_ir_net_clear[where(image_ir_net_clear gt -1500.)])

; Solar Heating

solar_heating_clear=image
solar_heating_clear[*,*]=0.

for j=0,n_elements(plot_xaxis)-1 do begin

	for k=1,n_elements(plot_yaxis)-2 do begin

		if image_solar_net_clear[j,k-1] ne -9999. and image_solar_net_clear[j,k+1] ne -9999. and density[j,k] lt 2.0 and $
		density[j,k] gt 0.0 then begin


			solar_heating_clear[j,k]=((86400.)/(density[j,k]*1004.))*((image_solar_net_clear[j,k+1]-image_solar_net_clear[j,k-1])/180.)


		endif


	endfor

endfor

xl=0.05 & yl=0.375 & xr=0.45 & yr=0.475 ; middle left 2


if window_4 ne 'no' then plot_2d_image_999_maxminin, solar_heating_clear, plot_xaxis, plot_yaxis/1000., $
n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Solar Clear Sky Heating (k/day)', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'K/day',$
	10.+max(solar_heating_clear[where(solar_heating_clear lt 100.)]), $
	min(solar_heating_clear[where(solar_heating_clear gt -100.)])


; IR Heating

IR_heating_clear=image
IR_heating_clear[*,*]=-9999.
net_heating_clear=IR_heating_clear

for j=0,n_elements(plot_xaxis)-1 do begin

	for k=1,n_elements(plot_yaxis)-2 do begin

		if image_ir_net_clear[j,k-1] ne -9999. and image_ir_net_clear[j,k+1] ne -9999. and density[j,k] lt 2.0 and $
		density[j,k] gt 0.0 then begin


			IR_heating_clear[j,k]=((86400.)/(density[j,k]*1004.))*((image_ir_net_clear[j,k+1]-image_ir_net_clear[j,k-1])/180.)
			net_heating_clear[j,k]=((86400.)/(density[j,k]*1004.))*(((image_ir_net_clear[j,k+1]-image_ir_net_clear[j,k-1])/180.)+((image_solar_net_clear[j,k+1]-image_solar_net_clear[j,k-1])/180.))

		endif


	endfor

endfor

xl=0.535 & yl=0.375 & xr=0.935 & yr=0.475 ; middle right 2


if window_4 ne 'no' then plot_2d_image_999_maxminin, IR_heating_clear, plot_xaxis, plot_yaxis/1000., $
n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'IR Clear Sky Heating (k/day)', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'K/day',$
	max(ir_heating_clear[where(ir_heating_clear lt 100.)]), $
	min(ir_heating_clear[where(ir_heating_clear gt -50.)])


xl=0.05 & yl=0.225 & xr=0.45 & yr=0.325 ; middle left 2


if window_4 ne 'no' then plot_2d_image_999_maxminin, net_heating_clear, plot_xaxis, plot_yaxis/1000., $
n_elements(plot_xaxis), n_elements(plot_yaxis), $
	 'Net Heating (k/day)', ' ', 'Height (km)', [xl,yl,xr,yr], $
	[xr+0.045,yl,xr+0.05,yr], 'K/day',$
	max(net_heating_clear[where(net_heating_clear lt 100.)]), $
	min(net_heating_clear[where(net_heating_clear gt -100.)])


xl=0.535 & yl=0.225 & xr=0.935 & yr=0.325 ; middle right 2

; broadbnd LW flux - Minnis vs Retrieval comparison

;plot, plot_xaxis, total_visible_optical_depth, psym=-3, pos=five_right_pos, ytitle='Tau', $
;	title='Visible Optical Depth (GOES -> *)',background=!d.n_colors-1, $
;	color=0, yrange=[0.0, max(total_visible_optical_depth)]

;if n_elements(Visible_Optical_Depth) gt 0 then begin

; find the coordinates of the site

;site_lat=36.6 & lat_index=0 & lat_index=where(abs(site_lat-latitude_M) eq min(abs(site_lat-latitude_M)) )
;site_lon=-97.5 & lon_index=0 & lon_index=where(abs(site_lon-longitude_M) eq min(abs(site_lon-longitude_M)) )

;oplot, hrfrac_M, Visible_Optical_Depth[lon_index[0],lat_index[0],3,*], psym=2, linestyle=2, color=0

;endif


xl=0.05 & yl=0.075 & xr=0.45 & yr=0.175 ; middle left 2


vector_olr_clear=fltarr(n_elements(plot_xaxis)) & vector_albedo_clear=vector_olr_clear
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

	vector_olr_clear[j]=TOA_IR_up_clear[jj]
	vector_albedo_clear[j]=(TOA_solar_up_clear[jj]/TOA_solar_down_clear[jj])

jj=jj+1
endfor

if window_4 ne 'no' then plot, plot_xaxis, vector_olr_clear, psym=-3, pos=six_left_pos, ytitle='OLR', $
	title='Clear Sky Outgoing Longwave Flux (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[100.0, 400.]

if n_elements(Visible_Optical_Depth) gt 0 then begin

; find the coordinates of the site

site_lat=36.6 & lat_index=0 & lat_index=where(abs(site_lat-latitude_M) eq min(abs(site_lat-latitude_M)) )
site_lon=-97.5 & lon_index=0 & lon_index=where(abs(site_lon-longitude_M) eq min(abs(site_lon-longitude_M)) )

if window_4 ne 'no' then oplot, hrfrac_M, Broadband_LW_Flux[lon_index[0],lat_index[0],1,*], psym=2, linestyle=2, color=0

endif

if window_4 ne 'no' then plot, plot_xaxis, vector_albedo_clear*100., psym=-3, pos=six_right_pos, ytitle='OLR', $
	title='Broadband Albedo (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[0.0, 100.]

if n_elements(Visible_Optical_Depth) gt 0 then begin

if window_4 ne 'no' then oplot, hrfrac_M, Broadband_SW_Albedo[lon_index[0],lat_index[0],1,*], psym=2, linestyle=2, color=0

endif


if window_4 ne 'no' then write_png, fname_g_prefix+'plot4'+fname_g_suffix, tvrd(), red, green, blue
if window_4 ne 'no' then spawn, 'convert '+fname_g_prefix+'plot4'+fname_g_suffix+' '+fname_g_prefix+'plot4.gif'


endif ; if window_4 ne 'no' then begin


; Now we want to plot the cloud forcing values,  TOA:  this will be the cloudy sky
; net minus the clear sky net solar, IR and total


;;;;;;;;;;;;;;;;;;;;; window 4 - clear stuff

if window_5 ne 'no' then begin

if window_5 ne 'no' then WINDOW, 5, XSIZE=1000, YSIZE=1000, TITLE='ARM Integrated Column Product Plots', /pixmap

!p.multi = [0,6,6]
;height=height/1000.
!p.charsize=2.0



;TOA solar cloud forcing

vector_scf_toa=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if TOA_solar_down[jj] gt 0 and TOA_solar_down[jj] lt 1500. $
		and TOA_solar_up[jj] gt 0 and TOA_solar_up[jj] lt 1500. then begin
	vector_scf_toa[j]=(TOA_solar_down[jj]-TOA_solar_up[jj])-(TOA_solar_down_clear[jj]-TOA_solar_up_clear[jj])	; down into the system is positive
		endif else begin
	vector_scf_toa[j]=-9999.
		endelse
jj=jj+1
endfor

if n_elements((where(vector_scf_toa lt 1500. and vector_scf_toa gt -1500.))) eq 1 then window_5='no'

;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then plot, plot_xaxis, vector_scf_toa , psym=-3, pos=one_left_pos, ytitle='Shortwave CRF (w/m^2)', $
	title='TOA Shortwave Cloud Radiative Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_scf_toa[where(vector_scf_toa gt -1500.)])-10., $
	10.+max(vector_scf_toa[where(vector_scf_toa lt 100.)])]


if n_elements(Visible_Optical_Depth) gt 0 then begin

;find the closest point in time of the arm data with the satellite data
vector_scf_toa_M=fltarr(n_elements(Visible_Optical_Depth[1,1,1,*]))
for kk=0,n_elements(hrfrac_M)-1 do begin

time_index=where( abs(hrfrac-hrfrac_M[kk]) eq min(abs(hrfrac-hrfrac_M[kk])))

if hrfrac_M[kk] gt 17. and hrfrac_M[kk] lt 18. then print, (TOA_solar_down[time_index[0]]*(Broadband_SW_Albedo[lon_index[0],lat_index[0],1,kk]/100.))

		if TOA_solar_down[time_index[0]] gt 0.0 and TOA_solar_down[time_index[0]] lt 1500. $
		and Broadband_SW_Albedo[lon_index[0],lat_index[0],1,kk] gt 0.0 $
		and Broadband_SW_Albedo[lon_index[0],lat_index[0],1,kk] le 100.0 $
		and TOA_solar_down[time_index[0]] gt 0.0 and TOA_solar_down[time_index[0]] lt 1500. $
		and TOA_solar_down_clear[time_index[0]] gt 0.0 and TOA_solar_up_clear[time_index[0]] gt 0.0 then begin
	vector_scf_toa_M[kk]=	(TOA_solar_down[time_index[0]]-$
	(TOA_solar_down[time_index[0]]*(Broadband_SW_Albedo[lon_index[0],lat_index[0],1,kk]/100.)))-$
	(TOA_solar_down_clear[time_index[0]]-(TOA_solar_down[time_index[0]]*(Broadband_SW_Albedo[lon_index[0],lat_index[0],0,kk]/100.)))
	;(TOA_solar_down_clear[time_index[0]]-TOA_solar_up_clear[time_index[0]])

if hrfrac[time_index[0]] gt 13. and hrfrac[time_index[0]] lt 14. then begin
	print, 'here'
endif

		endif else begin
	vector_scf_toa_M[kk]=-9999.
		endelse

endfor

if window_5 ne 'no' then oplot, hrfrac_M, vector_scf_toa_M, psym=2, linestyle=2, color=0

endif

;;;;; TOA Longwave Cloud Radiative Forcing

vector_lcf_toa=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if TOA_IR_up[jj] gt 0 and TOA_IR_up[jj] lt 1500. then begin
	vector_lcf_toa[j]=(TOA_IR_up_clear[jj])-(TOA_IR_up[jj])	; down into the system is positive
		endif else begin
	vector_lcf_toa[j]=-9999.
		endelse
jj=jj+1
endfor


;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then plot, plot_xaxis, vector_lcf_toa , psym=-3, pos=one_right_pos, ytitle='Longwave CRF (w/m^2)', $
	title='TOA Longwave Cloud Radiative Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_lcf_toa[where(vector_lcf_toa gt -1500.)])-10., $
	10.+max(vector_lcf_toa[where(vector_lcf_toa lt 100.)])]


if n_elements(Visible_Optical_Depth) gt 0 then begin

;find the closest point in time of the arm data with the satellite data
vector_lcf_toa_M=fltarr(n_elements(Visible_Optical_Depth[1,1,1,*]))
for kk=0,n_elements(hrfrac_M)-1 do begin

time_index=where( abs(hrfrac-hrfrac_M[kk]) eq min(abs(hrfrac-hrfrac_M[kk])))

		if Broadband_LW_Flux[lon_index[0],lat_index[0],1,kk] gt 0.0 $
		and Broadband_LW_Flux[lon_index[0],lat_index[0],1,kk] lt 1500. and $
		TOA_IR_up_clear[time_index[0]] gt 0.0 and TOA_IR_up_clear[time_index[0]] lt 1500. then begin

	vector_lcf_toa_M[kk]=Broadband_LW_Flux[lon_index[0],lat_index[0],0,kk]-Broadband_LW_Flux[lon_index[0],lat_index[0],1,kk]

		endif else begin
	vector_lcf_toa_M[kk]=-9999.
		endelse

endfor

if window_5 ne 'no' then oplot, hrfrac_M, vector_lcf_toa_M, psym=2, linestyle=2, color=0

endif

; net TOA CRF

vector_netcf_toa=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if vector_lcf_toa[j] gt -9999. and vector_scf_toa[j] gt -9999. then begin
	vector_netcf_toa[j]=vector_lcf_toa[j]+vector_scf_toa[j]	; down into the system is positive
			endif else begin
			if vector_lcf_toa[j] gt -9999. and vector_scf_toa[j] le -9999. then begin
			vector_netcf_toa[j]=vector_lcf_toa[j]
			endif else begin
	vector_netcf_toa[j]=-9999.
			endelse
		endelse
jj=jj+1
endfor


;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then plot, plot_xaxis, vector_netcf_toa , psym=-3, pos=two_left_pos, ytitle='Net CRF (w/m^2)', $
	title='TOA Net Cloud Radiative Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_netcf_toa[where(vector_netcf_toa gt -1500.)])-10., $
	10.+max(vector_netcf_toa[where(vector_netcf_toa lt 100.)])]


if n_elements(Visible_Optical_Depth) gt 0 then begin

;find the closest point in time of the arm data with the satellite data
vector_netcf_toa_M=fltarr(n_elements(Visible_Optical_Depth[1,1,1,*]))
for kk=0,n_elements(hrfrac_M)-1 do begin

		if vector_lcf_toa_M[kk] gt -9999. and vector_scf_toa_M[kk] gt -9999. $
		and Solar_Zenith_Angle[lon_index[0],lat_index[0],kk] gt 0.$
		and Solar_Zenith_Angle[lon_index[0],lat_index[0],kk] lt 70. then begin

	vector_netcf_toa_M[kk]=vector_lcf_toa_M[kk]+vector_scf_toa_M[kk]

		endif else begin

		if vector_lcf_toa_M[kk] gt -9999. and vector_scf_toa_M[kk] le -9999. $
		and Solar_Zenith_Angle[lon_index[0],lat_index[0],kk] le 0. then begin
			vector_netcf_toa_M[kk]=vector_lcf_toa_M[kk]
		endif else begin
	vector_netcf_toa_M[kk]=-9999.
		endelse
		endelse

endfor

if window_5 ne 'no' then oplot, hrfrac_M, vector_netcf_toa_M, psym=2, linestyle=2, color=0

endif


;;;;; Surface net Longwave radiation

cloud_free_flag=intarr(n_elements(reflectivity[0,*])) & cloud_free_flag[*]=1
for j=0,n_elements(cloud_free_flag)-1 do begin
	if max(reflectivity[*,j]) gt -9999. then cloud_free_flag[j]=0
endfor

vector_net_longwave_sfc=fltarr(n_elements(plot_xaxis))
vector_net_longwave_sfc_obs=fltarr(n_elements(plot_xaxis))
vector_net_longwave_sfc_clr=fltarr(n_elements(plot_xaxis))
vector_net_longwave_sfc_obs_clr=fltarr(n_elements(plot_xaxis))
longwave_sfc_cld_forcing=fltarr(n_elements(plot_xaxis))
longwave_sfc_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
;longwave_sfc_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
IR_offset=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if IR_flux_up[0,jj] gt 0. and IR_flux_up[0,jj] lt 1500. $
		and IR_flux_down[0,jj] and IR_flux_down[0,jj] lt 1500. then begin
	vector_net_longwave_sfc[j]=(IR_flux_down[1,jj])-(IR_flux_up[1,jj])	; down into the system is positive
		endif else begin
	vector_net_longwave_sfc[j]=-9999.
		endelse

		if IR_flux_up_clear[0,jj] gt 0. and IR_flux_up_clear[0,jj] lt 1500. $
		and IR_flux_down_clear[0,jj] and IR_flux_down_clear[0,jj] lt 1500. then begin
	vector_net_longwave_sfc_clr[j]=(IR_flux_down_clear[1,jj])-(IR_flux_up_clear[1,jj])	; down into the system is positive
		endif else begin
	vector_net_longwave_sfc_clr[j]=-9999.
		endelse



	;	if down_long[jj] gt 0. and up_long[jj] gt 0. and sfc_temp[jj] gt 0. and sfc_mxrat[jj] gt 0. and $
	;		vap[jj] gt 0 then begin



	;down_long_surface_clear=mlr_coeff_down[0]+(mlr_coeff_down[1]*sfc_temp[jj])+ $
	;	(mlr_coeff_down[2]*sfc_mxrat[jj])+(mlr_coeff_down[3]*hrfrac[jj])+ $
	;	(mlr_coeff_down[4]*vap[jj])

	down_long_surface_clear=IR_flux_down_clear[1,jj]*(down_long[jj]/IR_flux_down[1,jj])

	;up_long_surface_clear=mlr_coeff_up[0]+(mlr_coeff_up[1]*sfc_temp[jj])+ $
	;	(mlr_coeff_up[2]*sfc_mxrat[jj])+(mlr_coeff_up[3]*hrfrac[jj])+ $
	;	(mlr_coeff_up[4]*vap[jj])

	up_long_surface_clear=IR_flux_up_clear[1,jj]*(up_long[jj]/IR_flux_up[1,jj])

	vector_net_longwave_sfc_obs_clr[j]=down_long_surface_clear-up_long_surface_clear	; down into the system is positive
	;	endif else begin
	;vector_net_longwave_sfc_obs_clr[j]=-9999.
	;	endelse

	;vector_net_longwave_sfc_obs_clr[j]=vector_net_longwave_sfc_clr[j]-20.

		if down_long[jj] gt 0. and down_long[jj] lt 1500. $
		and up_long[jj] and up_long[jj] lt 1500. then begin
	vector_net_longwave_sfc_obs[j]=(down_long[jj])-(up_long[jj])	; down into the system is positive
		endif else begin
	vector_net_longwave_sfc_obs[j]=-9999.
		endelse


	if vector_net_longwave_sfc[j] gt -9999. and vector_net_longwave_sfc_clr[j] gt -9999. then begin
		longwave_sfc_cld_forcing[j]=vector_net_longwave_sfc[j]-vector_net_longwave_sfc_clr[j]
	endif else begin
		longwave_sfc_cld_forcing[j]=-9999.
	endelse


	if vector_net_longwave_sfc_obs_clr[j] gt -9999. and vector_net_longwave_sfc_obs_clr[j] gt -9999. then begin
		longwave_sfc_cld_forcing_obs[j]=vector_net_longwave_sfc_obs[j]-vector_net_longwave_sfc_obs_clr[j]
	endif else begin
		longwave_sfc_cld_forcing_obs[j]=-9999.
	endelse

;if max(cloud_free_flag) eq 1 then begin

;	index=where(cloud_free_flag eq 1)
;	index2=where(abs(hrfrac[index]-plot_xaxis[j]) eq min(abs(hrfrac[index]-plot_xaxis[j])))


;	IR_offset[jj]=vector_net_longwave_sfc_obs[index[index2[0]]]-vector_net_longwave_sfc[index[index2[0]]]

;	if abs(vector_net_longwave_sfc[jj]+IR_offset[jj]-vector_net_longwave_sfc_obs[jj]) gt abs(vector_net_longwave_sfc[jj]-20.0-vector_net_longwave_sfc_obs[jj]) then IR_offset[jj]=-20.

;endif else begin

;	IR_offset[j]=-20.

;endelse

jj=jj+1
endfor



;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then begin
	plot, plot_xaxis, vector_net_longwave_sfc , psym=-3, pos=two_right_pos, ytitle='Longwave CRF (w/m^2)', $
	title='Sfc Net Longwave Radiation (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_net_longwave_sfc[where(vector_net_longwave_sfc gt -1500.)]), $
	10.+max(vector_net_longwave_sfc[where(vector_net_longwave_sfc lt 1000.)])]

	oplot, plot_xaxis, vector_net_longwave_sfc_obs, psym=3, linestyle=2, color=35

endif

;;;;; Surface Longwave down radiation


vector_down_longwave_sfc=fltarr(n_elements(plot_xaxis))
vector_down_longwave_sfc_obs=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if IR_flux_down[0,jj] and IR_flux_down[0,jj] lt 1500. then begin
	vector_down_longwave_sfc[j]=(IR_flux_down[0,jj])	; down into the system is positive
		endif else begin
	vector_down_longwave_sfc[j]=-9999.
		endelse

		if down_long[jj] gt 0. and down_long[jj] lt 1500. then begin
	vector_down_longwave_sfc_obs[j]=(down_long[jj])	; down into the system is positive
		endif else begin
	vector_down_longwave_sfc_obs[j]=-9999.
		endelse

jj=jj+1
endfor



;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then begin
	plot, plot_xaxis, vector_down_longwave_sfc , psym=-3, pos=three_left_pos, ytitle='sfc longwave down (w/m^2)', $
	title='Sfc Longwave Down (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_down_longwave_sfc[where(vector_down_longwave_sfc gt -1500.)]), $
	10.+max(vector_down_longwave_sfc[where(vector_down_longwave_sfc lt 1000.)])]

	oplot, plot_xaxis, vector_down_longwave_sfc_obs, psym=3, linestyle=2, color=35

endif

;;;;; Surface Longwave up radiation


vector_up_longwave_sfc=fltarr(n_elements(plot_xaxis))
vector_up_longwave_sfc_obs=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if IR_flux_up[0,jj] gt 0. and IR_flux_up[0,jj] lt 1500. then begin
	vector_up_longwave_sfc[j]=(IR_flux_up[0,jj])	; down into the system is positive
		endif else begin
	vector_up_longwave_sfc[j]=-9999.
		endelse

		if up_long[jj] gt 0. and up_long[jj] lt 1500. then begin
	vector_up_longwave_sfc_obs[j]=(up_long[jj])	; down into the system is positive
		endif else begin
	vector_up_longwave_sfc_obs[j]=-9999.
		endelse

jj=jj+1
endfor

if window_5 ne 'no' then begin
	plot, plot_xaxis, vector_up_longwave_sfc, psym=-3, pos=three_right_pos, ytitle='sfc longwave up (w/m^2)', $
	title='Sfc Longwave up (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_up_longwave_sfc_obs[where(vector_up_longwave_sfc_obs gt -1500.)]), $
	10.+max(vector_up_longwave_sfc[where(vector_up_longwave_sfc lt 1000.)])]

	oplot, plot_xaxis, vector_up_longwave_sfc_obs, psym=3, linestyle=2, color=35

endif

; plot the longwave cloud forcing


if window_5 ne 'no' then begin
	plot, plot_xaxis, longwave_sfc_cld_forcing , psym=-3, pos=four_left_pos, ytitle='sfc longwave cloud forcing (w/m^2)', $
	title='Sfc Longwave cld Forcing (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(longwave_sfc_cld_forcing[where(longwave_sfc_cld_forcing gt -1500.)]), $
	10.+max(longwave_sfc_cld_forcing[where(longwave_sfc_cld_forcing lt 1000.)])]

	oplot, plot_xaxis, longwave_sfc_cld_forcing_obs, psym=-3, linestyle=2, color=35

endif
;;;;;;;;;;;;;;;;;;;;;;; shortwave

;;;;; Surface net shortwave radiation

vector_net_shortwave_sfc=fltarr(n_elements(plot_xaxis))
vector_net_shortwave_sfc_obs=fltarr(n_elements(plot_xaxis))
solar_offset=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if solar_flux_up[0,jj] gt 0. and solar_flux_up[0,jj] lt 1500. $
		and solar_flux_down[0,jj] and solar_flux_down[0,jj] lt 1500. then begin
	vector_net_shortwave_sfc[j]=(solar_flux_down[0,jj])-(solar_flux_up[0,jj])	; down into the system is positive
		endif else begin
	vector_net_shortwave_sfc[j]=-9999.
		endelse

		if down_short[jj] gt 0. and down_short[jj] lt 1500. $
		and up_short[jj] and up_short[jj] lt 1500. then begin
	vector_net_shortwave_sfc_obs[j]=(down_short[jj])-(up_short[jj])	; down into the system is positive
		endif else begin
	vector_net_shortwave_sfc_obs[j]=-9999.
		endelse

jj=jj+1
endfor



;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then begin
	plot, plot_xaxis, vector_net_shortwave_sfc , psym=-3, pos=four_right_pos, ytitle='shortwave CRF (w/m^2)', $
	title='Sfc Net shortwave Radiation (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_net_shortwave_sfc[where(vector_net_shortwave_sfc gt -1500.)]), $
	10.+max(vector_net_shortwave_sfc[where(vector_net_shortwave_sfc lt 1000.)])]

	oplot, plot_xaxis, vector_net_shortwave_sfc_obs, psym=3, linestyle=2, color=35

endif

;;;;; Surface shortwave down radiation


vector_down_shortwave_sfc=fltarr(n_elements(plot_xaxis))
vector_down_shortwave_sfc_obs=fltarr(n_elements(plot_xaxis))
vector_down_shortwave_sfc_clr=fltarr(n_elements(plot_xaxis))
vector_down_shortwave_sfc_obs_clr=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if solar_flux_down[0,jj] and solar_flux_down[0,jj] lt 1500. then begin
	vector_down_shortwave_sfc[j]=(solar_flux_down[0,jj])	; down into the system is positive
		endif else begin
	vector_down_shortwave_sfc[j]=-9999.
		endelse

		if solar_flux_down_clear[0,jj] and solar_flux_down_clear[0,jj] lt 1500. then begin
	vector_down_shortwave_sfc_clr[j]=(solar_flux_down_clear[0,jj])	; down into the system is positive
		endif else begin
	vector_down_shortwave_sfc_clr[j]=-9999.
		endelse

		if down_short[jj] gt 0. and down_short[jj] lt 1500. then begin
	vector_down_shortwave_sfc_obs[j]=(down_short[jj])	; down into the system is positive
		endif else begin
	vector_down_shortwave_sfc_obs[j]=-9999.
		endelse

		if fluxclear[jj] gt 0. and fluxclear[jj] lt 1500. then begin
	vector_down_shortwave_sfc_obs_clr[j]=(fluxclear[jj])	; down into the system is positive
		endif else begin
	vector_down_shortwave_sfc_obs_clr[j]=-9999.
		endelse

jj=jj+1
endfor



;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then begin
	plot, plot_xaxis, vector_down_shortwave_sfc , psym=-3, pos=five_left_pos, ytitle='sfc shortwave down (w/m^2)', $
	title='Sfc shortwave Down (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_down_shortwave_sfc[where(vector_down_shortwave_sfc gt -1500.)]), $
	10.+max(vector_down_shortwave_sfc[where(vector_down_shortwave_sfc lt 1000.)])]

	oplot, plot_xaxis, vector_down_shortwave_sfc_obs, psym=2, linestyle=2, color=35

endif

;;;;; Surface shortwave up radiation


vector_up_shortwave_sfc=fltarr(n_elements(plot_xaxis))
vector_up_shortwave_sfc_obs=fltarr(n_elements(plot_xaxis))
vector_up_shortwave_sfc_clr=fltarr(n_elements(plot_xaxis))
vector_up_shortwave_sfc_obs_clr=fltarr(n_elements(plot_xaxis))
vector_net_shortwave_sfc_clr=fltarr(n_elements(plot_xaxis))
vector_net_shortwave_sfc_obs_clr=fltarr(n_elements(plot_xaxis))
vector_net_shortwave_sfc=fltarr(n_elements(plot_xaxis))
vector_net_shortwave_sfc_obs=fltarr(n_elements(plot_xaxis))
shortwave_sfc_cld_forcing=fltarr(n_elements(plot_xaxis))
shortwave_sfc_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if solar_flux_up[0,jj] gt 0. and solar_flux_up[0,jj] lt 1500. then begin
	vector_up_shortwave_sfc[j]=(solar_flux_up[0,jj])	; down into the system is positive
		endif else begin
	vector_up_shortwave_sfc[j]=-9999.
		endelse

		if up_short[jj] gt 0. and up_short[jj] lt 1500. then begin
	vector_up_shortwave_sfc_obs[j]=(up_short[jj])	; down into the system is positive
		endif else begin
	vector_up_shortwave_sfc_obs[j]=-9999.
		endelse

		if solar_flux_up_clear[0,jj] gt 0. and solar_flux_up_clear[0,jj] lt 1500. then begin
	vector_up_shortwave_sfc_clr[j]=(solar_flux_up_clear[0,jj])	; down into the system is positive
		endif else begin
	vector_up_shortwave_sfc_clr[j]=-9999.
		endelse

		if up_short[jj] gt 0. and up_short[jj] lt 1500. $
		and vector_down_shortwave_sfc_obs_clr[j] gt 0. $
		and vector_up_shortwave_sfc_obs[j] gt 0. and $
		vector_down_shortwave_sfc_obs[j] gt 0. $
		then begin
	vector_up_shortwave_sfc_obs_clr[j]=vector_down_shortwave_sfc_obs_clr[j]*$
	(vector_up_shortwave_sfc_obs[j]/vector_down_shortwave_sfc_obs[j]) 	; down into the system is positive
		endif else begin
	vector_up_shortwave_sfc_obs_clr[j]=-9999.
		endelse

;net shortwave

	if vector_up_shortwave_sfc[j] gt 0. and vector_down_shortwave_sfc[j] gt 0. then begin
		vector_net_shortwave_sfc[j]=vector_down_shortwave_sfc[j]-vector_up_shortwave_sfc[j]
	endif else begin
		vector_net_shortwave_sfc[j]=-9999.
	endelse

	if vector_up_shortwave_sfc_obs[j] gt 0. and vector_down_shortwave_sfc_obs[j] gt 0. then begin
		vector_net_shortwave_sfc_obs[j]=vector_down_shortwave_sfc_obs[j]-vector_up_shortwave_sfc_obs[j]
	endif else begin
		vector_net_shortwave_sfc_obs[j]=-9999.
	endelse

	if vector_up_shortwave_sfc_clr[j] gt 0. and vector_down_shortwave_sfc_clr[j] gt 0. then begin
		vector_net_shortwave_sfc_clr[j]=vector_down_shortwave_sfc_clr[j]-vector_up_shortwave_sfc_clr[j]
	endif else begin
		vector_net_shortwave_sfc_clr[j]=-9999.
	endelse

	if vector_up_shortwave_sfc_obs_clr[j] gt 0. and vector_down_shortwave_sfc_obs_clr[j] gt 0. then begin
		vector_net_shortwave_sfc_obs_clr[j]=vector_down_shortwave_sfc_obs_clr[j]-vector_up_shortwave_sfc_obs_clr[j]
	endif else begin
		vector_net_shortwave_sfc_obs_clr[j]=-9999.
	endelse

	if vector_net_shortwave_sfc_clr[j] gt 0. and vector_net_shortwave_sfc[j] gt 0. then begin
		shortwave_sfc_cld_forcing[j]=vector_net_shortwave_sfc[j]-vector_net_shortwave_sfc_clr[j]
	endif else begin
		shortwave_sfc_cld_forcing[j]=-9999.
	endelse

	if vector_net_shortwave_sfc_obs_clr[j] gt 0. and vector_net_shortwave_sfc_obs[j] gt 0. then begin
		shortwave_sfc_cld_forcing_obs[j]=vector_net_shortwave_sfc_obs[j]-vector_net_shortwave_sfc_obs_clr[j]
	endif else begin
		shortwave_sfc_cld_forcing_obs[j]=-9999.
	endelse

jj=jj+1
endfor



;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then begin
	plot, plot_xaxis, vector_up_shortwave_sfc , psym=-3, pos=five_right_pos, ytitle='sfc shortwave up (w/m^2)', $
	title='Sfc shortwave up (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_up_shortwave_sfc_obs[where(vector_up_shortwave_sfc_obs gt -1500.)]), $
	10.+max(vector_up_shortwave_sfc[where(vector_up_shortwave_sfc lt 1000.)])]

	oplot, plot_xaxis, vector_up_shortwave_sfc_obs, psym=2, linestyle=2, color=35

endif

; now plot the shortwave cloud forcing

;TOA_solar_up, TOA_solar_down
if window_5 ne 'no' then begin
	plot, plot_xaxis, shortwave_sfc_cld_forcing , psym=-3, pos=six_left_pos, ytitle='sfc shortwave Cloud Forcing (w/m^2)', $
	title='Sfc shortwave Cloud Forcing (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(shortwave_sfc_cld_forcing[where(shortwave_sfc_cld_forcing gt -1500.)]), $
	10.+max(shortwave_sfc_cld_forcing[where(shortwave_sfc_cld_forcing lt 1000.)])]

	oplot, plot_xaxis, shortwave_sfc_cld_forcing_obs, psym=2, linestyle=2, color=35

endif



if window_5 ne 'no' then xyouts, 0.5,0.97, title_string, /normal, alignment=0.5, color=0

if window_5 ne 'no' then write_png, fname_g_prefix+'plot5'+fname_g_suffix, tvrd(), red, green, blue
if window_5 ne 'no' then spawn, 'convert '+fname_g_prefix+'plot5'+fname_g_suffix+' '+fname_g_prefix+'plot5.gif'


endif	; if window_5 ne 'no' then begin
if window_6 ne 'no' then begin
if window_6 ne 'no' then WINDOW, 6, XSIZE=1000, YSIZE=1000, TITLE='ARM Integrated Column Product Plots', /pixmap

!p.multi = [0,6,6]
;height=height/1000.
!p.charsize=2.0

if n_elements((where(vector_scf_toa lt 1500. and vector_scf_toa gt -1500.))) eq 1 then window_6='no'

if window_6 ne 'no' then begin
	plot, plot_xaxis, vector_scf_toa , psym=-3, pos=one_left_pos, ytitle='Shortwave CRF (w/m^2)', $
	title='TOA Shortwave Cloud Radiative Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_scf_toa[where(vector_scf_toa gt -1500.)]), $
	10.+max(vector_scf_toa[where(vector_scf_toa lt 100.)])]

	if n_elements(hrfrac_M) gt 1 then oplot, hrfrac_M, vector_scf_toa_M, psym=2, linestyle=2, color=35

endif

if window_6 ne 'no' then begin

	plot, plot_xaxis, vector_lcf_toa , psym=-3, pos=one_right_pos, ytitle='Longwave CRF (w/m^2)', $
	title='TOA Longwave Cloud Radiative Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_lcf_toa[where(vector_lcf_toa gt -1500.)]), $
	10.+max(vector_lcf_toa[where(vector_lcf_toa lt 100.)])]

	if n_elements(hrfrac_M) gt 1 then oplot, hrfrac_M, vector_lcf_toa_M, psym=2, linestyle=2, color=35

endif

if window_6 ne 'no' then begin
	plot, plot_xaxis, vector_netcf_toa , psym=-3, pos=two_left_pos, ytitle='Net CRF (w/m^2)', $
	title='TOA Net Cloud Radiative Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(vector_netcf_toa[where(vector_netcf_toa gt -1500.)]), $
	10.+max(vector_netcf_toa[where(vector_netcf_toa lt 100.)])]

	if n_elements(hrfrac_M) gt 1 then oplot, hrfrac_M, vector_netcf_toa_M, psym=2, linestyle=2, color=35

endif


if window_6 ne 'no' then begin
	plot, plot_xaxis, longwave_sfc_cld_forcing , psym=-3, pos=two_right_pos, ytitle='sfc longwave cloud forcing (w/m^2)', $
	title='Sfc Longwave cld Forcing (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(longwave_sfc_cld_forcing[where(longwave_sfc_cld_forcing gt -1500.)]), $
	10.+max(longwave_sfc_cld_forcing[where(longwave_sfc_cld_forcing lt 1000.)])]

	oplot, plot_xaxis, longwave_sfc_cld_forcing_obs, psym=-3, linestyle=2, color=35

endif

if window_6 ne 'no' then begin
	plot, plot_xaxis, shortwave_sfc_cld_forcing , psym=-3, pos=three_left_pos, ytitle='sfc shortwave Cloud Forcing (w/m^2)', $
	title='Sfc shortwave Cloud Forcing (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(shortwave_sfc_cld_forcing[where(shortwave_sfc_cld_forcing gt -1500.)]), $
	10.+max(shortwave_sfc_cld_forcing[where(shortwave_sfc_cld_forcing lt 1000.)])]

	oplot, plot_xaxis, shortwave_sfc_cld_forcing_obs, psym=-3, linestyle=2, color=35

endif

net_sfc_cld_forcing=fltarr(n_elements(plot_xaxis))
net_sfc_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if shortwave_sfc_cld_forcing[j] gt -9999. and longwave_sfc_cld_forcing[j] gt -9999. then begin
	net_sfc_cld_forcing[j]=shortwave_sfc_cld_forcing[j]+longwave_sfc_cld_forcing[j]	; down into the system is positive
		endif else begin
		if longwave_sfc_cld_forcing[j] gt -9999. then begin
	net_sfc_cld_forcing[j]=longwave_sfc_cld_forcing[j]
		endif else begin
	net_sfc_cld_forcing[j]=-9999.
		endelse
		endelse

		if shortwave_sfc_cld_forcing_obs[j] gt -9999. and longwave_sfc_cld_forcing_obs[j] gt -9999. then begin
	net_sfc_cld_forcing_obs[j]=shortwave_sfc_cld_forcing_obs[j]+longwave_sfc_cld_forcing_obs[j]	; down into the system is positive
		endif else begin
	if longwave_sfc_cld_forcing_obs[j] gt -9999. then begin
	net_sfc_cld_forcing_obs[j]=longwave_sfc_cld_forcing_obs[j]
		endif else begin
	net_sfc_cld_forcing_obs[j]=-9999.
		endelse
		endelse

jj=jj+1
endfor


if window_6 ne 'no' then begin
	plot, plot_xaxis, net_sfc_cld_forcing , psym=-3, pos=three_right_pos, ytitle='sfc net Cloud Forcing (w/m^2)', $
	title='Net Sfc Cloud Forcing (Obs-dashed)',background=!d.n_colors-1, $
	color=0, yrange=[min(net_sfc_cld_forcing[where(net_sfc_cld_forcing gt -1500.)]), $
	10.+max(net_sfc_cld_forcing[where(net_sfc_cld_forcing lt 1000.)])]

	oplot, plot_xaxis, net_sfc_cld_forcing_obs, psym=-3, linestyle=2, color=35

endif


shortwave_atm_cld_forcing=fltarr(n_elements(plot_xaxis))
shortwave_atm_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
longwave_atm_cld_forcing=fltarr(n_elements(plot_xaxis))
longwave_atm_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
net_atm_cld_forcing=fltarr(n_elements(plot_xaxis))
net_atm_cld_forcing_obs=fltarr(n_elements(plot_xaxis))
jj=start_time_index
for j=0,n_elements(plot_xaxis)-1 do begin

		if shortwave_sfc_cld_forcing[j] gt -9999. and vector_scf_toa[j] gt -9999. then begin
	shortwave_atm_cld_forcing[j]=vector_scf_toa[j]-shortwave_sfc_cld_forcing[j]	; down into the system is positive
		endif else begin
	shortwave_atm_cld_forcing[j]=-9999.
		endelse

		if longwave_sfc_cld_forcing[j] gt -9999. and vector_lcf_toa[j] gt -9999. then begin
	longwave_atm_cld_forcing[j]=vector_lcf_toa[j]-longwave_sfc_cld_forcing[j]	; down into the system is positive
		endif else begin
	longwave_atm_cld_forcing[j]=-9999.
		endelse

		if shortwave_sfc_cld_forcing[j] gt -9999. and longwave_sfc_cld_forcing[j] gt -9999. then begin
	net_atm_cld_forcing[j]=shortwave_atm_cld_forcing[j]-longwave_atm_cld_forcing[j]	; down into the system is positive
		endif else begin
		if shortwave_sfc_cld_forcing[j] le 0. and longwave_atm_cld_forcing[j] gt -9999. then begin
			net_atm_cld_forcing[j]=longwave_atm_cld_forcing[j]
		endif else begin
	net_atm_cld_forcing[j]=-9999.
		endelse
		endelse

if n_elements(hrfrac_M) gt 1 then begin 

time_index=where(abs(hrfrac_M-hrfrac[jj]) eq min(abs(hrfrac_M-hrfrac[jj])))

		if shortwave_sfc_cld_forcing_obs[j] gt -9999. and vector_scf_toa_M[time_index[0]] gt -9999. and $
		abs(hrfrac_M[time_index[0]]-hrfrac[jj]) le 5.d/1440.d then begin
	shortwave_atm_cld_forcing_obs[j]=vector_scf_toa_M[time_index[0]]-shortwave_sfc_cld_forcing_obs[j]	; down into the system is positive
		endif else begin
	shortwave_atm_cld_forcing_obs[j]=-9999.
		endelse

		if longwave_sfc_cld_forcing_obs[j] gt -9999. and vector_lcf_toa_M[time_index[0]] gt -9999. and $
		abs(hrfrac_M[time_index[0]]-hrfrac[jj]) le 5.d/1440.d then begin
	longwave_atm_cld_forcing_obs[j]=vector_lcf_toa_M[time_index[0]]-longwave_sfc_cld_forcing_obs[j]	; down into the system is positive
		endif else begin
	longwave_atm_cld_forcing_obs[j]=-9999.
		endelse

;vector_netcf_toa_M[kk]=vector_lcf_toa_M[kk]+vector_scf_toa_M[kk]

		if shortwave_sfc_cld_forcing_obs[j] gt -9999. and longwave_atm_cld_forcing_obs[j] gt -9999. then begin
	net_atm_cld_forcing_obs[j]=shortwave_atm_cld_forcing_obs[j]-longwave_atm_cld_forcing_obs[j]	; down into the system is positive
		endif else begin
		if shortwave_sfc_cld_forcing_obs[j] le 0. and longwave_atm_cld_forcing_obs[j] gt -9999. then begin
			net_atm_cld_forcing_obs[j]=longwave_atm_cld_forcing_obs[j]
		endif else begin
	net_atm_cld_forcing_obs[j]=-9999.
		endelse
		endelse

endif		; if n_elements(hrfrac_M) gt 1 then begin

jj=jj+1
endfor

if window_6 ne 'no' then begin
	plot, plot_xaxis, shortwave_atm_cld_forcing , psym=-3, pos=four_left_pos, ytitle='Atm Shortwave Cloud Forcing (w/m^2)', $
	title='Shortwave Atmospheric Cloud Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(shortwave_atm_cld_forcing[where(shortwave_atm_cld_forcing gt -1500.)]), $
	10.+max(shortwave_atm_cld_forcing[where(shortwave_atm_cld_forcing lt 1000.)])]

	if n_elements(hrfrac_M) gt 1 then oplot, plot_xaxis, shortwave_atm_cld_forcing_obs, psym=2, linestyle=2, color=35

	plot, plot_xaxis, longwave_atm_cld_forcing , psym=-3, pos=four_right_pos, ytitle='Atm Longwave Cloud Forcing (w/m^2)', $
	title='Longwave Atmospheric Cloud Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(longwave_atm_cld_forcing_obs[where(longwave_atm_cld_forcing_obs gt -1500.)]), $
	10.+max(longwave_atm_cld_forcing[where(longwave_atm_cld_forcing lt 1000.)])]

	if n_elements(hrfrac_M) gt 1 then oplot, plot_xaxis, longwave_atm_cld_forcing_obs, psym=2, linestyle=2, color=35

	plot, plot_xaxis, net_atm_cld_forcing , psym=-3, pos=five_left_pos, ytitle='Atm net Cloud Forcing (w/m^2)', $
	title='Net Atmospheric Cloud Forcing (GOES -> *)',background=!d.n_colors-1, $
	color=0, yrange=[min(net_atm_cld_forcing[where(net_atm_cld_forcing gt -1500.)]), $
	10.+max(net_atm_cld_forcing[where(net_atm_cld_forcing lt 1000.)])]

	if n_elements(hrfrac_M) gt 1 then oplot, plot_xaxis, net_atm_cld_forcing_obs, psym=2, linestyle=2, color=35

endif


if window_6 ne 'no' then xyouts, 0.5,0.97, title_string, /normal, alignment=0.5, color=0

if window_6 ne 'no' then write_png, fname_g_prefix+'plot6'+fname_g_suffix, tvrd(), red, green, blue
if window_6 ne 'no' then spawn, 'convert '+fname_g_prefix+'plot6'+fname_g_suffix+' '+fname_g_prefix+'plot6.gif'


endif ;if window_6 ne 'no' then begin
; plot the TOA IR and solar forcing

;endif	; window_5

endif ; if data_flag eq 1 then begin


endif	; fnf
;spawn, 'gzip '+fname+' &'
end