;concepts illustrated:
;   - Drawing contours over a polar stereographic map
;   - Creating a color map using named colors
;   - Adding a common labelbar to paneled plots
;   - Paneling four plots on a page
;   - Changing the minimum latitude for a polar stereographic map
;   - Subselecting an array using its coordinate arrays
;before running
;>csh
;setenv NCARG_COLORMAPS /home/evi/data/data0_5/nc_data/colormaps:$NCARG_ROOT/lib/ncarg/colormaps
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;************************************************
begin
;************************************************
; read in netCDF file
;************************************************
  filename_a00="T20081029_12.nc"
  gribfile_a00=addfile(filename_a00,"r")

;*****************************************
;plotting GRADPV
;*********************************
 pv1_00 = gribfile_a00->GRADPV(:,:)
 lat= gribfile_a00->lat
 lon= gribfile_a00->lon

;************************************************
; create plot
;************************************************
  wks = gsn_open_wks("png","GRADPV")               ; open a ps file
;***********************************************
; create array of named colors. Note that you have
; to include the forground (white), and background (black)
; as the first two colors
;***********************************************
data=pv1_00

data=0
nx=dimsizes(lat)-1
ny=dimsizes(lon)-1

 do ii=0,nx
  do jj=0,ny
if (pv1_00(ii,jj) .gt. 15)then
data(ii,jj)=pv1_00(ii,jj)
else
data(ii,jj)=0
end if
end do
end do

;########################
;#colors for gradpv
;#####################

 colors=(/"white","black","white","royalblue1","cyan","aquamarine"\
           ,"palegreen","olivedrab1","gold","yellow","orange"\
         , "chocolate1","brown1","brown4"/)



   gsn_define_colormap(wks, colors)


;*******************************************************************
;320 isentropic surface hour 00
;*************************************************************
  res = True                                     ; plot mods desired
                         
;res@gsnCenterString      = "GRADPV 320K isentropic surface" 

 
 res@gsnLeftString        = ("2008/10/29_12")
  res@gsnRightString       = "" 
  
  res@cnFillOn             = True               ; turns on the color
  res@mpFillOn             = False              ; turns off continent gray
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnLevelSelectionMode = "ExplicitLevels"  ; AutomaticLevels,ManualLevels,ExplicitLevels,EqualSpacedLevels
  res@cnLevels        = (/5,10,13,15,19,21,24,27,30,36/)  
  res@mpMinLatF                 =  20;10.         ; range to zoom in on
  res@mpMaxLatF                 = 70 ;80.
  res@mpMinLonF                 = -70 ;-80.
  res@mpMaxLonF                 = 50 ;60.
  res@gsnMaximize      = True
  res@gsnAddCyclic =False
  res@lbLabelAutoStride    =   True

;&&&&&&&&&&&&&&&&&&&PLOT&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


    plot = gsn_csm_contour_map(wks,data(:,:),res)

gsn_define_colormap(wks, colors)

;rame (wks)
; draw(plot)

end