台风路径图

YJane

           想用NCL做台风路径图，用自己的数据资料编辑，但是提示说是    “。。除以0，无法进行”
我读不懂这是具体指哪里出错，望大神能指点。

Copyright (C) 1995-2014 - All Rights Reserved
University Corporation for Atmospheric Research
NCAR Command Language Version 6.2.1
The use of this software is governed by a License Agreement.
See http://www.ncl.ucar.edu/ for more details.
(0)    Input is ARW type data
(0)    Doing a single time plot
(0)    Data interval is 0 hour(s)
(0)    The Track Plot Window is set to
(0)       Latitute  from 8 to 40
(0)       Longitute from -90 to 139.5
(0)    The buffer between the data and the plot edge is 2 degrees
fatal:divide: Division by 0, Can't continue
fatal:Div: operator failed, can't continue
fatal:["Execute.c":8578]:Execute: Error occurred at or near line 716 in file track

fatal:["Execute.c":8578]:Execute: Error occurred at or near line 1273 in file track

下面是我的输入文件：
ATCF 2003-07-18_06:00:00    10.0   137.5  995   20
ATCF 2003-07-18_12:00:00    10.2   137.1  990   23
ATCF 2003-07-18_18:00:00    10.6   136.5  990   25
ATCF 2003-07-19_00:00:00    10.6   135.6  985   28
ATCF 2003-07-19_06:00:00    10.5   134.8  980   30
ATCF 2003-07-19_12:00:00    10.6   134.2  980   30
ATCF 2003-07-19_18:00:00    10.9   133.4  975   33
ATCF 2003-07-20_00:00:00    11.5   132.7  970   35
ATCF 2003-07-20_06:00:00    12.1   131.6  960   40
ATCF 2003-07-20_12:00:00    12.5   130.7  950   45
ATCF 2003-07-20_18:00:00    13.3   129.6  950   45
ATCF 2003-07-21_00:00:00    13.6   128.2  950   45
ATCF 2003-07-21_06:00:00    14.0   127.1  945   50
ATCF 2003-07-21_12:00:00    15.0   125.9  945   50
ATCF 2003-07-21_18:00:00    15.8   124.5  945   50
ATCF 2003-07-22_00:00:00    16.4   123.0  950   45
ATCF 2003-07-22_06:00:00    16.9   121.4  950   45
ATCF 2003-07-22_12:00:00    17.7   119.6  950   45
ATCF 2003-07-22_18:00:00    18.0   118.4  950   45
ATCF 2003-07-23_00:00:00    18.4   117.0  955   40
ATCF 2003-07-23_06:00:00    18.8   115.5  955   40
ATCF 2003-07-23_12:00:00    19.4   114.2  955   40
ATCF 2003-07-23_18:00:00    20.0   112.9  950   45
ATCF 2003-07-24_00:00:00    21.1   111.9  950   45
ATCF 2003-07-24_06:00:00    22.1   110.1  975   30
ATCF 2003-07-24_12:00:00    22.8   108.9  985   23
ATCF 2003-07-24_18:00:00    23.4   107.3  995   20
ATCF 2003-07-25_00:00:00    23.5   106.0  996   15
ATCF 2003-07-25_06:00:00    23.1   105.0  996   15

设置的NCL 文件：

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


;_______________________________________________________________________________________
; White and Black are background colors
; LightBlue and SandyBrown are used for the map background
; Azure4 (darkgrey) is used for the map borders and frame
; Blue, DarkGreen, Yellow, Red and Purple3 are used for the hurricane symbols (5,6,7,8,9)
; Add colors below "Purple3" if more are needed

undef("TrackColors")
function TrackColors()
begin
  track_colors = (/"White","Black", \
                   "Azure4", \
                   "LightBlue", \
                   "SandyBrown", \
                   "Blue", \
                   "DarkGreen", \
                   "Yellow", \
                   "Red", \
                   "Purple3", \
                   "CadetBlue4", \
                   "Olivedrab1", \
                   "Green", \
                   "Orchid1", \
                   "Cyan", \
                   "CornFlowerBlue"/)
  return(track_colors)
end

;_______________________________________________________________________________________
; Make sure the final window is bigger than the data spread

undef("MapWindow")
procedure MapWindow(lat,lon,wnd,tkres)
begin

  do i=0,dimsizes(wnd)-1
    if ( wnd(i) .eq. -1.0 ) then
      lat(i) = 30.
      lon(i) = -80.
    end if
  end do
  min_lat = min(lat)
  max_lat = max(lat)
  min_lon = min(lon)
  max_lon = max(lon)

  if ( min_lat .lt. tkres@latMIN .and. .not. tkres@latMINset ) then
    min_lat = min_lat - tkres@disEDGE
  else
    min_lat = tkres@latMIN
  end if
  tkres@min_lat = min_lat

  if ( max_lat .gt. tkres@latMAX .and. .not. tkres@latMAXset ) then
    max_lat = max_lat + tkres@disEDGE
  else
    max_lat = tkres@latMAX
  end if
  tkres@max_lat = max_lat

  if ( min_lon .lt. tkres@lonMIN .and. .not. tkres@lonMINset ) then
    min_lon = min_lon - tkres@disEDGE
  else
    min_lon = tkres@lonMIN
  end if
  tkres@min_lon = min_lon

  if ( max_lon .gt. tkres@lonMAX .and. .not. tkres@lonMAXset ) then
    max_lon = max_lon + tkres@disEDGE
  else
    max_lon = tkres@lonMAX
  end if
  tkres@max_lon = max_lon

