;   - Changing the line thickness for multiple curves in an XY plot
;   - Drawing XY plot curves with both lines and markers
;   - Changing the default markers in an XY plot
;   - Making all curves in an XY plot solid

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

  wks_type = "pdf" ; figure format
  
;+++++++++++++++++++ Buoy data +++++++++++++++++++++++++
; read data
  lines = asciiread("buoywind_radii_all.dat",-1,"string")
  year  = stringtointeger(str_get_field(lines(0:),1," "))
  month = stringtointeger(str_get_field(lines(0:),2," "))
  day   = stringtointeger(str_get_field(lines(0:),3," "))
  hour  = stringtointeger(str_get_field(lines(0:),4," "))
  dist  = stringtofloat(str_get_field(lines(0:),6," "))   ; distance to TC center
  az    = stringtofloat(str_get_field(lines(0:),7," "))   ; azimuth angle to TC center
  wsd    = stringtofloat(str_get_field(lines(0:),8," "))   ; wind speed
  delete(lines)

  az = where(az.lt.0, az+360, az) ; convert azimuth angle from negative to positive
  hr = (day-day(0))*24+hour       ; compute hours from initial time

;+++++++++++++++++++ Land data +++++++++++++++++++++++++
; read data
  lines = asciiread("sfcwind_radii_all.dat",-1,"string")
  yearx = stringtointeger(str_get_field(lines(0:),1," "))
  monthx = stringtointeger(str_get_field(lines(0:),2," "))
  dayx   = stringtointeger(str_get_field(lines(0:),3," "))
  hourx  = stringtointeger(str_get_field(lines(0:),4," "))
  distx  = stringtofloat(str_get_field(lines(0:),5," "))   ; distance to TC center
  azx    = stringtofloat(str_get_field(lines(0:),6," "))   ; azimuth angle to TC center
  wsdx   = stringtofloat(str_get_field(lines(0:),7," "))   ; wind speed
  delete(lines)

  azx = where(azx.lt.0, azx+360, azx) ; convert azimuth angle from negative to positive
  hrx = (dayx-dayx(0))*24+hourx       ; compute hours from initial time

; data leading time attached to plot file name
  nt=dimsizes(year)
  ntx=dimsizes(yearx)
  ys=sprinti("%0.4i", year(nt-1))
  ms=sprinti("%0.2i", month(nt-1))
  dd=sprinti("%0.2i", day(nt-1))
  hs=sprinti("%0.2i", hour(nt-1))
  hsx=sprinti("%0.2i", hourx(ntx-1))

  if(hs.ge.hsx)
  leadt=ys+ms+dd+hs
  else
  leadt=ys+ms+dd+hsx
  end if
  print(leadt+"")

; define plot configuration
  res               = True
  res@gsnDraw     = False
  res@gsnFrame    = False
  res@gsnMaximize   = True
  res@tiYAxisString = "Radii (nm)" 
  res@tiXAxisString = "Time start ("+month(0)+"/"+day(0)+"/"+year(0)+")"
  res@tmYROn = False
  res@tmXTOn = False
  res@xyMarkLineMode = "Markers"
  res@xyMarkers      = 16
  res@xyMarkerColor  = "NavyBlue"
  res@xyMarkerSizeF  = 0.01

; Set up resources for the labelbar.
  lbres                    = True          ; labelbar only resources
  ;lbres@lbAutoManage       = True          ; Necessary to control sizes
  lbres@lbPerimOn          = True
  lbres@lgMonoDashIndex    = True
  lbres@lbMonoFillPattern  = True          ; Solid fill pattern
  ;lbres@lgLineColors       = (/"red","NavyBlue"/)
  ;lbres@lgLineThicknessF   = 2.
  lbres@lbPerimColor       = "green"
  lbres@lbFillColors       = (/"Red","NavyBlue"/) ; labelbar colors
  labels                   = (/"Land","Buoy"/)
  lbres@lbBoxMajorExtentF  = 0.1          ; puts space between color boxes
   lbres@lbFillDotSizeF = 0.1
  lbres@lbLabelFontHeightF = 0.01          ; font height. default is small
  lbres@vpWidthF           = 0.15
  lbres@vpHeightF          = 0.1 
  lbres@lbLabelJust        = "CenterLeft"  ; left justify labels

  amres                  = True
  amres@amJust           = "BottomRight"  ;"TopRight"
  amres@amParallelPosF   = 0.5
  amres@amOrthogonalPosF = 0.5

; set up land scatter plot
  pmres                  = True
  pmres@gsMarkerIndex    = 16
  pmres@gsMarkerSizeF    = 0.01
  pmres@gsMarkerColor    = "Red"
  pmres@gsLineThicknessF = 2.0

; wind speed window
  ws  = (/23,27,38,42,62,66/)
  wss = (/"25","40","64"/)

; quadrants
  qd  = (/0,90,180,270,360/)
  qdr = (/"NE","SE","SW","NW"/)

  ds0=dist
  ds0x=distx

  do i = 0, 2   ; loop 3 wind windows

; exclude wind speed outside window
  i1=(i+1)*2-2
  i2=(i+1)*2-1
  ds1  = where(wsd.lt.ws(i1) .or. wsd.gt.ws(i2), -999 , ds0)
  ds1x = where(wsdx.lt.ws(i1) .or. wsdx.gt.ws(i2), -999 , ds0x)

  do j = 0, 3    ; loop 4 quadrants

; exclude other quadrants
  ds2  = where(az.lt.qd(j) .or. az.gt.qd(j+1), -999 , ds1)
  ds2x = where(azx.lt.qd(j) .or. azx.gt.qd(j+1), -999 , ds1x)

; make radii plot within this wind window & quadrant
  if(max(ds2).gt.0.and.max(ds2x).gt.0)

  res@tiMainString  = "Radii "+wss(i)+"kts "+qdr(j)
  wks = gsn_open_wks(wks_type,"LandBuoy"+leadt+"_R"+wss(i)+"_"+qdr(j))

; ++++++ buoy wind radii ++++++++++++++++++++
  ds2@_FillValue = -999
  ds2x@_FillValue = -999

  if(min(ds2).le.min(ds2x))
  res@trYMinF = min(ds2)        ; set minimum Y-axis value
  else
  res@trYMinF = min(ds2x)        ; set minimum Y-axis value
  end if

  if(max(ds2).ge.max(ds2x))
  res@trYMaxF = max(ds2)        ; set maximum Y-axis value
  else
  res@trYMaxF = max(ds2x)        ; set maximum Y-axis value
  end if

  plot = gsn_xy(wks,hr,ds2,res) ; create plot
  delete(ds2)

; ++++++ Add land surface wind radii ++++++++++++++++++++
  dum = gsn_add_polymarker(wks, plot, hrx, ds2x, pmres)
  delete(ds2x)

; create labelbar
  lbid = gsn_create_labelbar(wks,2,labels,lbres)

; Insert labelbar on plot
  annoid = gsn_add_annotation(plot,lbid,amres)

  draw(plot)
  frame(wks)
  delete(plot)
  delete(wks)
  delete(dum)

  end if

  end do  ; loop 4 quadrants

  end do  ; loop 3 wind windows

  end