*----------------------------------------------------------------------
*
*  Adapted from plot_hurricane_center2a.gs, which was adapted from:
*
*  (1) "plot_hurricane_center2.gs" (from "plot_hurricane_center.gs"),
*      locally modified versions derived from Joe Covert's script
*      "plot_tc_shi.gs" (Joe.Covert@noaa.gov), which was posted to the
*      GrADS Listserv on 9/13/2004 (his script was named plot_tc_shi.gs)
*      on 4/27/06, found via a Google search for "plot_tc_shi.gs":
*      http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs
*
*  (2) code in "hurricane_tracking.txt" (taken from an email exchange
*      posted to the GrADS Listserv by Kun-Hsuan Chou and
*      Arturo Caracas Uribe in Oct, 2004)
*
*
*  To use: grads -l
*
*----------------------------------------------------------------------

function main()

  'reinit'

************************************************************************
*  Allow external GrADS functions
************************************************************************

  rc = gsfallow("on")

************************************************************************
*  Define dataset file information and plotting limits for base map...
************************************************************************

  model = 'isabel_ras'
  psfc_file = '/net/cjs/stg/'model'/'model'.nc'
  wind_file = '/net/cjs/stg/'model'/windr-'model'.ctl'

  lon1 = 298 ; lon2 = 308
  lat1 =  20 ; lat2 =  25

************************************************************************
*  Legend captions, colors, and marks for the tracks
************************************************************************

  ntracks = 3
  track.1 = 'psocean min'       ; open_mod_file.1 = 'sdfopen 'psfc_file''
    color.1 = 2 ; mark.1 = 3 ; line.1 = 1
  track.2 = 'rel vorticity max' ; open_mod_file.2 = 'open 'wind_file''
    color.2 = 4 ; mark.2 = 6 ; line.2 = 5
  track.3 = 'wind speed max'    ; open_mod_file.3 = 'open 'wind_file''
    color.3 = 3 ; mark.3 = 2 ; line.3 = 2
*
*  Define the legend plotter input...
*
  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
  n=1
  while ( n <= ntracks )
    titles = ''titles' "'track.n'"'
    colors = ''colors' 'color.n
    marks = ''marks' "'mark.n'"'
    lines = ''lines' "'line.n'"'
    n=n+1
  endwhile

  legend_info = ''colors' 'lines' 'marks' 'titles' -p'

*  ...other drawing primitives

  marksize = 0.1 ; lthick = 5

************************************************************************
*  Use "cbar_line" or "cbar_line_box"?
************************************************************************

  legend_plotter = 'cbar_line'
*  legend_plotter = 'cbar_line_box'

************************************************************************
*  Plot page limits
************************************************************************

  plot_vpage = '0 11 0 8.5'
  plot_area  = '1 10 1 7.5'

  xl = subwrd(plot_area,1) ; xr = subwrd(plot_area,2)
  yb = subwrd(plot_area,3) ; yt = subwrd(plot_area,4)

  xc = xl + (xr-xl)/2

************************************************************************
*  GrADS metafile output...
************************************************************************

  'enable print hurricane_track-'model'.gx'

*-----------------------------------------------------------------------
*  Define and plot tracks...
*-----------------------------------------------------------------------

  n=1

  while ( n <= ntracks )

*  Open and define dataset file information...

    open_mod_file.n
    'q file'
    rec3 = sublin(result,3) ; binary_file = subwrd(rec3,2)
    say
    say
    say '*** 'binary_file' ***'
    say
    rec5 = sublin(result,5)
    nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt = subwrd(rec5,12)

*  Define base map using specified lon/lat...

    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'q dims'
    say
    say result
    rec2 = sublin(result,2)
    l1 = subwrd(rec2,6) ; l2 = subwrd(rec2,8)
    say '    (l1,l2) = ('l1','l2')'
    x1r = subwrd(rec2,11) ; x2r = subwrd(rec2,13)
    say '    (x1r,x2r) = ('x1r','x2r')'
    x1 = math_int(x1r+0.5) ; x2 = math_int(x2r+0.5)
    say '    (x1,x2) = ('x1','x2')'
    nxp = x2-x1+1
    say '    nxp = 'nxp
    rec3 = sublin(result,3)
    l1 = subwrd(rec3,6) ; l2 = subwrd(rec3,8)
    say '    (l1,l2) = ('l1','l2')'
    y1r = subwrd(rec3,11) ; y2r = subwrd(rec3,13)
    say '    (y1r,y2r) = ('y1r','y2r')'
    y1 = math_int(y1r+0.5) ; y2 = math_int(y2r+0.5)
    say '    (y1,y2) = ('y1','y2')'
    nyp = y2-y1+1
    say '    nyp = 'nyp
    say

    if( n = 1 )

*  Plot base map...

      'set vpage 'plot_vpage
      'set parea 'plot_area
      'set mproj latlon'
      'set mpt * 1 1 5'
      'set mpdset hires'
      'set poli on'
*      'set grid on'
      'set grads off'
      'draw map'

*  add lon/lat lines...

      'set ccolor 1'
      'set cstyle 5' ; 'd lon'
      'set cstyle 5' ; 'd lat'