end

;_______________________________________________________________________________________
; Set up the basic map resources                             

undef("MapResources")
procedure MapResources(res,tkres)
begin

  res             = True

  res@gsnMaximize = True     ; Maximize plot size in the frame.
  res@gsnDraw     = False    ; Don't draw the plot just yet.
  res@gsnFrame    = False    ; Don't advance the frame just yet.

  res@pmTickMarkDisplayMode       = "Always"

  res@mpMinLatF                   =  tkres@min_lat
  res@mpMaxLatF                   =  tkres@max_lat
  res@mpMinLonF                   =  tkres@min_lon
  res@mpMaxLonF                   =  tkres@max_lon
  res@mpOutlineBoundarySets       = "GeophysicalAndUSStates"
  res@mpDataBaseVersion           = "MediumRes"
  res@mpFillColors                = (/"background","LightBlue","SandyBrown","LightBlue", \
                                      "transparent"/)
  res@mpGeophysicalLineColor      = "Azure4"
  res@mpGeophysicalLineThicknessF = 0.5
  res@mpGridLineColor             = "Azure4"
  res@mpGridLineThicknessF        = 0.5
  res@mpGridMaskMode              = 3
  res@mpGridSpacingF              = 5
  res@mpLimbLineColor             = "Azure4"
  res@mpLimbLineThicknessF        = 0.5
  res@mpNationalLineColor         = "Azure4"
  res@mpNationalLineThicknessF    = 0.5
  res@mpUSStateLineColor          = "Azure4"
  res@mpUSStateLineThicknessF     = 0.5
  res@mpGridAndLimbOn             = True
  res@mpOutlineOn                 = True

;
; Tickmark stuff.
;
  res@tmXBMajorOutwardLengthF     = 0.0
  res@tmXBMajorLengthF            = 0.0
  res@tmXBLabelFontHeightF    =     0.007

end

;_______________________________________________________________________________________
; Add a title header topleft

undef("PlotHeaderL")
procedure PlotHeaderL(wks,map,tkres)
begin

  ;Add header at top
  txt0 = create "MainPlotTille" textItemClass wks
    "txString"              : tkres@header
    "txFontHeightF"         : 0.018
    "txBackgroundFillColor" : "White"
    "txPerimOn"             : True
    "txPerimColor"          : "Black"
    "txFont"                : "helvetica-bold"
  end create
  anno = NhlAddAnnotation(map,txt0)
  setvalues anno
    "amZone"           : 0
    "amSide"           : "Top"
    "amJust"           : "TopLeft"
    "amParallelPosF"   : -0.5
    "amOrthogonalPosF" : 0.5
    "amResizeNotify"   : False
  end setvalues

end

;_______________________________________________________________________________________
; Add time information topright

undef("PlotHeaderR")
procedure PlotHeaderR(wks,map,header,tkres)
begin

  ; Set return character
    cr = "~C~"

  ; Placement of the time info
    slon = tkres@max_lon - 0.2
    slat = tkres@max_lat - 0.2

    txres                       = True
    txres@txFont                = "helvetica-bold"
    txres@txFontColor           = "Black"
    txres@txFontHeightF         = 0.01
    txres@txPerimOn             = True
    txres@txPerimColor          = "Black"
    txres@txJust                = "TopRight"
    txres@txBackgroundFillColor = "White"

    if (tkres@plotTYPE .eq. 0) then
      timestring = "TIME:"
      do i = 0, dimsizes(header)-1
        timestring = timestring + cr + header(i)
      end do
    else
      timestring = "Forecasts:"
      do i = 0, dimsizes(header)-1
        ii = i+1
        timestring = timestring + cr + ii + ". " + header(i)
      end do
    end if

    gsn_text(wks,map,timestring,slon,slat,txres)

end


;_______________________________________________________________________________________
; Functions to plot different symbols

undef("create_hurricane_symbol")
function create_hurricane_symbol(wks)
begin

  mstring = "p"
  fontnum = 37
  xoffset = 0.0
  yoffset = 0.0
  ratio   = 1.0
  size    = 1.0
  angle   = 0.0

  new_index = NhlNewMarker(wks, mstring, fontnum, xoffset, yoffset, \
                           ratio, size, angle)

  return(new_index)
end

undef("create_diamond_symbol")
function create_diamond_symbol(wks)
begin

  mstring = "q"
  fontnum = 35
  xoffset = 0.0
  yoffset = 0.0
  ratio   = 1.0
  size    = 1.0
  angle   = 0.0

  new_index = NhlNewMarker(wks, mstring, fontnum, xoffset, yoffset, \
                           ratio, size, angle)

  return(new_index)
end


;_______________________________________________________________________________________
; get lat / lon / lev / wnd from input data
; input filenames are hardcorded - may need to change this later
; return an array (4,*), where the second dim is the number of
;    input data to be plotted.
;    The first dim contain the fields lat-lon-lev-wnd

