*----------------------------------------------------------------------
*
*  Adapted from plot_hurricane_track.gs, which was adapted from:
*
*  (1) "plot_hurricane_center2[a,c].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 = 1
  track_type.1 = 'psfc min'       ; open_mod_file.1 = 'sdfopen 'psfc_file''
    color.1 = 2 ; mark.1 = 3 ; line.1 = 1
    track_units.1 = 'hPa' ; label_scale.1 = 0.01 ; label_format.1 = '%6.1f'
  track_type.2 = 'rel vorticity max' ; open_mod_file.2 = 'open 'wind_file''
    color.2 = 4 ; mark.2 = 6 ; line.2 = 5
    track_units.2 = 's`a-1`n' ; label_scale.2 = 1 ; label_format.2 = '%8.3e'
  track_type.3 = 'wind speed max'    ; open_mod_file.3 = 'open 'wind_file''
    color.3 = 3 ; mark.3 = 2 ; line.3 = 2
    track_units.3 = 'm s`a-1`n' ; label_scale.3 = 1 ; label_format.3 = '%6.1f'
*
*  Define the legend plotter input...
*
  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
  n=1
  while ( n <= ntracks )
    titles = ''titles' "'track_type.n' ('track_units.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 cyclone_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

*  Plot base map...

    plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2,model)

* Define locations of "cyclone" for plotting points

    if( track_type.n = 'psfc min' )
      say
      say '--- 'track_type.n ' ---'
      say
      ts=1
      tt=ts
      while ( tt <= nt )
        'set t 'tt
        say
        say '---------'
        say ' tt = 'tt
        say '---------'
        say
        'set lon 'lon1' 'lon2
        'set lat 'lat1' 'lat2
        'set z 1'
        'define psfc = ps'
*
*  display surface pressure variable for interactive identification
*  of cyclone center... (substitute name of your surface pressure variable; a
*  variable that is defined over land and ocean...)
*
        'set gxout shaded'
        'd psfc'
*
*  following interactive "box location" code adapted from "cbar_line_box.gs"...
*
        say 'Click where you want the left upper corner of the box'
        'query bpos'
        xb1o = subwrd(result,3)
        yb2o = subwrd(result,4)
        say 'Click where you want the right lower corner of the box'
        'query bpos'
        xb2o = subwrd(result,3)
        yb1o = subwrd(result,4)

        say
        say '...left upper corner of the box at X Y: 'xb1o' 'yb2o
        say '...right lower corner of the box at X Y: 'xb2o' 'yb1o
*
*  draw box to show specified area for defining cyclone center...
*
        'set line 1 1'
        'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
*
*  convert box "xy" coordinates to "grid" coordinates
*  for finding cyclone center...
*
        'q xy2gr 'xb1o' 'yb2o
        x1o = subwrd(result,3) ; x1 = math_int(x1o+0.5)
        y2o = subwrd(result,6) ; y2 = math_int(y2o+0.5)
        'q xy2gr 'xb2o' 'yb1o
        x2o = subwrd(result,3) ; x2 = math_int(x2o+0.5)
        y1o = subwrd(result,6) ; y1 = math_int(y1o+0.5)
        say
        say 'Specified grid coordinates (x1o,x2o), (y1o,y2o) from box,'
        say 'rounded to (x1,x2), (y1,y2) for finding minimum "psfc"...'
        say
        say ' (x1o,x2o) = ('x1o','x2o')  -->  (x1,x2) = ('x1','x2')'
        say ' (y1o,y2o) = ('y1o','y2o')  -->  (y1,y2) = ('y1','y2')'
*
*  find minimum "psfc" within box area (x1,x2), (y1,y2)...
*
        'set x 'x1
        'set y 'y1
        'set z 1'
        'set t 'tt
        nxp = x2-x1+1
*        say '    nxp = 'nxp
        nyp = y2-y1+1
*        say '    nyp = 'nyp
        'set gxout print'
        'd minloc(min(psfc,y='y1',y='y2'),x='x1',x='x2')'
          rec=sublin(result,nxp+3)
          xc=subwrd(rec,1)
        'd minloc(min(psfc,x='x1',x='x2'),y='y1',y='y2')'
          rec=sublin(result,nyp+3)
          yc=subwrd(rec,1)
        say
        say 'location and value of "'track_type.n'"...'
        say
        say '    (xc,yc) = ('xc','yc')'