*  add title, x- and y-labels...

      'draw title \\ 'model
      'draw xlab Longitude \\'
      'draw ylab \\ Latitude'

    endif

* Define locations of "hurricane" for plotting points

    if( track.n = 'psocean min' )
      ts=1
      tt=ts
      while ( tt <= nt )
        'set x 'x1' 'x2
        'set y 'y1' 'y2
        'set z 1'
        'set t 'tt
        'define psocean = maskout(ps,-topog)'
        'set x 'x1
        'set y 'y1
        'set z 1'
        'set t 'tt
        say ' tt = 'tt
        'set gxout print'
        'd minloc(min(psocean,y='y1',y='y2'),x='x1',x='x2')'
          rec=sublin(result,nxp+3)
          xc=subwrd(rec,1)
        'd minloc(min(psocean,x='x1',x='x2'),y='y1',y='y2')'
          rec=sublin(result,nyp+3)
          yc=subwrd(rec,1)
        say '    (xc,yc) = ('xc','yc')'
        'set x 'xc
        lonval = subwrd(result,4)
        'set y 'yc
        latval = subwrd(result,4)
        say '    (lonval,latval) = ('lonval','latval')'
        'q w2xy 'lonval' 'latval
        xpos.tt = subwrd(result,3)
        ypos.tt = subwrd(result,6)
        say '    (xpos,ypos) = ('xpos.tt','ypos.tt')'
***        pause()
        'undefine psocean'
        tt=tt+1
      endwhile
    endif
    if( track.n = 'rvor max' )
      ts=1
      tt=ts
      while ( tt <= nt )
        'set x 'x1' 'x2
        'set y 'y1' 'y2
        'set z 1'
        'set t 'tt
        'define rvor = hcurl(uur,vvr)'
        'set x 'x1
        'set y 'y1
        'set z 1'
        'set t 'tt
        say ' tt = 'tt
        'set gxout print'
        'd maxloc(max(rvor,y='y1',y='y2'),x='x1',x='x2')'
          rec=sublin(result,nxp+3)
          xc=subwrd(rec,1)
        'd maxloc(max(rvor,x='x1',x='x2'),y='y1',y='y2')'
          rec=sublin(result,nyp+3)
          yc=subwrd(rec,1)
        say '    (xc,yc) = ('xc','yc')'
        'set x 'xc
        lonval = subwrd(result,4)
        'set y 'yc
        latval = subwrd(result,4)
        say '    (lonval,latval) = ('lonval','latval')'
        'q w2xy 'lonval' 'latval
        xpos.tt = subwrd(result,3)
        ypos.tt = subwrd(result,6)
        say '    (xpos,ypos) = ('xpos.tt','ypos.tt')'
***        pause()
        'undefine rvor'
        tt=tt+1
      endwhile
    endif
    if( track.n = 'wind speed max' )
      ts=1
      tt=ts
      while ( tt <= nt )
        'set x 'x1' 'x2
        'set y 'y1' 'y2
        'set z 1'
        'set t 'tt
        'define wspd = mag(uur,vvr)'
        'set x 'x1
        'set y 'y1
        'set z 1'
        'set t 'tt
        say ' tt = 'tt
        'set gxout print'
        'd maxloc(max(wspd,y='y1',y='y2'),x='x1',x='x2')'
          rec=sublin(result,nxp+3)
          xc=subwrd(rec,1)
        'd maxloc(max(wspd,x='x1',x='x2'),y='y1',y='y2')'
          rec=sublin(result,nyp+3)
          yc=subwrd(rec,1)
        say '    (xc,yc) = ('xc','yc')'
        'set x 'xc
        lonval = subwrd(result,4)
        'set y 'yc
        latval = subwrd(result,4)
        say '    (lonval,latval) = ('lonval','latval')'
        'q w2xy 'lonval' 'latval
        xpos.tt = subwrd(result,3)
        ypos.tt = subwrd(result,6)
        say '    (xpos,ypos) = ('xpos.tt','ypos.tt')'
***        pause()
        'undefine wspd'
        tt=tt+1
      endwhile
    endif

*  Plotting track...

*  Specify drawing primitives

    marktype = mark.n
    lcolor = color.n
    lstyle = line.n

    say ' marktype, marksize, lcolor, lstyle and lthick:'
    say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick

    tt=ts
    xc = xpos.tt
    yc = ypos.tt
    'set line 'lcolor' 'lstyle' 'lthick
    'draw mark 'marktype' 'xc' 'yc' 'marksize
    tt = tt+1
    while ( tt <= nt )
      xn = xpos.tt
      yn = ypos.tt
      say ' 'xc' 'yc' 'xn' 'yn
      'draw line 'xc' 'yc' 'xn' 'yn
      'draw mark 'marktype' 'xn' 'yn' 'marksize
      tt=tt+1
      xc = xn
      yc = yn
    endwhile

    'close 1'

    pause()

    n=n+1
  endwhile

  say
  say '...plot legend'
  ''legend_plotter' -x 3 -y 6 'legend_info
  say

  'print'
  'disable print'

  say
  say '************************************************'
  say ' Finished with this script.'
  say '************************************************'
  say

return