undef("GetTrackData")
function GetTrackData(tkres)
begin

    if ( tkres@dataSOURCE .eq. "RIP4" ) then
      filename = tkres@filename

      ; Read lat/lon speed and pressure
      tmp = systemfunc("cut -c4-10 " + filename)
        lat = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData))
      tmp = systemfunc("cut -c13-20 " + filename)
        lon = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData))
      tmp = systemfunc("cut -c53-60 " + filename)
        lev = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData))
      tmp = systemfunc("cut -c65-71 " + filename)
        wnd = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData))
      ; need to know how many points we have
      fno = lat
      tkres@numGridPoints = dimsizes(lat)
      do it = 0,tkres@numGridPoints-1
        fno(it) = it
      end do

    end if

    if ( tkres@dataSOURCE .eq. "ARW" ) then
      filename = tkres@filename

      ; Read lat/lon speed and pressure
      lat = stringtofloat(systemfunc("cut -c29-34 " + filename))
      lon = stringtofloat(systemfunc("cut -c35-43 " + filename))
      lev = stringtofloat(systemfunc("cut -c44-50 " + filename))
      wnd = stringtofloat(systemfunc("cut -c51-57 " + filename))
      dummy = dimsizes(lat)
      fno = new( dummy , float )
      fno = 0.0
  
    end if

    npoints = dimsizes(lat)
    ll_lev_wnd = new( (/5,npoints/) , float )
    ll_lev_wnd(0,:) = lat
    ll_lev_wnd(1,:) = lon
    ll_lev_wnd(2,:) = lev
    ll_lev_wnd(3,:) = wnd
    ll_lev_wnd(4,:) = fno
    return (ll_lev_wnd)

end


;_______________________________________________________________________________________
; get date and time information for lables on plot

undef("Get_dd_hh")
function Get_dd_hh(fno,time_start_end,tkres)
begin

  if ( tkres@dataSOURCE .eq. "RIP4" ) then
    tm  = new (dimsizes(fno), string)

    MonthEnd = new ( 12 , integer )
    MonthEnd = (/ 32,29,32,31,32,31,32,32,31,32,31,32 /)
   
    yy = (systemfunc("echo " + time_start_end(0) + " | cut -c1-4" ) )
    mm = stringtointeger(systemfunc("echo " + time_start_end(0) + " | cut -c5-6" ) )
    dd = stringtointeger(systemfunc("echo " + time_start_end(0) + " | cut -c7-8" ) )
    hh = stringtointeger(systemfunc("echo " + time_start_end(0) + " | cut -c9-10" ) )

    fno = fno*tkres@interval
    fnoi = floattoint(fno)
    do i = 0,dimsizes(fnoi)-1
      if ( i .eq. 0 ) then
        hh = hh + fnoi(0)
      else
        hh = hh + fnoi(i) - fnoi(i-1)
      end if

      if ( hh .eq. 24 ) then
        dd = dd + 1
        hh = 0
        if ( dd .eq. MonthEnd(mm) ) then
          mm = mm + 1
          dd = 1
        end if
      end if

      if ( mm .lt. 10 ) then
        time_start_end(1) = yy + "0" + mm
      else
        time_start_end(1) = yy + mm
      end if
      if ( dd .lt. 10 ) then
        tm(i) = "0" + dd
        time_start_end(1) = time_start_end(1) + "0" + dd
      else
        tm(i) = dd
        time_start_end(1) = time_start_end(1) + dd
      end if
      if ( hh .lt. 10 ) then
        tm(i) = tm(i) + "/0" + hh
        time_start_end(1) = time_start_end(1) + "0" + hh
      else
        tm(i) = tm(i) + "/" + hh
        time_start_end(1) = time_start_end(1) + hh
      end if
    end do


  end if

  if ( tkres@dataSOURCE .eq. "ARW" ) then

    ; Read Time stamps
    dd = (systemfunc("cut -c15-16 " + tkres@filename))
    hh = (systemfunc("cut -c18-19 " + tkres@filename))
    tm = dd + "/" + hh
  
  end if

  return (tm)

end


;_______________________________________________________________________________________
; get interval of the input data            

undef("GetInterval")
function GetInterval(tkres)
begin
  if ( tkres@dataSOURCE .eq. "RIP4" ) then
      interval = tkres@dataINT
  end if
  if ( tkres@dataSOURCE .eq. "ARW" ) then
    ; This dataset has a full date, so we can calculate the interval of
    ;      incoming data
      filename = tkres@filename
      mins = stringtofloat(systemfunc("cut -c21-22 " + filename))
    ; As the first time may be wrong, we use the 2nd and 3rd times for calculations
      interval = mins(3) - mins(2)
    ; We need to know the answer as a fraction of an hour
      interval = interval / 60.
  end if
  print("Data interval is " + interval + " hour(s)" )
  return(interval)
end


;_______________________________________________________________________________________
; get start and end times to place on the map as information
; return an array (2), where the the first string is the start
;    time and the second one the end time

undef("GetTimeInfo")
function GetTimeInfo(tkres)
begin
  if ( tkres@dataSOURCE .eq. "RIP4" ) then
    tstart = tkres@startTIME
    tend = tkres@startTIME
  end if
  if ( tkres@dataSOURCE .eq. "ARW" ) then
    tstart= systemfunc("head -n 1 " + tkres@filename  + " | cut -c7-19 ")
    tend  = systemfunc("tail -n 1 " + tkres@filename  + " | cut -c7-19 ")
  end if
  time_start_end = new(2,string)
  time_start_end(0) = tstart
  time_start_end(1) = tend
  return(time_start_end)
end


;_______________________________________________________________________________________
; get DataType and Hurricane Categories
; return an array (6), where the first dimension is the DataType
;    and the other 5 is cat1 to cat5