*
*  find "world" coordinates of (xc,yc) and convert "world" coordinates
*  to "xy" coordinates for plotting track of cyclone center...
*
        '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')'
        'd psfc' ; rec2 = sublin(result,2)
        tval.tt = subwrd(rec2,1)
        say '    'track_type.n' = 'tval.tt
        pause()
        'undefine psfc'
        tt=tt+1
      endwhile
    endif
    if( track_type.n = 'rel vorticity max' )
      say
      say '--- 'track_type.n ' ---'
      say
      ts=1
      tt=ts
      while ( tt <= nt )
        'set t 'tt
        say
        say '---------'
        say ' tt = 'tt
        say '---------'
        say
        'set lon 'lon1' 'lon2
        'set lat 'lat1' 'lat2
        'set z 1'
        'define rvor = hcurl(uur,vvr)'
*
*  display surface vorticity variable for interactive identification
*  of cyclone center... (substitute name of your surface vorticity variable; a
*  variable that is defined over land and ocean...)
*
        'set gxout shaded'
        'd rvor'
*
*  following interactive "box location" code adapted from "cbar_line_box.gs"...
*
        say 'Click where you want the left upper corner of the box'
        'query bpos'
        xb1o = subwrd(result,3)
        yb2o = subwrd(result,4)
        say 'Click where you want the right lower corner of the box'
        'query bpos'
        xb2o = subwrd(result,3)
        yb1o = subwrd(result,4)

        say
        say '...left upper corner of the box at X Y: 'xb1o' 'yb2o
        say '...right lower corner of the box at X Y: 'xb2o' 'yb1o
*
*  draw box to show specified area for defining cyclone center...
*
        'set line 1 1'
        'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
*
*  convert box "xy" coordinates to "grid" coordinates
*  for finding cyclone center...
*
        'q xy2gr 'xb1o' 'yb2o
        x1o = subwrd(result,3) ; x1 = math_int(x1o+0.5)
        y2o = subwrd(result,6) ; y2 = math_int(y2o+0.5)
        'q xy2gr 'xb2o' 'yb1o
        x2o = subwrd(result,3) ; x2 = math_int(x2o+0.5)
        y1o = subwrd(result,6) ; y1 = math_int(y1o+0.5)
        say
        say 'Specified grid coordinates (x1o,x2o), (y1o,y2o) from box,'
        say 'rounded to (x1,x2), (y1,y2) for finding maximum "rvor"...'
        say
        say ' (x1o,x2o) = ('x1o','x2o')  -->  (x1,x2) = ('x1','x2')'
        say ' (y1o,y2o) = ('y1o','y2o')  -->  (y1,y2) = ('y1','y2')'
*
*  find maximum "rvor" within box area (x1,x2), (y1,y2)...
*
        'set x 'x1
        'set y 'y1
        'set z 1'
        'set t 'tt
        nxp = x2-x1+1
*        say '    nxp = 'nxp
        nyp = y2-y1+1
*        say '    nyp = 'nyp
        '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
        say 'location and value of "'track_type.n'"...'
        say
        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')'
        'd rvor' ; rec2 = sublin(result,2)
        tval.tt = subwrd(rec2,1)
        say '    'track_type.n' = 'tval.tt
        pause()
        'undefine rvor'
        tt=tt+1
      endwhile
    endif
    if( track_type.n = 'wind speed max' )
      say
      say '--- 'track_type.n ' ---'
      say
      ts=1
      tt=ts
      while ( tt <= nt )
        'set t 'tt
        say
        say '---------'
        say ' tt = 'tt
        say '---------'
        say
        'set lon 'lon1' 'lon2
        'set lat 'lat1' 'lat2
        'set z 1'
        'define wspd = mag(uur,vvr)'
*
*  display surface wind speed variable for interactive identification
*  of cyclone center... (substitute name of your surface wind speed variable; a
*  variable that is defined over land and ocean...)
*
        'set gxout shaded'
        'd wspd'
*
*  following interactive "box location" code adapted from "cbar_line_box.gs"...
*
        say 'Click where you want the left upper corner of the box'
        'query bpos'
        xb1o = subwrd(result,3)
        yb2o = subwrd(result,4)
        say 'Click where you want the right lower corner of the box'
        'query bpos'
        xb2o = subwrd(result,3)
        yb1o = subwrd(result,4)

        say
        say '...left upper corner of the box at X Y: 'xb1o' 'yb2o
        say '...right lower corner of the box at X Y: 'xb2o' 'yb1o
