;*************************************************
; station_1.ncl
;*************************************************
;
; Concepts illustrated:
;   - Contouring 1D arrays of X,Y,Z station data
;   - Reading an ASCII file with several columns of data
;   - Drawing filled dots on a map
;   - Drawing filled contours over a Lambert Conformal map
;   - Turning off map fill
;   - Controlling which contour lines get drawn
;   - Zooming in on a particular area on a Lambert Conformal map
;   - Reversing a color map
;
; This example reads in station data represented by
; 1D arrays, and generates a filled contour plot
; via the triangular mesh capability in NCL.
;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
;
; Data is stored in four columns:  station_name  lat  lon  pwv
; Read in each line as a string, and use "str_get_field" to 
; read in the fields of interest.
;
  fname = "pw.dat"
  lines = asciiread(fname,-1,"string")

;
; Use "str_get_field" to indicate which fields to read in. Each field
; is separated by spaces.
;
  lat = stringtofloat(str_get_field(lines(1:),2," "))
  lon = stringtofloat(str_get_field(lines(1:),3," "))
  pwv = stringtofloat(str_get_field(lines(1:),4," "))

;
; This second file is not so tricky.  The 2D lat/lon data is sorted with 
; lat values first, and then lon values.
;
  nlat     = 70
  nlon     = 70
  latlon2d = asciiread("stn_latlon.dat",(/2,nlat,nlon/),"float")
  lat2d    = latlon2d(0,:,:)
  lon2d    = latlon2d(1,:,:)

  delete(lines)        ; Remove arrays we don't need anymore.
  delete(latlon2d)

  tlat1 = 60.0
  tlat2 = 30.0
  clon  = -98.5
  clat  =  36.3

  wks = gsn_open_wks("x11","station")
  gsn_define_colormap(wks,"WhViBlGrYeOrRe")
  gsn_reverse_colormap(wks)

  res                             = True
  res@gsnFrame                    = False   ; So we can draw markers
  res@gsnMaximize                 = True

  res@gsnSpreadColors             = True
  res@gsnSpreadColorEnd           = 90
  res@pmTickMarkDisplayMode       = "Always"

  res@trGridType                  = "TriangularMesh"    ; The default if you
                                                        ; have 1D data

  res@cnLineLabelPlacementMode    = "Constant"
  res@cnLineDashSegLenF           = 0.3

  res@cnLevelSelectionMode        = "ManualLevels"
  res@cnMinLevelValF              = 15         	; 15.25
  res@cnMaxLevelValF              = 50		; 49.75
  res@cnLevelSpacingF             = 0.25

  res@cnFillOn                    = True
  res@cnLinesOn                   = True
  res@cnLineLabelsOn              = True
  res@cnLevelFlags		  = new(139,"string")
  res@cnLevelFlags(:)		  = "NoLine"
  res@cnLevelFlags(0::20)         = "LineAndLabel"

  res@lbLabelStride               = 20  
  res@lbBoxLinesOn                = False
  res@tiMainString                = "GPS PWV (18Z)"

  res@sfXArray                    = lon
  res@sfYArray                    = lat

  res@mpProjection                = "LambertConformal"
  res@mpLambertParallel1F         = tlat1
  res@mpLambertParallel2F         = tlat2
  res@mpLambertMeridianF          = clon

  res@mpLimitMode                 = "Corners"
  res@mpLeftCornerLatF            = lat2d(20,18)
  res@mpLeftCornerLonF            = lon2d(20,18)
  res@mpRightCornerLatF           = lat2d(58,52)
  res@mpRightCornerLonF           = lon2d(58,52)

  res@mpFillOn                    = False
  res@mpOutlineDrawOrder          = "PostDraw"
  res@mpFillDrawOrder             = "PreDraw"
  res@mpOutlineBoundarySets       = "GeophysicalAndUSStates"
  res@mpUSStateLineColor          = "Gray10"
  res@mpUSStateLineDashPattern    = 2

;
; Draw the map (frame won't get advanced because gsnFrame was set to False).
;
  map = gsn_csm_contour_map(wks,pwv,res)

;
; Draw markers on the plot in the lat/lon locations.
;
  mkres               = True
  mkres@gsMarkerIndex = 17     ; Filled circle
  mkres@gsMarkerSizeF = 0.03

  gsn_polymarker(wks,map,lon,lat,mkres)

  frame(wks)   ; Now advance the frame.

end