undef("GetCatInfo")
function GetCatInfo(tkres)
begin
  if ( tkres@dataSOURCE .eq. "RIP4" ) then
    ; we know this data is in m/s
    DataType = "m/s"
    cat1 = 33.0
    cat2 = 42.0
    cat3 = 51.0
    cat4 = 59.0
    cat5 = 69.0
  end if
  if ( tkres@dataSOURCE .eq. "ARW" ) then
    ; we know this data is in kts
    DataType = "kts"
    cat1 = 64.0
    cat2 = 84.0
    cat3 = 97.0
    cat4 = 114.0
    cat5 = 135.0
  end if
  cat = new(6,string)
  cat(0) = DataType
  cat(1) = cat1
  cat(2) = cat2
  cat(3) = cat3
  cat(4) = cat4
  cat(5) = cat5
  return(cat)
end


;_______________________________________________________________________________________
; Plot best track if asked for
; currently only a standard plot, but this can be changed

undef("PlotBestTrack")
procedure PlotBestTrack(wks,map,tkres)
begin

    print("Now attempting to plot bestTRACK from data file " + tkres@btFILE)

    lat  = stringtofloat(systemfunc("cut -c18-21 " + tkres@btFILE))
    latsign  = systemfunc("cut -c22-22 " + tkres@btFILE)
    lon  = stringtofloat(systemfunc("cut -c24-28 " + tkres@btFILE))
    lonsign  = systemfunc("cut -c29-29 " + tkres@btFILE)
    bthh  = systemfunc("cut -c1-2 " + tkres@btFILE)
    btdd  = systemfunc("cut -c11-12 " + tkres@btFILE)
    btdate = btdd + "/" + bthh
    btcat = stringtoint(systemfunc("cut -c33-35 " + tkres@btFILE))

    ; Ckect to see if we have S/N ; W/E data
    checkDATA = 0
    do isign=0,dimsizes(lat)-1
       if (latsign(isign) .eq. "S") then
         lat(isign) = -1.0 * lat(isign)
       else
         if (latsign(isign) .ne. "N") then
           checkDATA = checkDATA + 1
         end if
       end if
       if (lonsign(isign) .eq. "W") then
         lon(isign) = -1.0 * lon(isign)
       else
         if (lonsign(isign) .ne. "E") then
           checkDATA = checkDATA + 1
         end if
       end if
    end do
    if ( checkDATA .gt. 0 ) then
       print("###########################################################")
       print("Looks like some of your best TRACK data is not in columns.")
       print("Please ensure all data is in columns and re-run the script.")
       print("###########################################################")
    end if


    lnres = True
    lnres@gsLineThicknessF = 2.6
    do i=0,dimsizes(lat)-2
         lnres@gsLineColor    = tkres@btCOLOR
      ;if ( lon(i) .lt. -75 ) then
         gsn_polyline(wks,map,(/lon(i),lon(i+1)/),(/lat(i),lat(i+1)/),lnres)
      ;end if
    end do

    print("   We have " + dimsizes(lat) + " data points")
    print("   Track start point (lon/lat) = "+ lon(0) +" "+ lat(0))
    print("   Track end point (lon/lat)   = "+ lon(dimsizes(lat)-1) +" "+ lat(dimsizes(lat)-1))

    txres2               = True
    txres2@txFont        = "helvetica-bold"
    txres2@txFontHeightF = 0.008
    txres2@txJust        = "CenterLeft"
    txres2@txAngleF        = 45.00
    lnres2                  = True
    lnres2@gsLineColor      = "Black"
    lnres2@gsLineThicknessF = 0.7
    pmres = True
    pmres@gsMarkerIndex = create_hurricane_symbol(wks)
    do i=0,dimsizes(lat)-1
      pmres@gsMarkerColor = "Black"
      pmres@gsMarkerIndex = 16
      pmres@gsMarkerSizeF = 0.005
      if (btcat(i) .ge. 74) then
        pmres@gsMarkerColor = 5
      end if
      if (btcat(i) .ge. 96) then
        pmres@gsMarkerColor = 6
      end if
      if (btcat(i) .ge. 111) then
        pmres@gsMarkerColor = 7
      end if
      if (btcat(i) .ge. 131) then
        pmres@gsMarkerColor = 8
      end if
      if (btcat(i) .ge. 155) then
        pmres@gsMarkerColor = 9
      end if

      if (btcat(i) .ge. 74) then
        pmres@gsMarkerIndex = create_hurricane_symbol(wks)
        pmres@gsMarkerSizeF = 0.01
      end if

      gsn_polyline(wks,map,(/lon(i),lon(i)+1.0/),(/lat(i),lat(i)+1.0/),lnres2)
      gsn_polymarker(wks,map,lon(i)+1.0,lat(i)+1.0,pmres)
      if ( bthh(i) .eq. "03" )then
        gsn_text(wks,map,btdate(i),lon(i)+1.2,lat(i)+1.2,txres2)
      end if

    end do  

end


;_______________________________________________________________________________________
; Get some line number if we are dealing with RIP4 data
; This is to strip out all junk from the input data
;
; line_breaks will hold 3 pieces of information
;    start of good data         
;     end of good data            
;     number of lines in file
;  

