;******************************************************************************;

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/popRemap.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

;***************************************************
; READ IN DATA
;***************************************************

in1 = addfile("siconc_SIday_BCC-CSM2-MR_historical_r1i1p1f1_gn_19000101-19091231.nc","r")
siconc = in1->siconc
lat_bcc = in1->latitude
lon_bcc = in1->longitude

dims=dimsizes(siconc)

j_bcc=dims(1)
i_bcc=dims(2)

;************************************
; AVERAGE OVER SUBSET OF DATA IGNORING ZEROS
;************************************

; Subset data 
startDay=0
endDay=364
numDays=365

subset = siconc(startDay:endDay,:,:)

; Sum values over time 
sum_subset=dim_sum_n_Wrap(subset,0)

; Find # of days where the subset is non-zero for each point 
zeros_subset=where(subset.eq.0.0,0.0,1)
sum_zeros_subset=dim_sum_n_Wrap(zeros_subset,0)
sum_zeros_subset=where(sum_zeros_subset.eq.0.0,1,sum_zeros_subset)

; Find average of non-zero values  
avg_subset=sum_subset/sum_zeros_subset

; Set lat and lon 
avg_subset@lon2d = lon_bcc
avg_subset@lat2d = lat_bcc

;***************************************************************************
; SET PLOT FEATURES           
;***************************************************************************
 wks = gsn_open_wks("pdf","siconc_test")   

 plot1 = new(1,graphic)    

	res                       = True
	res@gsnPolar              = "NH"               
    res@gsnAddCyclic		  = True

	res@cnFillOn              = True
	res@cnLinesOn			  = False
	res@cnMonoLevelFlag       = False
	res@cnLevelFlags          = "NoLine" 
    res@cnLineLabelInterval   = 0
    res@cnLevelSelectionMode  = "ExplicitLevels"   	
	res@cnLevels              = (/1,10,20,30,40,50,60,70,80,90,100/)    
	
    res@mpOutlineOn           = True          
	res@mpLandFillColor       = "white"           
	res@mpMinLatF             = 57                 	
	res@mpGridLonSpacingF     = 90.                 

	res@cnMissingValPerimOn     = False      
	res@cnMissingValFillColor   = "white"
	res@cnMissingValFillPattern =  0     
	res@cnFillPalette 		    =   "CBR_wet"

                                    
	res@lbLabelBarOn = True
	res@lbTitleString    		= "ice concentration (%)"
	res@gsnDraw                 = True   
	res@gsnFrame                 = True          
 

;************************************
; CREATE PLOT
;************************************
 
plot1(0)=gsn_csm_contour_map_polar(wks,avg_subset,res)

end