*
*  draw box to show specified area for defining cyclone center...
*
        'set line 1 1'
        'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
*
*  convert box "xy" coordinates to "grid" coordinates
*  for finding cyclone center...
*
        'q xy2gr 'xb1o' 'yb2o
        x1o = subwrd(result,3) ; x1 = math_int(x1o+0.5)
        y2o = subwrd(result,6) ; y2 = math_int(y2o+0.5)
        'q xy2gr 'xb2o' 'yb1o
        x2o = subwrd(result,3) ; x2 = math_int(x2o+0.5)
        y1o = subwrd(result,6) ; y1 = math_int(y1o+0.5)
        say
        say 'Specified grid coordinates (x1o,x2o), (y1o,y2o) from box,'
        say 'rounded to (x1,x2), (y1,y2) for finding maximum "wspd"...'
        say
        say ' (x1o,x2o) = ('x1o','x2o')  -->  (x1,x2) = ('x1','x2')'
        say ' (y1o,y2o) = ('y1o','y2o')  -->  (y1,y2) = ('y1','y2')'
*
*  find maximum "wspd" within box area (x1,x2), (y1,y2)...
*
        'set x 'x1
        'set y 'y1
        'set z 1'
        'set t 'tt
        nxp = x2-x1+1
*        say '    nxp = 'nxp
        nyp = y2-y1+1
*        say '    nyp = 'nyp
        '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
        say 'location and value of "'track_type.n'"...'
        say
        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')'
        'd wspd' ; rec2 = sublin(result,2)
        tval.tt = subwrd(rec2,1)
        say '    'track_type.n' = 'tval.tt
        pause()
        'undefine wspd'
        tt=tt+1
      endwhile
    endif

*  Clear frame and re-plot base map...

    plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2,model)

*  Plotting track...

*  Specify drawing primitives

    marktype = mark.n
    lcolor = color.n
    lstyle = line.n
    lformat = label_format.n
    lscale = label_scale.n

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

    tt=ts
    xc = xpos.tt
    yc = ypos.tt
    'set line 'lcolor' 'lstyle' 'lthick
    'draw mark 'marktype' 'xc' 'yc' 'marksize
    'set string 1 l 4'
    'set strsiz 0.1'
    say 'Point: 'xc' 'yc
    say 'Click where you want the left edge of the track value label'
    'query bpos'
    xpl = subwrd(result,3)
    ypl = subwrd(result,4)
    'draw string 'xpl' 'ypl' 'math_format(lformat,lscale*tval.tt)
    tt = tt+1
    while ( tt <= nt )
      plot_label = true
      xn = xpos.tt
      yn = ypos.tt
      'draw line 'xc' 'yc' 'xn' 'yn
      'draw mark 'marktype' 'xn' 'yn' 'marksize
      say
      say 'Point: 'xn' 'yn
      if( xn = xc & yn = yc )
        say '*** overlapping point ***'
        say 'Do you want to plot a(nother) label at this point? y or n'
        pull answer
        if( answer = 'n' )
          plot_label = false
        endif
      else
        say 'Do you want to plot a label at this point? y or n'
        pull answer
        if( answer = 'n' )
          plot_label = false
        endif
      endif
      if( plot_label = true )
        'set string 1 l 4'
        'set strsiz 0.1'
        say 'Click where you want the left edge of the label'
        'query bpos'
        xpl = subwrd(result,3)
        ypl = subwrd(result,4)
        'draw string 'xpl' 'ypl' 'math_format(lformat,lscale*tval.tt)
      else
        say '...skipping label'
      endif
      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

function plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2,model)

*  Plots base map...

  'c'
  'set vpage 'plot_vpage
  'set parea 'plot_area
  'set lon 'lon1' 'lon2
  'set lat 'lat1' 'lat2
  'set mproj latlon'
  'set mpt * 1 1 5'
  'set mpdset hires'
  'set poli on'
*  'set grid on'
  'set grads off'
  'set gxout contour'
  '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 \\ Cyclone Track ("'model'")'
  'draw xlab Longitude \\'
  'draw ylab \\ Latitude'

return