undef("GetLineNums")
function GetLineNums(tkres)
begin

    if ( tkres@dataSOURCE .eq. "RIP4" ) then
      filename = tkres@filename

      s_tester = (systemfunc("cut -c4-10 " + filename))
      tester = new(dimsizes(s_tester),float)

      line_breaks = new(3,integer)
      line_breaks    = 0
      line_breaks(1) = dimsizes(s_tester) - 1
      line_breaks(2) = dimsizes(s_tester)
      tkres@bad_lines = False
      found_start = False
      do i =0,dimsizes(s_tester)-1

         tester(i) = stringtofloat(s_tester(i))
         if ( .not. tkres@bad_lines .and. ismissing(tester(i)) ) then
           tkres@bad_lines = True
         end if
         if ( .not. found_start .and. .not. ismissing(tester(i)) ) then
           line_breaks(0) = i
           found_start = True
         end if
         if ( found_start .and. ismissing(tester(i)) ) then
           line_breaks(1) = i - 1
           return (line_breaks)
         end if
      end do

    end if

    return (line_breaks)

end


;_______________________________________________________________________________________
; Draw the track line on the map that we have created previously

undef("PlotTrackLine")
procedure PlotTrackLine(wks,map,lat,lon,wnd,cat,tkres)
begin
  
  ; We are going to draw the line every 3 hour, or every timestep,
  ; Depending on frequency of input data
    line_int = floattoint(3/tkres@interval)   
    if ( line_int .le. 1 )
      line_int = 1
    end if

  ; Line resources
    lnres                  = True
    lnres@gsLineColor      = "Black"
    lnres@gsLineThicknessF = 2.7

    if (tkres@plotTYPE .eq. 4)
      lnres@gsLineColor      = tkres@StartLineColor + tkres@iNumFil
      lnres@gsLineThicknessF = 3.0
    end if

  ; Draw the line
    do i=0,dimsizes(wnd)-line_int-1,line_int
      if ( wnd(i) .ne. -1.0 .and. wnd(i+line_int) .ne. -1.0 ) then
        if (tkres@plotTYPE .eq. 3)
          lnres@gsLineColor      = "Azure4"
          if (wnd(i) .gt. stringtofloat(cat(1)))
            lnres@gsLineColor    = 5
          end if
          if (wnd(i) .gt. stringtofloat(cat(2)))
            lnres@gsLineColor    = 6
          end if
          if (wnd(i) .gt. stringtofloat(cat(3)))
            lnres@gsLineColor    = 7
          end if
          if (wnd(i) .gt. stringtofloat(cat(4)))
            lnres@gsLineColor    = 8
          end if
          if (wnd(i) .gt. stringtofloat(cat(5)))
            lnres@gsLineColor    =9
          end if
        end if
        gsn_polyline(wks,map,(/lon(i),lon(i+line_int)/),(/lat(i),lat(i+line_int)/),lnres)
      end if
    end do

end


;_______________________________________________________________________________________
; Draw pressure and wind information to left of track line      
; Plot date information to the right of the track line

undef("PlotTrackInfo")
procedure PlotTrackInfo(wks,map,lat,lon,wnd,wndlev,dd_hh,tkres)
begin
  
  ; We are going to plot the information every 12 hours, or every timestep,
  ; Depending on frequency of input data
    label_int= floattoint(12/tkres@interval)   
    if ( label_int .le. 1 )
      label_int = 1
    end if

    if (tkres@plotTYPE .le. 1)
    ; Plot the Pressure/Wind Info to left of track line
      txres               = True
      txres@txFont        = "helvetica-bold"
      txres@txFontHeightF = 0.009
      txres@txJust        = "CenterRight"
      txres@txAngleF      = tkres@textANGLE
      do i=0,dimsizes(wnd)-1,label_int
        if ( wnd(i) .ne. -1.0) then
          gsn_text(wks,map,wndlev(i),lon(i),lat(i),txres)
        end if
      end do
      delete(txres)
    end if

  ; Plot date information to right of track line
    txres               = True
    txres@txFont        = "helvetica-bold"
    txres@txJust        = "CenterLeft"
    txres@txAngleF      = tkres@textANGLE

    dis_marker = 0.0
    if (tkres@plotTYPE .le. 1)
      txres@txBackgroundFillColor = "White"
      txres@txPerimOn             = True
      txres@txPerimColor          = "Black"
      txres@txFontHeightF         = 0.008
      dis_marker = 0.3   ; how far must time stamp be drawn from marker
      if (tkres@textANGLE .gt. 25.)
        dis_marker = 0.1   ; how far must time stamp be drawn from marker
      end if
    end if

    if (tkres@plotTYPE .ge. 2)
      txres@txFontHeightF         = 0.007
      dis_marker = 0.2   ; how far must time stamp be drawn from marker
      label_int= floattoint(24/tkres@interval)   ; we will only plot these every 24 hours
    end if

    lon_dd_hh = lon + dis_marker
    lat_dd_hh = lat + dis_marker

    do i=0,dimsizes(wnd)-1,label_int
      if ( wnd(i) .ne. -1.0) then
        gsn_text(wks,map,dd_hh(i),lon_dd_hh(i),lat_dd_hh(i),txres)
      end if
    end do
    delete(txres)

  ; Add forecast numbers
    if (tkres@plotTYPE .ge. 2)
      txres                       = True
      txres@txFont                = "helvetica-bold"
      txres@txJust                = "BottomLeft"
      txres@txBackgroundFillColor = "White"
      txres@txPerimOn             = True
      txres@txPerimColor          = "Black"
      txres@txFontHeightF         = 0.008
      forecastnumber = 0
      do i=0,dimsizes(wnd)-1,floattoint(3/tkres@interval) ; this is to find end of track line
        itest = dimsizes(wnd)-1-i
        if ( forecastnumber .eq. 0 .and. wnd(itest) .ne. -1.0 ) then
          gsn_text(wks,map,tkres@iNumFil+1,lon(itest),lat(itest),txres)
          forecastnumber = 1
        end if
      end do

    end if  

