[gradsusr] how to find typhoon center via minloc function

sushant puranik sushantpuranik at gmail.com
Sat Apr 30 07:25:09 EDT 2011


Hello Charles
I tried maxwind.gs it works fine as per your suggestion. But i want to know
at which lat lon maximum value is instead of x and y location. So can you
please tell me how to convert x y in lat lon.

Thanks


Sushant

On Sat, Apr 30, 2011 at 11:43 AM, sushant puranik
<sushantpuranik at gmail.com>wrote:

> Thank you for help i will try it. Please don't mind if error occurs i'll
> mail you.
> Thanks again
>
>
> On Sat, Apr 30, 2011 at 1:32 AM, Charles Seman <Charles.Seman at noaa.gov>wrote:
>
>> Hello Sushant,
>>
>> I think it should be possible.  Please find attached an email exchange to
>> this List from several years ago "hurricane_tracking.txt" in which Kun-Hsuan
>> Chou gave a method to track the maximum wind speed using maxloc and max.
>>  This technique works on a 2D field.  To find the location of the max wind
>> in 4D wind dataset, a short script "maxwind.gs" was coded that blends the
>> technique used in a local script "plot_hurricane_center3.gs" and the
>> technique in Kun-Hsuan Chou's "hurricane_tracking.txt". ("
>> plot_hurricane_center3.gs" itself was developed using Joe Covert's script
>> "plot_tc_shi.gs" and Kun-Hsuan Chou's "hurricane_tracking.txt".)  "
>> maxwind.gs" loops over the time levels and finds the value of the max
>> wind and its (x,y,z) location for each time level.  These values are printed
>> out to the screen.  Here is a copy of what the "maxwind.gs" script
>> printed out using GrADS 1.9b4 for a local model "ucomp" climatology dataset:
>> ---
>> ...
>> (nx,ny,nz,nt) = (144,90,23,12)
>> tt = 1
>>  max ucomp = 81.833
>>  xc = 137
>>  yc = 71
>>  zc = 23
>>  ucomp(x=137,y=71,z=23) = 81.833
>> tt = 2
>>  max ucomp = 68.9826
>>  xc = 57
>>  yc = 61
>>  zc = 10
>>  ucomp(x=57,y=61,z=10) = 68.9826
>> tt = 3
>>  max ucomp = 62.7262
>>  xc = 56
>>  yc = 60
>>  zc = 11
>>  ucomp(x=56,y=60,z=11) = 62.7262
>> tt = 4
>>  max ucomp = 61.6225
>>  xc = 10
>>  yc = 15
>>  zc = 23
>>  ucomp(x=10,y=15,z=23) = 61.6225
>> tt = 5
>>  max ucomp = 85.5501
>>  xc = 4
>>  yc = 18
>>  zc = 23
>>  ucomp(x=4,y=18,z=23) = 85.5501
>> tt = 6
>>  max ucomp = 109.233
>>  xc = 4
>>  yc = 19
>>  zc = 23
>>  ucomp(x=4,y=19,z=23) = 109.233
>> tt = 7
>>  max ucomp = 119.698
>>  xc = 143
>>  yc = 17
>>  zc = 23
>>  ucomp(x=143,y=17,z=23) = 119.698
>> tt = 8
>>  max ucomp = 117.018
>>  xc = 2
>>  yc = 16
>>  zc = 22
>>  ucomp(x=2,y=16,z=22) = 117.018
>> tt = 9
>>  max ucomp = 98.7713
>>  xc = 63
>>  yc = 9
>>  zc = 20
>>  ucomp(x=63,y=9,z=20) = 98.7713
>> tt = 10
>>  max ucomp = 74.1137
>>  xc = 65
>>  yc = 7
>>  zc = 19
>>  ucomp(x=65,y=7,z=19) = 74.1137
>> tt = 11
>>  max ucomp = 77.7408
>>  xc = 137
>>  yc = 67
>>  zc = 23
>>  ucomp(x=137,y=67,z=23) = 77.7408
>> tt = 12
>>  max ucomp = 84.3092
>>  xc = 16
>>  yc = 70
>>  zc = 23
>>  ucomp(x=16,y=70,z=23) = 84.3092
>> ...finished.
>> ---
>>
>> It seems the script is working OK.  You can give it a try for one of your
>> datasets.  To run the script as it is coded, run it in batch mode as
>> follows:
>> % grads -blcx "maxwind.gs nc_file wind"
>> or on the GrADS command line
>> http://grads.iges.org/grads/gadoc/gradcomdrun.html as (no quotes as with
>> the batch command):
>> ga-> run maxwind.gs nc_file wind
>> where "nc_file" is the netCDF file containing the "wind" variable name
>> (the "nc_file" and "wind" are generic names; please replace them with the
>> actual netCDF file and wind variable in the netCDF file).  Note, please also
>> include the quotes surrounding "maxwind.gs nc_file wind" for the batch
>> command or else the batch command won't work.
>>
>> I hope this helps.  Please let me know if you have any questions.
>>
>> Thanks,
>> Chuck
>>
>> sushant puranik wrote:
>>
>>> Hello Charles
>>> I want to make a jet stream plot. Can i make it with the help of
>>> hurricane center script by changing minloc and min to maxloc and max.
>>> any suggestion?
>>>
>>> thanks
>>>
>>> Sushant
>>> On Wed, Jan 5, 2011 at 1:01 AM, Charles Seman <Charles.Seman at noaa.gov<mailto:
>>> Charles.Seman at noaa.gov>> wrote:
>>>
>>>    Jie,
>>>
>>>    Attached is a file File_Error_1_Grads_1.9-Plot_tc_shi.gs.txt
>>>    containing the text from an email exchange to the GrADS Listserv
>>>    between Joe Covert and Diane Stokes on 9/13/2004 -- on 4/27/06,
>>>    found via a Google search for "plot_tc_shi.gs
>>>    <http://plot_tc_shi.gs>":
>>> http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs
>>>
>>>
>>>    The attached GrADS script plot_hurricane_center3.gs
>>>    <http://plot_hurricane_center3.gs> (along with pause.gsf, aGrADS
>>>
>>>    script function used by the plot script) is a locally modified
>>>    version derived from Joe Covert's script "plot_tc_shi.gs
>>>    <http://plot_tc_shi.gs>" (Joe.Covert at noaa.gov
>>>    <mailto:Joe.Covert at noaa.gov>) which was posted to the GrADS Listserv
>>>
>>>    on 9/13/2004 (his script was named plot_tc_shi.gs
>>>    <http://plot_tc_shi.gs>) and code from hurricane_tracking.txt (taken
>>>
>>>    from an email exchange posted to the GrADS Listserv by Kun-Hsuan
>>>    Chou and Arturo Caracas Uribe in Oct, 2004)...
>>>
>>>    Please find below code from hurricane_tracking.txt (taken from an
>>>    email exchange posted to the GrADS Listserv by Kun-Hsuan Chou and
>>>    Arturo Caracas Uribe in Oct, 2004) illustrating a technique for
>>>    finding max wind in a 2D field... and below that some code from
>>>    attached script plot_hurricane_center3.gs
>>>    <http://plot_hurricane_center3.gs> to find the min sea-level
>>>
>>>    pressure in a 2D field...
>>>
>>>    code from Kun-Hsuan Chou to find maximum wind speed from
>>>    hurricane_tracking.txt:
>>>    ---
>>>    'd maxloc(max(mag(u,v),lon=120,lon=130),lat=15,lat=25)'
>>>     line=sublin(result,2)
>>>     ygrd=subwrd(line,4)
>>>    'd maxloc(max(mag(u,v),lat=15,lat=25),lon=120,lon=130)'
>>>     line=sublin(result,2)
>>>     xgrd=subwrd(line,4)
>>>    'set x 'xgrd
>>>    lonval = subwrd(result,4)
>>>    'set y 'ygrd
>>>    latval = subwrd(result,4)
>>>    'q w2xy 'lonval' 'latval
>>>    xpos = subwrd(result,3)
>>>    ypos = subwrd(result,6)
>>>    'draw mark 1 'xpos' 'ypos' .2'
>>>    ---
>>>
>>>    sample code from plot_hurricane_center3.gs
>>>    <http://plot_hurricane_center3.gs> to find hurricane center (see
>>>
>>>    script for supporting code):
>>>    ---
>>>    *
>>>    * find minimum "pressure" within box area (x1,x2), (y1,y2)...
>>>    *
>>>        'set x 'x1
>>>        'set y 'y1
>>>        'set z 1'
>>>        'set t 'tt
>>>        'set gxout print'
>>>        nxp = x2-x1+1
>>>    *      say '    nxp = 'nxp
>>>        nyp = y2-y1+1
>>>    *      say '    nyp = 'nyp
>>>        'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'
>>>          rec=sublin(result,nxp+3)
>>>          xc=subwrd(rec,1)
>>>        'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'
>>>          rec=sublin(result,nyp+3)
>>>          yc=subwrd(rec,1)
>>>        say
>>>        say 'location of minimum "pressure"...'
>>>        say
>>>        say ' (xc,yc) = ('xc','yc')'
>>>    *
>>>    * find "world" coordinates of (xc,yc) and convert "world" coordinates
>>>    * to "xy" coordinates for plotting track of hurricane 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 = subwrd(result,3)
>>>        ypos = subwrd(result,6)
>>>        say ' (xpos,ypos) = ('xpos','ypos')'
>>>    *
>>>    * write (xpos,ypos) to output file...
>>>    *
>>>        res = write (''hurricane.nf <http://hurricane.nf>'','t = 'tt'
>>>
>>>    xpos = 'xpos' ypos = 'ypos'')
>>>
>>>        pause()
>>>    ---
>>>
>>>    Please let me know if you have any questions.
>>>
>>>    Hope this helps,
>>>    Chuck
>>>
>>>    Jie TANG wrote:
>>>
>>>
>>>        hi,grads folks ,
>>>         I am using grads v1.9, trying to find the center of typhoon via
>>>        minloc function.
>>>        the key script is shown as below:
>>>        slplat=min(slp,lat=19,lat=26)
>>>  tclon=minloc(slplat,lon=119,lon=126)
>>>  slplon=min(slp,lon=119,lon=126)
>>> tclat=minloc(slplon,lat=19,lat=26)
>>>
>>>        but grads tell me that "function not found min " and when i
>>>        changed my script to be :
>>>        tclat=min(minloc(slpt,lon=119,lon=126), lat=19,lat=26)
>>>        tclon=min(minloc(slpt, lat=19,lat=26 )lon=119,lon=126)
>>>
>>>
>>>        how can i finx my scrpit ? thank you .:)
>>>        --
>>>        TANG Jie
>>>        Email: totangjie at gmail.com <mailto:totangjie at gmail.com>
>>>        <mailto:totangjie at gmail.com <mailto:totangjie at gmail.com>>
>>>
>>>
>>>        Tel: 0086-2154896104
>>>        Shanghai Typhoon Institute,China
>>>
>>>  ------------------------------------------------------------------------
>>>
>>>        _______________________________________________
>>>        gradsusr mailing list
>>>        gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
>>>
>>>        http://gradsusr.org/mailman/listinfo/gradsusr
>>>
>>>
>>>    --
>>>    Please note that Charles.Seman at noaa.gov
>>>    <mailto:Charles.Seman at noaa.gov> should be considered my NOAA
>>>    email address, not cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>.
>>>
>>>
>>>    ********************************************************************
>>>    Charles Seman                                Charles.Seman at noaa.gov
>>>    <mailto:Charles.Seman at noaa.gov>
>>>
>>>    U.S. Department of Commerce / NOAA / OAR
>>>    Geophysical Fluid Dynamics Laboratory         voice: (609) 452-6547
>>>    201 Forrestal Road                              fax: (609) 987-5063
>>>    Princeton, NJ  08540-6649            http://www.gfdl.noaa.gov/~cjs/
>>>    ********************************************************************
>>>
>>>    "The contents of this message are mine personally and do not reflect
>>> any
>>>    official or unofficial position of the United States Federal
>>> Government,
>>>    the United States Department of Commerce, or NOAA."
>>>
>>>
>>>
>>>    worked great ... thanks.  joe
>>>
>>>    Diane Stokes wrote:
>>>
>>>     > Joe,
>>>     >
>>>     > Read works fine for me in 1.9.
>>>     >
>>>     > I think this is the issue where 'pull' in 1.9 throws in a carriage
>>>     > return.  The read statement is not getting the exact filename.
>>>     >
>>>     > Try adding:
>>>     >   file=sublin(file,1)
>>>     > after each:
>>>     >   pull file
>>>     >
>>>     > Diane
>>>     >
>>>     > Joe Covert wrote:
>>>     >
>>>     >> There appears to be a problem with version 1.9 of GrADS "read"
>>>    function.
>>>     >>
>>>     >> The plot_tc_shi.gs <http://plot_tc_shi.gs> script (see below)
>>>
>>>    uses the read function to open the
>>>     >> track data file and the program gives the subject error message
>>>    "File
>>>     >> Error 1".
>>>     >>
>>>     >> jc
>>>     >>
>>>     >>
>>>     >> *  Script to draw an hurricane-track plot.
>>>     >> *  Does little error checking on the input file.
>>>     >> *  Assumes the input file is set up as follows:
>>>     >> *
>>>     >> *    Line 1:  Title
>>>     >> *    Line 2:  Drawing primitives for marks: marktype size
>>>     >> *    Line 3:  Drawing primitives for lines: color style thickness
>>>     >> *    Line 4:  Starting hour and the interval of plotting points
>>>     >> *             e.g., 0 6 means that track starts at 0 hour and mark
>>>     >> *                   will be plotted every 6 hours.
>>>     >> *    Rest of lines:  hour  long.  lat.
>>>     >> *             e.g.,   0    -70.5  25.0
>>>     >> *                     6    -71.8  25.2
>>>     >> *                            :
>>>     >> *                            :
>>>     >> *
>>>     >> *  Also assumes that a file has been opened (any file, doesn't
>>>     >> *  matter -- the set command doesn't work until a file has been
>>>     >> *  opened).
>>>     >> *
>>>     >>
>>>     >> function main()
>>>     >>
>>>     >> *  'clear'
>>>     >>
>>>     >>   'open dummy.ctl'
>>>     >>   'set lat 20 50'
>>>     >>   'set lon -90 -30'
>>>     >>   'set mpdset hires'
>>>     >>   'set poli on'
>>>     >>   'draw map'
>>>     >>
>>>     >>
>>>     >> say 'Enter File Name: (type q to stop)'
>>>     >> pull fname
>>>     >>
>>>     >> while (fname != 'q')
>>>     >>
>>>     >> *  Read the 1st record: Title
>>>     >>
>>>     >>   ret = read(fname)
>>>     >>   rc = sublin(ret,1)
>>>     >>   if (rc>0)
>>>     >>       say 'File Error 1'
>>>     >>       return
>>>     >>   endif
>>>     >>   title = sublin(ret,2)
>>>     >>   say title
>>>     >>
>>>     >> *  Read the drawing primitives
>>>     >>
>>>     >>   ret = read(fname)
>>>     >>   rc = sublin(ret,1)
>>>     >>   if (rc>0)
>>>     >>      say 'File Error 2'
>>>     >>      return
>>>     >>   endif
>>>     >>   dpline = sublin(ret,2)
>>>     >>   marktype = subwrd(dpline,1)
>>>     >>   marksize = subwrd(dpline,2)
>>>     >>   ret = read(fname)
>>>     >>   rc = sublin(ret,1)
>>>     >>   if (rc>0)
>>>     >>      say 'File Error 3'
>>>     >>      return
>>>     >>   endif
>>>     >>   dpline = sublin(ret,2)
>>>     >>   lcolor = subwrd(dpline,1)
>>>     >>   lstyle = subwrd(dpline,2)
>>>     >>   lthick = subwrd(dpline,3)
>>>     >>   say ' marktype, marksize, lcolor, lstyle and lthick:'
>>>     >>   say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick
>>>     >>
>>>     >> * Read starting hour and the interval hours of plotting points
>>>     >>
>>>     >>   ret = read(fname)
>>>     >>   rc = sublin(ret,1)
>>>     >>   if (rc>0)
>>>     >>      say 'File Error 4'
>>>     >>      return
>>>     >>   endif
>>>     >>   dhour = sublin(ret,2)
>>>     >>   start = subwrd(dhour,1)
>>>     >>   jump = subwrd(dhour,2)
>>>     >>   say ' starting hour and the interval hours of plotting points:'
>>>     >>   say '  'start' 'jump
>>>     >>
>>>     >> *  Read all data points
>>>     >>
>>>     >>   ret = read(fname)
>>>     >>   rc = sublin(ret,1)
>>>     >>   while (rc = 0)
>>>     >>       loc = sublin(ret,2)
>>>     >>       hour = subwrd(loc,1)
>>>     >>       dlong.hour = subwrd(loc,2)
>>>     >>       dlat.hour = subwrd(loc,3)
>>>     >> *      prompt ' hour ' hour ' are read,'
>>>     >> *      say ' long=' dlong.hour '    lat=' dlat.hour
>>>     >>       ret = read(fname)
>>>     >>       rc = sublin(ret,1)
>>>     >>   endwhile
>>>     >>
>>>     >>   if (rc!=2 & rc!=0)
>>>     >>          say 'File Error 5, rc=' rc
>>>     >>          return
>>>     >>   endif
>>>     >>
>>>     >>   endhour = hour
>>>     >>   say ' endhour=' endhour
>>>     >>
>>>     >> * Plotting
>>>     >>
>>>     >>   'set line 'lcolor' 'lstyle' 'lthick
>>>     >>   'query w2xy 'dlong.start' 'dlat.start
>>>     >>   xprev = subwrd(result,3)
>>>     >>   yprev = subwrd(result,6)
>>>     >>   'draw mark 'marktype' 'xprev' 'yprev' 'marksize
>>>     >>   next = start+jump
>>>     >>   while (next <= endhour)
>>>     >> *      say ' 'dlong.start' 'dlat.start
>>>     >>       'query w2xy 'dlong.next' 'dlat.next
>>>     >>       xnext = subwrd(result,3)
>>>     >>       ynext = subwrd(result,6)
>>>     >>       'draw line 'xprev' 'yprev' 'xnext' 'ynext
>>>     >> *      say ' 'xprev' 'yprev' 'xnext' 'ynext
>>>     >>       'draw mark 'marktype' 'xnext' 'ynext' 'marksize
>>>     >>       next = next+jump
>>>     >>       xprev = xnext
>>>     >>       yprev = ynext
>>>     >>   endwhile
>>>     >>
>>>     >>   say fname ' is working fine.'
>>>     >>
>>>     >> * read in the filename to be plotted next
>>>     >>
>>>     >>   say 'Enter File Name: (type q to stop)'
>>>     >>   pull fname
>>>     >>
>>>     >> endwhile
>>>     >>
>>>
>>>
>>>
>>>  *----------------------------------------------------------------------
>>>    *
>>>    *  Adapted from plot_hurricane_center2c.gs
>>>    <http://plot_hurricane_center2c.gs>, which was adapted from:
>>>
>>>    *
>>>    *  (1) "plot_hurricane_center2[a][b].gs" (from
>>>    "plot_hurricane_center.gs <http://plot_hurricane_center.gs>"),
>>>
>>>    *      locally modified versions derived from Joe Covert's script
>>>    *      "plot_tc_shi.gs <http://plot_tc_shi.gs>" (Joe.Covert at noaa.gov
>>>    <mailto:Joe.Covert at noaa.gov>), which was posted to the
>>>
>>>    *      GrADS Listserv on 9/13/2004 (his script was named
>>>    plot_tc_shi.gs <http://plot_tc_shi.gs>)
>>>
>>>    *      on 4/27/06, found via a Google search for "plot_tc_shi.gs
>>>    <http://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...
>>>
>>>  ************************************************************************
>>>
>>>    **  nc_file1 = ''
>>>    **  nc_file2 = ''
>>>    **  nc_file3 = ''
>>>    **  nc_file = ''nc_file1' 'nc_file2' 'nc_file3''
>>>     nc_file = 'your_netCDF_file'
>>>
>>>     lon1 = 298 ; lon2 = 308
>>>     lat1 =  20 ; lat2 =  25
>>>
>>>
>>>  ************************************************************************
>>>    *  Legend captions, colors, and marks for the plots
>>>
>>>  ************************************************************************
>>>
>>>    **  nmodels = 1
>>>    **  model.1 = 'isabel_ras'
>>>    **  model.2 = 'isabel_isotke'
>>>    **  model.3 = 'isabel_my25'
>>>
>>>    **  titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2
>>>    **  titles.2 = model.2 ; colors.2 = 3 ; marks.2 = 2
>>>    **  titles.3 = model.3 ; colors.3 = 4 ; marks.3 = 2
>>>    ***
>>>    ***  Define the legend plotter input...
>>>    ***
>>>    **  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
>>>    **  n=1
>>>    **  while ( n <= nmodels )
>>>    **    titles = ''titles' "'titles.n'"'
>>>    **    colors = ''colors' 'colors.n
>>>    **    marks = ''marks' "'marks.n'"'
>>>    **    lines = ''lines' 1'
>>>    **    n=n+1
>>>    **  endwhile
>>>
>>>    **  legend_info = ''colors' 'lines' 'marks' 'titles' -p'
>>>
>>>     nmodels = 1
>>>     model.1 = 'title_name_for_your_netCDF_dataset'
>>>
>>>     titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2
>>>    *
>>>    *  Define the legend plotter input...
>>>    *
>>>     titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
>>>     n=1
>>>     while ( n <= nmodels )
>>>       titles = ''titles' "'titles.n'"'
>>>       colors = ''colors' 'colors.n
>>>       marks = ''marks' "'marks.n'"'
>>>       lines = ''lines' 1'
>>>       n=n+1
>>>     endwhile
>>>
>>>     legend_info = ''colors' 'lines' 'marks' 'titles' -p'
>>>
>>>    *  ...other drawing primitives
>>>
>>>     marksize = 0.1
>>>
>>>     lstyle = 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.gx'
>>>
>>>
>>>  *-----------------------------------------------------------------------
>>>    *  Define and plot track of minimum surface pressure for each
>>> dataset...
>>>
>>>  *-----------------------------------------------------------------------
>>>
>>>     nf=1
>>>     while ( nf <= nmodels )
>>>
>>>    * Open dataset file and define parameters...
>>>
>>>       fname = subwrd(nc_file,nf)
>>>       'sdfopen 'fname
>>>       'set dfile 'nf
>>>       '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)
>>>
>>>       if( nf = 1 )
>>>
>>>    * Plot base map...
>>>
>>>         plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>>>
>>>       endif
>>>
>>>    * Define locations of "hurricane center" for plotting points, and
>>> write
>>>    * locations out to ASCII text file for read-in for plotting track
>>>    later...
>>>
>>>       hurricane.nf <http://hurricane.nf> = 'hurricane_'nf'.txt'
>>>
>>>
>>>       ts=2
>>>       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'
>>>    *
>>>    * here is where you put in the name of your sea-level pressure
>>> variable:
>>>    *
>>>         'define psl = name_of_your_sea_level_pressure_variable'
>>>         'set gxout shaded'
>>>         'd psl'
>>>    *
>>>    * following interactive "box location" code adapted from
>>>    "cbar_line_box.gs <http://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 hurricane center...
>>>    *
>>>         'set line 1 1'
>>>         'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
>>>    *
>>>    * convert box "xy" coordinates to "grid" coordinates
>>>    * for finding minimum "pressure" (hurricane 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
>>> "pressure"...'
>>>         say
>>>         say ' (x1o,x2o) = ('x1o','x2o')  -->  (x1,x2) = ('x1','x2')'
>>>         say ' (y1o,y2o) = ('y1o','y2o')  -->  (y1,y2) = ('y1','y2')'
>>>    *
>>>    * find minimum "pressure" within box area (x1,x2), (y1,y2)...
>>>    *
>>>         'set x 'x1
>>>         'set y 'y1
>>>         'set z 1'
>>>         'set t 'tt
>>>         'set gxout print'
>>>         nxp = x2-x1+1
>>>    *      say '    nxp = 'nxp
>>>         nyp = y2-y1+1
>>>    *      say '    nyp = 'nyp
>>>         'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'
>>>           rec=sublin(result,nxp+3)
>>>           xc=subwrd(rec,1)
>>>         'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'
>>>           rec=sublin(result,nyp+3)
>>>           yc=subwrd(rec,1)
>>>         say
>>>         say 'location of minimum "pressure"...'
>>>         say
>>>         say ' (xc,yc) = ('xc','yc')'
>>>    *
>>>    * find "world" coordinates of (xc,yc) and convert "world" coordinates
>>>    * to "xy" coordinates for plotting track of hurricane 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 = subwrd(result,3)
>>>         ypos = subwrd(result,6)
>>>         say ' (xpos,ypos) = ('xpos','ypos')'
>>>    *
>>>    * write (xpos,ypos) to output file...
>>>    *
>>>         res = write (''hurricane.nf <http://hurricane.nf>'','t = 'tt'
>>>
>>>    xpos = 'xpos' ypos = 'ypos'')
>>>
>>>         pause()
>>>    *
>>>    * erase box showing area for defining hurricane center...
>>>    *
>>>    *      'set line 0 1'
>>>    *      'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
>>>         tt=tt+1
>>>       endwhile
>>>
>>>       res = close (''hurricane.nf <http://hurricane.nf>'')
>>>
>>>
>>>       nf=nf+1
>>>     endwhile
>>>
>>>    *  Clear frame and re-plot base map...
>>>
>>>     plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>>>
>>>    *  Plot track(s)...
>>>
>>>     nf=1
>>>     while ( nf <= nmodels )
>>>
>>>    *  Specify drawing primitives
>>>
>>>       marktype = marks.nf <http://marks.nf>
>>>       lcolor = colors.nf <http://colors.nf>
>>>
>>>
>>>       say ' marktype, marksize, lcolor, lstyle and lthick:'
>>>       say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick
>>>
>>>       tt=ts
>>>    *
>>>    *  read (xpos,ypos) from ASCII text input file...
>>>    *
>>>       res = read(''hurricane.nf <http://hurricane.nf>'') ; rc =
>>>
>>>    sublin(res,1)
>>>       if( rc != 0 )
>>>         say '"read" status = 'rc
>>>         exit
>>>       endif
>>>       rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec
>>>       xc = subwrd(rec,6) ; yc = subwrd(rec,9)
>>>
>>>       'set line 'lcolor' 'lstyle' 'lthick
>>>       'draw mark 'marktype' 'xc' 'yc' 'marksize
>>>       tt = tt+1
>>>       while ( tt <= nt )
>>>         res = read(''hurricane.nf <http://hurricane.nf>'') ; rc =
>>>
>>>    sublin(res,1)
>>>         if( rc != 0 )
>>>           say '"read" status = 'rc
>>>           exit
>>>         endif
>>>         rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec
>>>         xn = subwrd(rec,6) ; yn = subwrd(rec,9)
>>>         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
>>>
>>>       nf=nf+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)
>>>
>>>    *  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 \\ Track of Minimum Surface Pressure'
>>>     'draw xlab Longitude \\'
>>>     'draw ylab \\ Latitude'
>>>
>>>    return
>>>
>>>    _______________________________________________
>>>    gradsusr mailing list
>>>    gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
>>>
>>>    http://gradsusr.org/mailman/listinfo/gradsusr
>>>
>>>
>>>
>>>
>>>
>>>
>>> ------------------------------------------------------------------------
>>>
>>> _______________________________________________
>>> gradsusr mailing list
>>> gradsusr at gradsusr.org
>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>
>>
>> --
>>
>> Please note that Charles.Seman at noaa.gov should be considered my NOAA
>> email address, not cjs at gfdl.noaa.gov.
>>
>> ********************************************************************
>>  Charles Seman                                Charles.Seman at noaa.gov
>>  U.S. Department of Commerce / NOAA / OAR
>>  Geophysical Fluid Dynamics Laboratory         voice: (609) 452-6547
>>  201 Forrestal Road                              fax: (609) 987-5063
>>  Princeton, NJ  08540-6649            http://www.gfdl.noaa.gov/~cjs/
>> ********************************************************************
>>
>> "The contents of this message are mine personally and do not reflect any
>> official or unofficial position of the United States Federal Government,
>> the United States Department of Commerce, or NOAA."
>>
>> *
>> *
>> * Subject: Re: Tracking Hurricanes
>> * From: Kun-Hsuan Chou <cwujou at TYPHOON.AS.NTU.EDU.TW>
>> * Date: Fri, 8 Oct 2004 08:03:20 +0800
>> * To: GRADSUSR at LIST.CINECA.IT
>> *
>> * Hi,
>> *
>> * I write a sample to track the maxminum wind speed, it should be the
>> script that you want.
>>
>> *
>> 'd maxloc(max(mag(u,v),lon=120,lon=130),lat=15,lat=25)'
>>  line=sublin(result,2)
>>  ygrd=subwrd(line,4)
>> 'd maxloc(max(mag(u,v),lat=15,lat=25),lon=120,lon=130)'
>>  line=sublin(result,2)
>>  xgrd=subwrd(line,4)
>> 'set x 'xgrd
>> lonval = subwrd(result,4)
>> 'set y 'ygrd
>> latval = subwrd(result,4)
>> 'q w2xy 'lonval' 'latval
>> xpos = subwrd(result,3)
>> ypos = subwrd(result,6)
>> 'draw mark 1 'xpos' 'ypos' .2'
>>
>>
>> *============================================================================
>> *Kun-Hsuan Chou
>> *Postdoctoral Researcher
>> *Department of Atmospheric Sciences
>> *National Taiwan University
>> *No. 1, Sec. 4, Roosevelt Rd., Taipei 106, Taiwan.
>> *Phone:   886-2-23625896 ext. 243
>> *Mobile:  0920079247
>> *E-mail:   cwujou at typhoon.as.ntu.edu.tw
>> *----- Original Message -----
>> *From: Arturo Caracas Uribe
>> *To: GRADSUSR at LIST.CINECA.IT
>> *Sent: Thursday, October 07, 2004 11:48 PM
>> *Subject: Tracking Hurricanes
>> *
>> *Does anyone know how to get the maximum values of vorticity of a domain?
>> (the value, lat and lon).I need the positions in order to track a hurricane.
>> *Thanks in advance
>> *
>> *Arturo Caracas
>> *Atmospheric Sciences, México.
>>
>>
>> *----------------------------------------------------------------------
>> *
>> *  Adapted from plot_hurricane_center2c.gs, which was adapted from:
>> *
>> *  (1) "plot_hurricane_center2[a][b].gs" (from "plot_hurricane_center.gs
>> "),
>>
>> *      locally modified versions derived from Joe Covert's script
>> *      "plot_tc_shi.gs" (Joe.Covert at 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...
>> ************************************************************************
>>
>> **  nc_file1 = ''
>> **  nc_file2 = ''
>> **  nc_file3 = ''
>> **  nc_file = ''nc_file1' 'nc_file2' 'nc_file3''
>>  nc_file = 'your_netCDF_file'
>>
>>  lon1 = 298 ; lon2 = 308
>>  lat1 =  20 ; lat2 =  25
>>
>> ************************************************************************
>> *  Legend captions, colors, and marks for the plots
>> ************************************************************************
>>
>> **  nmodels = 1
>> **  model.1 = 'isabel_ras'
>> **  model.2 = 'isabel_isotke'
>> **  model.3 = 'isabel_my25'
>>
>> **  titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2
>> **  titles.2 = model.2 ; colors.2 = 3 ; marks.2 = 2
>> **  titles.3 = model.3 ; colors.3 = 4 ; marks.3 = 2
>> ***
>> ***  Define the legend plotter input...
>> ***
>> **  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
>> **  n=1
>> **  while ( n <= nmodels )
>> **    titles = ''titles' "'titles.n'"'
>> **    colors = ''colors' 'colors.n
>> **    marks = ''marks' "'marks.n'"'
>> **    lines = ''lines' 1'
>> **    n=n+1
>> **  endwhile
>>
>> **  legend_info = ''colors' 'lines' 'marks' 'titles' -p'
>>
>>  nmodels = 1
>>  model.1 = 'title_name_for_your_netCDF_dataset'
>>
>>  titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2
>> *
>> *  Define the legend plotter input...
>> *
>>  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
>>  n=1
>>  while ( n <= nmodels )
>>    titles = ''titles' "'titles.n'"'
>>    colors = ''colors' 'colors.n
>>    marks = ''marks' "'marks.n'"'
>>    lines = ''lines' 1'
>>    n=n+1
>>  endwhile
>>
>>  legend_info = ''colors' 'lines' 'marks' 'titles' -p'
>>
>> *  ...other drawing primitives
>>
>>  marksize = 0.1
>>
>>  lstyle = 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.gx'
>>
>> *-----------------------------------------------------------------------
>> *  Define and plot track of minimum surface pressure for each dataset...
>> *-----------------------------------------------------------------------
>>
>>  nf=1
>>  while ( nf <= nmodels )
>>
>> * Open dataset file and define parameters...
>>
>>    fname = subwrd(nc_file,nf)
>>    'sdfopen 'fname
>>    'set dfile 'nf
>>    '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)
>>
>>    if( nf = 1 )
>>
>> * Plot base map...
>>
>>      plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>>
>>    endif
>>
>> * Define locations of "hurricane center" for plotting points, and write
>> * locations out to ASCII text file for read-in for plotting track later...
>>
>>    hurricane.nf = 'hurricane_'nf'.txt'
>>
>>
>>    ts=2
>>    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'
>> *
>> * here is where you put in the name of your sea-level pressure variable:
>> *
>>      'define psl = name_of_your_sea_level_pressure_variable'
>>      'set gxout shaded'
>>      'd psl'
>> *
>> * 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 hurricane center...
>> *
>>      'set line 1 1'
>>      'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
>> *
>> * convert box "xy" coordinates to "grid" coordinates
>> * for finding minimum "pressure" (hurricane 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 "pressure"...'
>>      say
>>      say ' (x1o,x2o) = ('x1o','x2o')  -->  (x1,x2) = ('x1','x2')'
>>      say ' (y1o,y2o) = ('y1o','y2o')  -->  (y1,y2) = ('y1','y2')'
>> *
>> * find minimum "pressure" within box area (x1,x2), (y1,y2)...
>> *
>>      'set x 'x1
>>      'set y 'y1
>>      'set z 1'
>>      'set t 'tt
>>      'set gxout print'
>>      nxp = x2-x1+1
>> *      say '    nxp = 'nxp
>>      nyp = y2-y1+1
>> *      say '    nyp = 'nyp
>>      'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'
>>        rec=sublin(result,nxp+3)
>>        xc=subwrd(rec,1)
>>      'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'
>>        rec=sublin(result,nyp+3)
>>        yc=subwrd(rec,1)
>>      say
>>      say 'location of minimum "pressure"...'
>>      say
>>      say ' (xc,yc) = ('xc','yc')'
>> *
>> * find "world" coordinates of (xc,yc) and convert "world" coordinates
>> * to "xy" coordinates for plotting track of hurricane 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 = subwrd(result,3)
>>      ypos = subwrd(result,6)
>>      say ' (xpos,ypos) = ('xpos','ypos')'
>> *
>> * write (xpos,ypos) to output file...
>> *
>>      res = write (''hurricane.nf'','t = 'tt' xpos = 'xpos' ypos =
>> 'ypos'')
>>
>>
>>      pause()
>> *
>> * erase box showing area for defining hurricane center...
>> *
>> *      'set line 0 1'
>> *      'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
>>      tt=tt+1
>>    endwhile
>>
>>    res = close (''hurricane.nf'')
>>
>>
>>    nf=nf+1
>>  endwhile
>>
>> *  Clear frame and re-plot base map...
>>
>>  plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>>
>> *  Plot track(s)...
>>
>>  nf=1
>>  while ( nf <= nmodels )
>>
>> *  Specify drawing primitives
>>
>>    marktype = marks.nf
>>    lcolor = colors.nf
>>
>>
>>    say ' marktype, marksize, lcolor, lstyle and lthick:'
>>    say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick
>>
>>    tt=ts
>> *
>> *  read (xpos,ypos) from ASCII text input file...
>> *
>>    res = read(''hurricane.nf'') ; rc = sublin(res,1)
>>
>>    if( rc != 0 )
>>      say '"read" status = 'rc
>>      exit
>>    endif
>>    rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec
>>    xc = subwrd(rec,6) ; yc = subwrd(rec,9)
>>
>>    'set line 'lcolor' 'lstyle' 'lthick
>>    'draw mark 'marktype' 'xc' 'yc' 'marksize
>>    tt = tt+1
>>    while ( tt <= nt )
>>      res = read(''hurricane.nf'') ; rc = sublin(res,1)
>>
>>      if( rc != 0 )
>>        say '"read" status = 'rc
>>        exit
>>      endif
>>      rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec
>>      xn = subwrd(rec,6) ; yn = subwrd(rec,9)
>>      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
>>
>>    nf=nf+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)
>>
>> *  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 \\ Track of Minimum Surface Pressure'
>>  'draw xlab Longitude \\'
>>  'draw ylab \\ Latitude'
>>
>> return
>>
>> function maxwind(args)
>>  'reinit'
>>  nc_file = subwrd(args,1)
>>  wind = subwrd(args,2)
>>  say nc_file
>>  say wind
>>  'sdfopen 'nc_file
>>  'q file'
>>  say result
>>
>>  rec5 = sublin(result,5)
>>  nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt =
>> subwrd(rec5,12)
>>  say '(nx,ny,nz,nt) = ('nx','ny','nz','nt')'
>>  nxp = nx
>>  nyp = ny
>>  nzp = nz
>>  'set gxout print'
>>
>>  tt = 1
>>  while ( tt <= nt )
>>    'set x 1'
>>    'set y 1'
>>
>>    'set z 1'
>>    'set t 'tt
>>    say 'tt = 'tt
>>    'd max(max(max('wind',x=1,x='nx'),y=1,y='ny'),z=1,z='nz')'
>>      rec = sublin(result,(nyp+1)*nzp+3)
>>      max_wind = subwrd(rec,1)
>>      say ' max 'wind' = 'max_wind
>>    'd maxloc(max(max('wind',y=1,y='ny'),z=1,z='nz'),x=1,x='nx')'
>>      rec = sublin(result,(nzp+1)*nxp+3)
>>      xc = subwrd(rec,1)
>>      say ' xc = 'xc
>>    'd maxloc(max(max('wind',x=1,x='nx'),z=1,z='nz'),y=1,y='ny')'
>>      rec = sublin(result,(nzp+1)*nyp+3)
>>      yc = subwrd(rec,1)
>>      say ' yc = 'yc
>>    'd maxloc(max(max('wind',x=1,x='nx'),y=1,y='ny'),z=1,z='nz')'
>>      rec = sublin(result,(nyp+1)*nzp+3)
>>      zc = subwrd(rec,1)
>>      say ' zc = 'zc
>>    'd 'wind'(x='xc',y='yc',z='zc')'
>>      rec = sublin(result,2)
>>      wind_xc_yc_zc = subwrd(rec,1)
>>      say ' 'wind'(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc
>>    tt=tt+1
>>  endwhile
>>  say '...finished.'
>>
>> return
>>
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org
>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>>
>
>
> --
>
>
>
>


-- 
सुशांत पुराणिक
वातावरणशास्त्र व अवकाशशास्त्र विभाग,
पुणे विद्यापीठ
पुणे-०७, महाराष्ट्र
भारत.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20110430/80d5c6fb/attachment-0001.html 


More information about the gradsusr mailing list