;********************************************************************
; rcm_4.ncl
;
; Concepts illustrated:
;   - Plotting RCM precipitation data
;   - Drawing line contours over a lambert conformal map
;   - Using a filter to convert Cray COS block binary data to netCDF 
;********************************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; 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"
;
; This file still has to be loaded manually
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
;*************************
; read in precip data 
; (units cm)
;*************************
  a = addfile("prc.1982112500-1983030500.nc","r") 
  x = a->prc                                                ; read in data
  dims=dimsizes(x)                                          ; get dim sizes
  nlat=dims(1)                            
  nlon=dims(2)
  
  LAT2D = a->lat2d
  LON2D = a->lon2d
;*************************
; set some parameters and
; read in date
;************************* 
  date_start    = 1982120100    ; pick start and end dates
  date_end      = 1983010100
  data_per_day  = 4
  date = a->date                        ; read in date                         
  index_start = ind(date.eq.date_start) ; get array index of desired times
  index_end   = ind(date.eq.date_end  )
  
  nday = (index_end-index_start)/data_per_day
;*************************
; compute monthly average
; units: mm/day               
;*************************
  y = new((/nlat,nlon/),float)                            ; allocate memory
  y = (/x(index_end,:,:) - x(index_start,:,:)/)*10./nday  ; compute
  y!0 ="lat"                                              ; name dimensions
  y!1 ="lon"
;*************************
; prepare data for ploting
; pull out all data except 
; last value in x and y 
; (bad points)
;************************
  precip = y(lat|0:nlat-2,lon|0:nlon-2)             
  precip@long_name="Precipitation"   ; give data long_name and units,
  precip@units    ="mm/day"          ; so plotted automatically. 
;***********************
; plot
;***********************
  wks = gsn_open_wks("png","rcm")          ; send graphics to PNG file
  
  res            = True                    ; plot mods desired
  
  res@gsnMaximize       = True

  res@mpLimitMode       = "Corners"        ; choose range of map
  res@mpLeftCornerLatF  = LAT2D(0,0)
  res@mpLeftCornerLonF  = LON2D(0,0)
  res@mpRightCornerLatF = LAT2D(nlat-2,nlon-2)
  res@mpRightCornerLonF = LON2D(nlat-2,nlon-2)
  
; The following 4 pieces of information are REQUIRED to properly display
; data on a native lambert conformal grid.  This data should be specified
; somewhere in the model itself.
  
; WARNING: our local RCM users could not provide us with this information,
; so this is our best guess as to the correct numbers. Use at your own risk.
  
  res@mpProjection        = "LambertConformal"
  res@mpLambertParallel1F = 30
  res@mpLambertParallel2F = 58
  res@mpLambertMeridianF  = 260
  
;
; Usually, when data is placed onto a map, it is transformed to the specified
; projection. Since this model is already on a native lambert conformal grid,
; we want to turn OFF the transformation.
;  
  res@tfDoNDCOverlay        = True
;  res@tfDoNDCOverlay        = "NDCViewport"  ; NCL V6.5.0 or later
  
  res@cnLineLabelFontHeightF     = .013    ; label font height
  res@cnLineThicknessF           = 1.5     ; make contour lines thicker
  
  res@cnHighLabelsOn             = True           ; high labels on
  res@cnHighLabelFontHeightF     = .013   
  res@cnLowLabelsOn              = True           ; low label on
  res@cnLowLabelFontHeightF      = .013   

  res@mpPerimOn              = True        ; draw box around map
  res@mpGridLineDashPattern  = 2           ; lat/lon lines as dashed
  res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries

  res@tiMainString      = "RCM2-82nmc03,Starting 19821121" ; title
                                     
  plot = gsn_contour_map(wks,precip,res)     ; Draw contours over a map.

  opt  = True                                ; set shading resources
  opt@gsnShadeFillType = "pattern"           ; select pattern fill
  opt@gsnShadeHigh     = 17                  ; select stipple pattern for fill

  plot = gsn_contour_shade(plot,-18,3.99,opt)        ; shade all contours > 4 contour

end