end


;_______________________________________________________________________________________
; Draw Sysmbols over the track line        

undef("PlotTrackSymbols")
procedure PlotTrackSymbols(wks,map,lat,lon,wnd,cat,tkres)
begin

  ; We are going to draw the symbols every 6 hours, or every timestep,
  ; Depending on frequency of input data
    symbol_int= floattoint(6/tkres@interval)   
    if ( symbol_int .le. 1 )
      symbol_int = 1
    end if

    pmres               = True

    if (tkres@plotTYPE .ge. 3 ) then
      do i=0,dimsizes(wnd)-1,symbol_int
        if ( wnd(i) .ne. -1.0) then
          pmres@gsMarkerColor = "White"
          pmres@gsMarkerIndex = 16
          pmres@gsMarkerSizeF = 0.003
          gsn_polymarker(wks,map,lon(i),lat(i),pmres)
        end if
      end do
      return
    end if
  
    do i=0,dimsizes(wnd)-1,symbol_int
      if ( wnd(i) .ne. -1.0) then
        pmres@gsMarkerColor = "Black"
        if ( wnd(i) .gt. stringtofloat(cat(1))) then
          pmres@gsMarkerIndex = create_hurricane_symbol(wks)
          pmres@gsMarkerSizeF = 0.01
        else
          ;pmres@gsMarkerIndex = 4
          pmres@gsMarkerIndex = 16
          pmres@gsMarkerSizeF = 0.006
        end if
        if (wnd(i) .gt. stringtofloat(cat(1))) then
          pmres@gsMarkerColor = "Blue"
        end if
        if (wnd(i) .gt. stringtofloat(cat(2))) then
          pmres@gsMarkerColor = "DarkGreen"
        end if
        if (wnd(i) .gt. stringtofloat(cat(3))) then
          pmres@gsMarkerColor = "Yellow"
        end if
        if (wnd(i) .gt. stringtofloat(cat(4))) then
          pmres@gsMarkerColor = "Red"
        end if
        if (wnd(i) .gt. stringtofloat(cat(5))) then
          pmres@gsMarkerColor = "Purple3"
        end if
        gsn_polymarker(wks,map,lon(i),lat(i),pmres)
      end if
    end do

end


;_______________________________________________________________________________________
; Add Hurricane sysmbols (color) and cat's bottomleft

undef("HurricaneLegend")
procedure HurricaneLegend(wks,map,tkres)
begin

    if (tkres@plotTYPE .eq. 4)
      return
    end if

  ; Placement of the legdends
    slon = tkres@min_lon + 0.5
    slat = tkres@min_lat + 0.5

    pmres = True
    pmres@gsMarkerSizeF = 0.013
    pmres@gsMarkerIndex = create_hurricane_symbol(wks)

    txres               = True
    txres@txFont        = "helvetica-bold"
    txres@txFontColor   = "White"
    txres@txFontHeightF = 0.008

    start_color = 5   ; our color table has Blue as color number 5
    do j=0,4          ; plot 5 color symbols
      pmres@gsMarkerColor = start_color + j
      gsn_polymarker(wks,map,slon,slat,pmres)
      gsn_text(wks,map,1+j,slon,slat,txres)
      slon = slon + 0.5
    end do

    ; Draw a box around the legend
    lnres                  = True
    lnres@gsLineColor      = "Black"
    lnres@gsLineThicknessF = 1.7
    slat = slat + 0.4
    gsn_polyline(wks,map,(/tkres@min_lon,slon/),(/slat+0.4,slat+0.4/),lnres)
    gsn_polyline(wks,map,(/slon,slon/),(/slat+0.4,tkres@min_lat/),lnres)
    slon = slon - 0.5 - 0.5 - 0.5

    ; Add some text
    txres@txFontColor   = "Black"
    txres@txFontHeightF = 0.009
    txres@txJust        = "BottomCenter"
    gsn_text(wks,map,"hur cat:",slon,slat,txres)

end


;_______________________________________________________________________________________
; Add pressure/wind and dd/hh legend to bottomright

undef("SymbolLegend")
procedure SymbolLegend(wks,map,tkres)
begin

    if ( tkres@plotTYPE .ge. 2)
      return
    end if
  ; Set return character
    cr = "~C~"

  ; Placement of the ledgends
    slon = tkres@max_lon - 2.0
    slat = tkres@min_lat + 0.5

    pmres               = True
    pmres@gsMarkerSizeF = 0.008
    pmres@gsMarkerIndex = create_hurricane_symbol(wks)
    pmres@gsMarkerColor = "Black"
    gsn_polymarker(wks,map,slon,slat,pmres)

    txres               = True
    txres@txFont        = "helvetica-bold"
    txres@txFontHeightF = 0.008
    txres@txJust        = "CenterRight"
    sslon = slon - 0.5
    left = "pressure" + cr + "max wnd ("+tkres@DataType+")"
    gsn_text(wks,map,left,sslon,slat,txres)

    txres@txBackgroundFillColor = "White"
    txres@txPerimOn             = True
    txres@txJust                = "CenterLeft"
    right = "dd/hh"
    slon = slon + 0.5
    gsn_text(wks,map,right,slon,slat,txres)

end


;_______________________________________________________________________________________
;_______________________________________________________________________________________
;_______________________________________________________________________________________
; Main script

begin

  ;----------------------------------------------------------------------
  ; Get dataSOURCE and file info from command line
  ;----------------------------------------------------------------------
      tkres = True

      ; Make sure we have a datafile to work with
      if (.not. isvar("inputFILE") ) then
        print(" ")
        print(" ### MUST SUPPLY a inputFILE ### ")
        print(" ")
        print("     Something like: ")
        print("     ncl CreateTracks.ncl inputFILE=track_2005082800_12k stormNAME=Katrina dataSOURCE=RIP4 wksTYPE=X11 ")
        print("          REMEMBER TO ADD QUOTES" )
        print("          Refer to the information at the top od this file for more info and syntax" )
        exit
      end if
      tkres@inputFILE = inputFILE

    ; dataSOURCE = ARW (data from WRF-ARW model)
    ; dataSOURCE = RIP4 (data from RIP)
      if (.not. isvar("dataSOURCE") ) then
        dataSOURCE = "ARW"
      end if
      if (dataSOURCE .ne. "RIP4" .and. dataSOURCE .ne. "ARW") then
        print("'" + dataSOURCE + "'" + " is not a valid dataSOURCE")
        print("  Current options are :")
        print("  RIP4 for data generated by rip4")
        print("  ARW for data generated by the WRF-ARW model")
        exit
      end if
      tkres@dataSOURCE = dataSOURCE
      print("Input is " + tkres@dataSOURCE + " type data")

    ; Do we have any special data that comes in at a different interval from default
    ; As we only care about this for RIP4 data (ARW has all the info in the file)
    ; We are going to make the default 3h
      if (.not. isvar("dataINT") .or. dataINT .eq. 0.0) then
         dataINT = 3.0
        if (dataSOURCE .eq. "RIP4" ) then
          print("   Using a default of " +dataINT+ " hours between data points" )
          print("   If this is not correct, use dataINT on the command line, to override" )
        end if
      end if
      tkres@dataINT = dataINT

  ;----------------------------------------------------------------------
  ; Are we dealing with SUMMARIES
  ;----------------------------------------------------------------------
      if (.not. isvar("plotTYPE") ) then
         plotTYPE = 0
      end if
      if (plotTYPE .eq. 0) then
         filename = inputFILE
         tkres@plotTYPE = 0
         print("Doing a single time plot")
      else
         filename   = (systemfunc("cut -d':' -f1 " + inputFILE))
         filedates   = (systemfunc("cut -d':' -f2 " + inputFILE))
         TRHeader = new(dimsizes(filename),string)
         tkres@plotTYPE = plotTYPE
         print("Doing summary plot type " + tkres@plotTYPE)
         if ( dimsizes(filename) .gt. 10 ) then
            print("########################################################")
            print("Number of input files for summary is: " + dimsizes(filename) )
            print("If code is hanging or this looks wrong,")
            print("make sure you are not pointing to a track data file,")
            print("but to a file containing a list if track data files.")
            print("This will be true even for one track plot.")
            print("########################################################")
         end if
      end if

    ; Do we have any special data that does not contain date info
    ; As we only care about this for RIP4 data (ARW has all the info in the file)
    ; No default - we need this set        
      if (dataSOURCE .eq. "RIP4" .and. plotTYPE .eq. 0 .and. .not. isvar("startTIME") ) then
        print(" ")
        print(" ### MUST SUPPLY a startTIME for RIP4 data ### ")
        print(" ")
        exit
      end if
      if (isvar("startTIME") ) then
        tkres@startTIME = startTIME
      end if

  ;----------------------------------------------------------------------
  ; Basic Setup                  
  ;----------------------------------------------------------------------

  ; open the workstation
    if (.not. isvar("wksTYPE") ) then
      wksTYPE = "pdf"
    end if
    tkres@wksTYPE = wksTYPE
    if (.not. isvar("outputFILE") ) then
      outputFILE = "hur_track"
    end if
    tkres@outputFILE = outputFILE
    wks = gsn_open_wks(wksTYPE,outputFILE)  

  ; set up header
    if (.not. isvar("stormNAME") ) then
      header = "ARW Forecast"
      stormNAME = "UnKnown"
    else
      header = "ARW Forecast: " + stormNAME
    end if
    tkres@header = header
    tkres@stormNAME = stormNAME

    ; are we doing the best track
    if (.not. isvar("bestTRACK") ) then
      bestTRACK = False
    end if
    tkres@bestTRACK = bestTRACK
    if (.not. isvar("btFILE") ) then
      btFILE = "BestTrack"
    end if
    tkres@btFILE = btFILE
    if (.not. isvar("btCOLOR") ) then
      btCOLOR = "Red"
    end if
    tkres@btCOLOR = btCOLOR

    ; less common settings
    ; rotate the text angle on the track line
    if (.not. isvar("textANGLE") ) then
      textANGLE=0.0
    end if
    tkres@textANGLE = textANGLE
    ; rotate the output .pdf file
    if (.not. isvar("rot") ) then
      rot=-90
    end if
    tkres@rotPDF = rot

    ; change the start number for summary 4 plots
    if (.not. isvar("StartLineColor") ) then
      StartLineColor=5
    end if
    tkres@StartLineColor = StartLineColor


  ; Set return character
    cr = "~C~"


  ;----------------------------------------------------------------------
  ; Read storm track and info
  ;----------------------------------------------------------------------

    do iNumFil=0,dimsizes(filename)-1

      tkres@iNumFil  = iNumFil
      tkres@filename = filename(iNumFil)
      if ( plotTYPE .ne. 0 ) then
        tkres@startTIME = filedates(iNumFil)
      end if

      if ( dataSOURCE .eq. "RIP4" ) then   ; need to clean up the datafile
       ; First figure out some line numbers
         line_breaks = GetLineNums(tkres)
         tkres@StartGoodData = line_breaks(0)
         tkres@EndGoodData   = line_breaks(1)
         if ( tkres@bad_lines ) then
           print(" ")
           print("Warning messages above this point is fine, and can be ignored")
           print("=============================================================")
           print(" ")
         end if
      end if
      ll_lev_wnd     = GetTrackData(tkres)
                       lat = ll_lev_wnd(0,:)
                       lon = ll_lev_wnd(1,:)
                       lev = ll_lev_wnd(2,:)
                       wnd = ll_lev_wnd(3,:)
                       fno = ll_lev_wnd(4,:)
      interval       = GetInterval(tkres)
                       tkres@interval = interval
      time_start_end = GetTimeInfo(tkres)
      dd_hh          = Get_dd_hh(fno,time_start_end,tkres)
      cat            = GetCatInfo(tkres)
                       tkres@DataType = cat(0)

  ; set up the data the way we want to plot it
    slev   = sprintf("  %4.0f",lev)
    swnd   = sprintf("  %5.1f",wnd)
    wndlev = slev + "    " + cr + swnd + "    "  

    ; Create top right header
    if (tkres@plotTYPE .eq. 0) then
      TRHeader = time_start_end
    else
      TRHeader(iNumFil) = time_start_end(0)
    end if

  ;----------------------------------------------------------------------
  ; Set up the window
  ;----------------------------------------------------------------------

    ; Set user over ride of plot domain to true. If user does not use
    ;     the override commands, this will be set to False again
    tkres@latMINset = True
    tkres@latMAXset = True
    tkres@lonMINset = True
    tkres@lonMAXset = True
    if (iNumFil .eq. 0) then
      ; defaults if not overwritten on the command line
      if (.not. isvar("latMIN") ) then
        latMIN = 25.
        tkres@latMINset = False
      end if
      tkres@latMIN = latMIN
      if (.not. isvar("latMAX") ) then
        latMAX = 40.
        tkres@latMAXset = False
      end if
      tkres@latMAX = latMAX
      if (.not. isvar("lonMIN") ) then
        lonMIN = -90.
        tkres@lonMINset = False
      end if
      tkres@lonMIN = lonMIN
      if (.not. isvar("lonMAX") ) then
        lonMAX = -75.
        tkres@lonMAXset = False
      end if
      tkres@lonMAX = lonMAX
      if (.not. isvar("disEDGE") ) then
        disEDGE = 2.
      end if
      tkres@disEDGE= disEDGE

      ; Make sure the final window is bigger than the data spread
        MapWindow(lat,lon,wnd,tkres)
  
      print("The Track Plot Window is set to")
      print("   Latitute  from " + tkres@min_lat + " to " + tkres@max_lat)
      print("   Longitute from " + tkres@min_lon + " to " + tkres@max_lon)
      print("The buffer between the data and the plot edge is " + tkres@disEDGE + " degrees")

    ;----------------------------------------------------------------------
    ; Set up resources and create map
    ;----------------------------------------------------------------------

    ; Get colors to use
      track_colors = TrackColors()
      gsn_define_colormap(wks,track_colors)
      ;gsn_draw_colormap(wks)

    ; Set some map resources and zoom in           
      res = False
      MapResources(res,tkres)

    ; Create the map.
      map = gsn_csm_map(wks,res)
    ; Add Header top left
      PlotHeaderL(wks,map,tkres)

      draw(map)

    end if

  ;----------------------------------------------------------------------
  ; The plotting section.
  ;----------------------------------------------------------------------

  ; Draw the track line
    PlotTrackLine(wks,map,lat,lon,wnd,cat,tkres)

  ; Add Symbols over line
    PlotTrackSymbols(wks,map,lat,lon,wnd,cat,tkres)

  ; Add pressure/wind information on track line
  ; Add date information on track line
  ; Only when needed - determine by plotTYPE
    PlotTrackInfo(wks,map,lat,lon,wnd,wndlev,dd_hh,tkres)

    delete(ll_lev_wnd)
    delete(interval)
    delete(fno)
    delete(cat)
    delete(lat)
    delete(lon)
    delete(lev)
    delete(wnd)
    delete(slev)
    delete(swnd)
    delete(wndlev)
    delete(dd_hh)

  end do


  ; Add time information top right
    PlotHeaderR(wks,map,TRHeader,tkres)

  ; Add hurricane legend in bottom left
    HurricaneLegend(wks,map,tkres)

  ; Add pressure/wind and dd/hh legend to bottom right
    SymbolLegend(wks,map,tkres)


  ; Finally - are we doing a best track
    if (tkres@bestTRACK) then
      PlotBestTrack(wks,map,tkres)
    end if


  ; We are done
  frame(wks)
  delete(wks)


end