[gradsusr] how to find typhoon center via minloc function
Charles Seman
Charles.Seman at noaa.gov
Wed May 4 13:53:43 EDT 2011
Hi Sushant,
Sorry for the error in the "Run commands" listing I sent you
yesterday... Today, I downloaded your files into a subdirectory, and
then did the set-up for the GrADS command and found there was a
difference in the variable name in the GrADS command I sent you
("wind_nc_file") vs that given in the "set wind_file = (...)" command I
sent you... I don't have a record of the command history, but apparently
the "wind_nc_file" had been set earlier in the session, and I made a
mistake when re-issuing the commands for verifying and for sending them
to you... Sorry about that!
Try with these commands:
% set wind_nc_file = ( uwind05.nc vwind05.nc )
% set wind_var = ( uwnd vwnd )
% set wind_units = m/s
% set lev_units = hPa
% grads -pcx "maxwind2.gs $wind_nc_file[1] $wind_var[1] $wind_nc_file[2]
$wind_var[2] $wind_units $lev_units"
Note, it seemed that the plot looked better in "portrait" mode ("-p"
option above). Also, the parameters were changed to search the full
"x", "y", and "t" domains, and to limit the search in the "z" domain to
between 1000 hPa and 500 hPa (using output from an "\ncdump -c" command
to see what levels to include):
x1 = 1 ; x2 = nx
y1 = 1 ; y2 = ny
z1 = 1 ; z2 = 6
t1 = 1 ; t2 = nt
The revised script "maxwind2.gs" and a plot file maxwind2.pdf file are
attached. According to the plot, all of the max wind values were found
at 500 hPa... If you wish to limit the search to between 1000 hPa and
600 hPa, try "z1 = 1 ; z2 = 5". Also, If you are interested in the
low-level jet nearer to India, you can change the x1, x2, y1, and y2
parameters to change the horizontal domain for the search (it seems that
the max wind values found at 500 hPa are associated with the S
Hemisphere jet stream where it is going into their winter for the time
period of the datasets).
Please let me know if you find any problems and/or have any questions.
Thanks,
Chuck
sushant puranik wrote:
> Hello Charles
> Thank you for the script. Script for individual data file runs fine. But
> I am not able to run maxwind2.gs <http://maxwind2.gs> using two separate
> data files uwind05.nc <http://uwind05.nc> and vwind05.nc
> <http://vwind05.nc> (which are attached here with this mail).
>
> One thing can we check maximum wind at low level. I want to check
> maximum wind from 1000 upto 600 or 500hPa only. My study is related with
> low level jet for that purpose i need maximum wind at low level only.
>
>
> Thank you again
>
> On Wed, May 4, 2011 at 1:17 AM, Charles Seman <Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>> wrote:
>
> Sushant,
>
> Glad to hear the script was helpful.
>
> An updated script "maxwind.gs <http://maxwind.gs>" (attached) for a
> single wind file and wind variable (e.g. ucomp or vcomp) is now able
> to search over a user-specified domain.
>
> Another script "maxwind2.gs <http://maxwind2.gs>" (attached) has
> been coded based on "maxwind.gs <http://maxwind.gs>" to find the max
> wind magnitude of 2 wind variables (e.g. ucomp, vcomp) over a
> user-specified domain (although the wind magnitude is defined over
> all (x,y,z) points at each time level).
>
> The scripts were run over a subset of the total (x,y,z,t) domain to
> verify they work for a limited domain (note, you can set these in
> the scripts yourself):
> x1 = 2 ; x2 = nx-1
> y1 = 2 ; y2 = ny-1
> z1 = 2 ; z2 = nz-1
> t1 = 2 ; t2 = nt-1
>
> Run commands:
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set wind_file = (
> ucomp.1981-2000.climatology.nc
> <http://ucomp.1981-2000.climatology.nc>
> vcomp.1981-2000.climatology.nc <http://vcomp.1981-2000.climatology.nc> )
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set wind_var = ( ucomp vcomp )
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set wind_units = m/s
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set lev_units = hPa
> cjs: /home/cjs/grads/Sushant_Puranik/ --> grads -lcx "maxwind.gs
> <http://maxwind.gs> $wind_nc_file[1] $wind_var[1] $wind_units
> $lev_units"
> ...
> cjs: /home/cjs/grads/Sushant_Puranik/ --> grads -lcx "maxwind.gs
> <http://maxwind.gs> $wind_nc_file[2] $wind_var[2] $wind_units
> $lev_units"
> ...
> cjs: /home/cjs/grads/Sushant_Puranik/ --> grads -lcx "maxwind2.gs
> <http://maxwind2.gs> $wind_nc_file[1] $wind_var[1] $wind_nc_file[2]
> $wind_var[2] $wind_units $lev_units"
> ...
>
> Please find attached sample plot files: maxwind_ucomp.pdf,
> maxwind_vcomp.pdf, maxwind2.pdf
>
> The scripts were tested and appear to be working OK. Please let me
> know if you find any problems and/or have any questions.
>
> Thanks,
> Chuck
>
> sushant puranik wrote:
>
> Thanks again It works fine. Since i am interested in low level
> wind so can you please tell me how to make it to check for low
> levels only(i have data at 17 pressure levels). Also i have
> uwind and vwind data in separate files, how to run script with
> two files?
>
>
> thank you
> On Mon, May 2, 2011 at 11:04 PM, Charles Seman
> <Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>>>
> wrote:
>
> Sushant,
>
> Following Kun-Hsuan Chou in the "hurricane_tracking.txt" file (or
> plot_hurricane_center3.gs <http://plot_hurricane_center3.gs>
> <http://plot_hurricane_center3.gs>),
> please find attached a revised version "maxwind.gs
> <http://maxwind.gs>
> <http://maxwind.gs>" which plots the wind and its max for
> each time
> level.
>
> To run:
>
> % grads -lcx "maxwind.gs <http://maxwind.gs>
> <http://maxwind.gs> nc_file wind wind_units"
>
> or
>
> ga-> run maxwind.gs <http://maxwind.gs> <http://maxwind.gs>
> nc_file wind wind_units
>
> Where a new argument "wind_units" has been added which is
> plotted in
> a title. Be sure there are no blank spaces in the
> "wind_units" so
> that GrADS gets the correct string to plot (e.g. m/s or knots,
> etc.). Note, the script as it is clears the screen before
> each time
> level (it plots the wind field itself at the level of maximum
> wind
> and an "x" at the location of the max wind). I'm not sure if
> this
> is what you want? If not, I would suggest adapting the code from
> maxwind.gs <http://maxwind.gs> <http://maxwind.gs> into a
> modified version of
> plot_hurricane_center3.gs <http://plot_hurricane_center3.gs>
> <http://plot_hurricane_center3.gs> (and
> rename it to something more appropriate like plot_maxwind.gs
> <http://plot_maxwind.gs>
> <http://plot_maxwind.gs> or whatever you choose). One thing
> I found
> when running the script several times on the GrADS command
> line is
> that at some point the script exits with a "Memory allocation
> error!" condition. I'm not sure why this happens. In an
> attempt to
> prevent this, 'c events' and 'close 1' statements were added in
> addition to the previous 'reinit' command, but this did not
> prevent
> the problem.
>
> I hope this helps,
> Chuck
>
> sushant puranik wrote:
>
> Hello Charles
> I tried maxwind.gs <http://maxwind.gs>
> <http://maxwind.gs> <http://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
> <mailto:sushantpuranik at gmail.com>
> <mailto:sushantpuranik at gmail.com <mailto:sushantpuranik at gmail.com>>
> <mailto:sushantpuranik at gmail.com
> <mailto:sushantpuranik at gmail.com>
> <mailto:sushantpuranik at gmail.com
> <mailto: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
> <mailto:Charles.Seman at noaa.gov> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov> <mailto:Charles.Seman at noaa.gov
> <mailto: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 <http://maxwind.gs>
> <http://maxwind.gs>
> <http://maxwind.gs>" was coded that blends the
> technique
> used in
>
> a local script "plot_hurricane_center3.gs
> <http://plot_hurricane_center3.gs>
> <http://plot_hurricane_center3.gs>
> <http://plot_hurricane_center3.gs>" and the
> technique in
>
> Kun-Hsuan Chou's "hurricane_tracking.txt".
> ("plot_hurricane_center3.gs
> <http://plot_hurricane_center3.gs>
> <http://plot_hurricane_center3.gs>
> <http://plot_hurricane_center3.gs>"
>
> itself was developed using Joe Covert's script
> "plot_tc_shi.gs <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>" and Kun-Hsuan Chou's
> "hurricane_tracking.txt".) "maxwind.gs
> <http://maxwind.gs>
> <http://maxwind.gs> <http://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 <http://maxwind.gs>
> <http://maxwind.gs> <http://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 <http://maxwind.gs>
> <http://maxwind.gs>
> <http://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 <http://maxwind.gs>
> <http://maxwind.gs>
> <http://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 <http://maxwind.gs>
> <http://maxwind.gs>
> <http://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>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto: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://plot_tc_shi.gs>
> <http://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>
> <http://plot_hurricane_center3.gs>
> <http://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>
> <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>"
> (Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>
> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>>
> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>>>
> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>
> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>>
> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>
> <mailto: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>
> <http://plot_tc_shi.gs>
> <http://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>
> <http://plot_hurricane_center3.gs>
> <http://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>
> <http://plot_hurricane_center3.gs>
> <http://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>
> <http://hurricane.nf> <http://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>>
> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>
> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>>> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>
> <mailto:totangjie at gmail.com <mailto:totangjie at gmail.com>>
> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>
> <mailto:totangjie at gmail.com <mailto:totangjie at gmail.com>>>>
> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>
> <mailto:totangjie at gmail.com <mailto:totangjie at gmail.com>>
> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>
> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>>> <mailto:totangjie at gmail.com
> <mailto:totangjie at gmail.com>
> <mailto:totangjie at gmail.com <mailto:totangjie at gmail.com>>
> <mailto: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>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>>>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>
> <mailto: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>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto: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>
> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>>
> <mailto:cjs at gfdl.noaa.gov
> <mailto:cjs at gfdl.noaa.gov> <mailto:cjs at gfdl.noaa.gov
> <mailto:cjs at gfdl.noaa.gov>>>
> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>
> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>>
>
> <mailto:cjs at gfdl.noaa.gov
> <mailto:cjs at gfdl.noaa.gov> <mailto:cjs at gfdl.noaa.gov
> <mailto:cjs at gfdl.noaa.gov>>>>.
>
>
>
> ********************************************************************
> Charles Seman
> Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto: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> <http://plot_tc_shi.gs>
> <http://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>
> <http://plot_hurricane_center2c.gs>
> <http://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>
> <http://plot_hurricane_center.gs>
> <http://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> <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>" (Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>
> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>>
> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>>>
> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>
> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>>
> <mailto:Joe.Covert at noaa.gov
> <mailto:Joe.Covert at noaa.gov>
> <mailto: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>
> <http://plot_tc_shi.gs>
> <http://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://plot_tc_shi.gs>
> <http://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>
> <http://hurricane.nf>
> <http://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>
> <http://cbar_line_box.gs>
> <http://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>
> <http://hurricane.nf> <http://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>
> <http://hurricane.nf> <http://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>
> <http://marks.nf>
> <http://marks.nf> <http://marks.nf>
> lcolor = colors.nf <http://colors.nf>
> <http://colors.nf>
> <http://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> <http://hurricane.nf>
> <http://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>
> <http://hurricane.nf> <http://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>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>>>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>>>>
>
>
> http://gradsusr.org/mailman/listinfo/gradsusr
>
>
>
>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org>>
> <mailto:gradsusr at gradsusr.org
> <mailto:gradsusr at gradsusr.org> <mailto: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>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov>
> <mailto: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>
> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>>
> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>
> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>>>.
>
>
> ********************************************************************
> Charles Seman
> Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>
> <mailto:Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>>
> <mailto:Charles.Seman at noaa.gov
> <mailto:Charles.Seman at noaa.gov> <mailto: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."
>
> *
> *
> * Subject: Re: Tracking Hurricanes
> * From: Kun-Hsuan Chou
> <cwujou at TYPHOON.AS.NTU.EDU.TW <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW>
> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW
> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW>>
> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW
> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW>
> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW
> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW>>>>
>
> * Date: Fri, 8 Oct 2004 08:03:20 +0800
> * To: GRADSUSR at LIST.CINECA.IT
> <mailto:GRADSUSR at LIST.CINECA.IT>
> <mailto:GRADSUSR at LIST.CINECA.IT
> <mailto:GRADSUSR at LIST.CINECA.IT>>
> <mailto:GRADSUSR at LIST.CINECA.IT <mailto:GRADSUSR at LIST.CINECA.IT>
> <mailto:GRADSUSR at LIST.CINECA.IT
> <mailto: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
> <mailto:cwujou at typhoon.as.ntu.edu.tw>
> <mailto:cwujou at typhoon.as.ntu.edu.tw
> <mailto:cwujou at typhoon.as.ntu.edu.tw>>
> <mailto:cwujou at typhoon.as.ntu.edu.tw
> <mailto:cwujou at typhoon.as.ntu.edu.tw>
> <mailto:cwujou at typhoon.as.ntu.edu.tw
> <mailto:cwujou at typhoon.as.ntu.edu.tw>>>
>
> *----- Original Message -----
> *From: Arturo Caracas Uribe
> *To: GRADSUSR at LIST.CINECA.IT
> <mailto:GRADSUSR at LIST.CINECA.IT>
> <mailto:GRADSUSR at LIST.CINECA.IT
> <mailto:GRADSUSR at LIST.CINECA.IT>>
> <mailto:GRADSUSR at LIST.CINECA.IT <mailto:GRADSUSR at LIST.CINECA.IT>
> <mailto:GRADSUSR at LIST.CINECA.IT
> <mailto: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
> <http://plot_hurricane_center2c.gs>
> <http://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>
> <http://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>
> <http://plot_tc_shi.gs>
> <http://plot_tc_shi.gs>"
> (Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>
> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>>
> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>
> <mailto: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>
> <http://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://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<
>
> ...
>
> [Message clipped] function maxwind(args)
> 'reinit'
> nc_file = subwrd(args,1)
> wind = subwrd(args,2)
> wind_units = subwrd(args,3)
> say
> say nc_file
> say wind' ('wind_units')'
> say
> 'sdfopen 'nc_file
> 'q file'
> say result
> say
> rec5 = sublin(result,5)
> nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz =
> subwrd(rec5,9) ;
> nt = subwrd(rec5,12)
> say
> say '(nx,ny,nz,nt) = ('nx','ny','nz','nt')'
> say
> nxp = nx
> nyp = ny
> nzp = nz
> *
> * get the x- and y-limits for virtual page
> *
> 'q gxinfo'
> rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu =
> subwrd(rec3,6)
> rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu =
> subwrd(rec4,6)
> tt = 1
> while ( tt <= nt )
> 'set x 1'
> 'set y 1'
> 'set z 1'
> 'set t 'tt
> say 'tt = 'tt
> 'set gxout print'
> '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)
> '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)
> '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 ' (xc,yc,zc) = ('xc','yc','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
> *
> * find "world" coordinates of (xc,yc) and convert "world"
> coordinates
> * to "xy" coordinates for plotting track of hurricane center...
> *
> 'set x 'xc
> lonc = subwrd(result,4)
> 'set y 'yc
> latc = subwrd(result,4)
> say ' (lonc,latc,zc) = ('lonc','latc','zc')'
> 'd 'wind'(lon='lonc',lat='latc',z='zc')'
> rec = sublin(result,2)
> wind_lnc_ltc_zc = subwrd(rec,1)
> say ' 'wind'(lon='lonc',lat='latc',z='zc') =
> 'wind_lnc_ltc_zc
> *
> * plot the wind and place a mark at the max
> *
> 'c'
> 'set vpage 'xpl' 'xpu' 'ypl' 'ypu
> 'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
> 'set x 1 'nx
> 'set y 1 'ny
> 'set z 'zc
> 'set gxout shaded'
> 'set grads off'
> 'd 'wind
> 'run cbar 0.7'
> 'draw title \ Max 'wind' = 'wind_lnc_ltc_zc' ('wind_units')\at
> (lonc,latc,zc)=('lonc','latc','zc')'
> 'q w2xy 'lonc' 'latc
> xpos = subwrd(result,3)
> ypos = subwrd(result,6)
> say ' (xpos,ypos) = ('xpos','ypos')'
> 'set line 1 1 6'
> 'draw mark 6 'xpos' 'ypos' 0.2'
> 'c events'
> pause()
> tt=tt+1
> endwhile
> say '...finished.'
> 'close 1'
> return
>
> function pause()
> *
> * Pauses till ENTER key is pressed
> *
> say
> say '*********************************'
> say 'Press "ENTER" key to continue...'
> say '*********************************'
> pull answer
> *
> return
>
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
> <mailto:gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>>
> http://gradsusr.org/mailman/listinfo/gradsusr
>
>
>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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."
>
> function maxwind(args)
> 'reinit'
> wind_nc_file = subwrd(args,1)
> wind_var = subwrd(args,2)
> wind_units = subwrd(args,3)
> lev_units = subwrd(args,4)
> say
> say wind_nc_file
> say wind_var' ('wind_units')'
> say
>
> * open "wind_nc_file" and define variables for (x,y,z,t) dimension
> limits
>
> 'sdfopen 'wind_nc_file
> 'q file'
> rec5 = sublin(result,5)
> nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ;
> nt = subwrd(rec5,12)
> say result
> say '(nx,ny,nz,nt) = ('nx','ny','nz','nt')'
> say
>
> * define (x,y,z,t) lower and upper limits for finding max wind
>
> x1 = 2 ; x2 = nx-1
> y1 = 2 ; y2 = ny-1
> z1 = 2 ; z2 = nz-1
> t1 = 2 ; t2 = nt-1
> say '(x1,x2) = ('x1','x2')'
> say '(y1,y2) = ('y1','y2')'
> say '(z1,z2) = ('z1','z2')'
> say '(t1,t2) = ('t1','t2')'
> say
>
> * get the x- and y-limits for virtual page
>
> 'q gxinfo'
> rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu = subwrd(rec3,6)
> rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu = subwrd(rec4,6)
> xptime = (xpu-xpl)/2 ; yptime = ypu-0.4
>
> * print plots into an output file?
>
> do_print = true
>
> * plot output
>
> if ( do_print = true )
> gx_file = 'maxwind_'wind_var'.gx'
> ps_file = 'maxwind_'wind_var'.ps'
> 'enable print 'gx_file
> endif
>
> * loop over time
>
> tt = t1
> while ( tt <= t2 )
>
> * find max wind and its location
>
> 'set x 1'
> 'set y 1'
> 'set z 1'
> 'set t 'tt
> 'q time' ; ttime = subwrd(result,3)
> say 'tt = 'tt' ('ttime')'
> 'set gxout print'
> 'd
> max(max(max('wind_var',x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> max_wind = subwrd(rec,1)
> say ' max 'wind_var' = 'max_wind
> 'd
> maxloc(max(max('wind_var',y='y1',y='y2'),z='z1',z='z2'),x='x1',x='x2')'
> rec = sublin(result,(z2-z1+2)*(x2-x1+1)+3)
> xc = subwrd(rec,1)
> 'd
> maxloc(max(max('wind_var',x='x1',x='x2'),z='z1',z='z2'),y='y1',y='y2')'
> rec = sublin(result,(z2-z1+2)*(y2-y1+1)+3)
> yc = subwrd(rec,1)
> 'd
> maxloc(max(max('wind_var',x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> zc = subwrd(rec,1)
> say ' (xc,yc,zc) = ('xc','yc','zc')'
> 'd 'wind_var'(x='xc',y='yc',z='zc')'
> rec = sublin(result,2)
> wind_xc_yc_zc = subwrd(rec,1)
> say ' 'wind_var'(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc
>
> * find "world" coordinates of (xc,yc) and convert "world" coordinates
> * to "xy" coordinates for plotting track of hurricane center...
>
> 'set x 'xc
> lonc = subwrd(result,4)
> 'set y 'yc
> latc = subwrd(result,4)
> 'set z 'zc
> levc = subwrd(result,4)
> say ' (lonc,latc,levc) = ('lonc','latc','levc')'
> 'd 'wind_var'(lon='lonc',lat='latc',lev='levc')'
> rec = sublin(result,2)
> wind_lnc_ltc_lvc = subwrd(rec,1)
> say ' 'wind_var'(lon='lonc',lat='latc',lev='levc') =
> 'wind_lnc_ltc_lvc
>
> * plot the wind at the level of the max and place a mark at the max
>
> 'c'
> 'set vpage 'xpl' 'xpu' 'ypl' 'ypu
> 'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
> 'set x 'x1' 'x2 ; lon1 = subwrd(result,4) ; lon2 = subwrd(result,5)
> 'set y 'y1' 'y2 ; lat1 = subwrd(result,4) ; lat2 = subwrd(result,5)
> 'set z 'zc
> 'set gxout shaded'
> 'set grads off'
> 'd 'wind_var
> 'run cbar 0.7'
> 'set strsiz 0.16'
> 'set string 1 c 6'
> 'draw string 'xptime' 'yptime' 'ttime
> 'draw title 'wind_var' ('wind_units') \ max='wind_lnc_ltc_lvc'
> at: \ (x,y,z)=('xc','yc','zc'),
> (lon,lat,lev)=('lonc'E,'latc'N,'levc''lev_units')'
> 'draw xlab plot: 'lon1'E to 'lon2'E \'
> 'draw ylab \ plot: 'lat1'N to 'lat2'N'
> 'q w2xy 'lonc' 'latc
> xpos = subwrd(result,3)
> ypos = subwrd(result,6)
> say ' (xpos,ypos) = ('xpos','ypos')'
> 'set line 1 1 6'
> 'draw mark 6 'xpos' 'ypos' 0.2'
> if ( do_print = true )
> 'print'
> endif
> 'c events'
> *** pause()
> tt=tt+1
> endwhile
> say '...finished.'
> 'close 1'
> if ( do_print = true )
> say
> 'disable print'
> '!gxps -c -i 'gx_file' -o 'ps_file
> '!ps2pdf 'ps_file
> '!rm 'gx_file' 'ps_file
> endif
> return
>
> function pause()
>
> * Pauses till ENTER key is pressed
>
> say
> say '*********************************'
> say 'Press "ENTER" key to continue...'
> say '*********************************'
> pull answer
> return
>
>
> function maxwind2(args)
> 'reinit'
> wind1_nc_file = subwrd(args,1) ; wind1_var = subwrd(args,2)
> wind2_nc_file = subwrd(args,3) ; wind2_var = subwrd(args,4)
> wind_units = subwrd(args,5)
> lev_units = subwrd(args,6)
> say
> say wind1_nc_file
> say wind1_var
> say wind2_nc_file
> say wind2_var
> say 'wind units: ('wind_units')'
> say ' lev units: ('lev_units')'
> say
>
> * open "wind1_nc_file" and "wind2_nc_file" and
> * define variables for (x,y,z,t) dimension limits
>
> 'sdfopen 'wind1_nc_file
> 'q file 1'
> rec5 = sublin(result,5)
> nx1 = subwrd(rec5,3) ; ny1 = subwrd(rec5,6) ; nz1 = subwrd(rec5,9)
> ; nt1 = subwrd(rec5,12)
> say result
> say '(nx1,ny1,nz1,nt1) = ('nx1','ny1','nz1','nt1')'
> say
> 'sdfopen 'wind2_nc_file
> 'q file 2'
> rec5 = sublin(result,5)
> nx2 = subwrd(rec5,3) ; ny2 = subwrd(rec5,6) ; nz2 = subwrd(rec5,9)
> ; nt2 = subwrd(rec5,12)
> say result
> say '(nx2,ny2,nz2,nt2) = ('nx2','ny2','nz2','nt2')'
> say
>
> * make sure (x,y,z,t) dimensions are the same for "wind1_nc_file"
> and "wind2_nc_file"
> * and, if so, define general variable names for (x,y,z,t) dimension
> limits
>
> if( nx2 = nx1 )
> nx = nx1
> else
> say
> say 'nx1 = 'nx1
> say 'nx2 = 'nx2
> say '...are not equal'
> say
> exit
> endif
> if( ny2 = ny1 )
> ny = ny1
> else
> say
> say 'ny1 = 'ny1
> say 'ny2 = 'ny2
> say '...are not equal'
> say
> exit
> endif
> if( nz2 = nz1 )
> nz = nz1
> else
> say
> say 'nz1 = 'nz1
> say 'nz2 = 'nz2
> say '...are not equal'
> say
> exit
> endif
> if( nt2 = nt1 )
> nt = nt1
> else
> say
> say 'nt1 = 'nt1
> say 'nt2 = 'nt2
> say '...are not equal'
> say
> exit
> endif
>
> * define (x,y,z,t) lower and upper limits for finding max wind
>
> x1 = 2 ; x2 = nx-1
> y1 = 2 ; y2 = ny-1
> z1 = 2 ; z2 = nz-1
> t1 = 2 ; t2 = nt-1
> say '(x1,x2) = ('x1','x2')'
> say '(y1,y2) = ('y1','y2')'
> say '(z1,z2) = ('z1','z2')'
> say '(t1,t2) = ('t1','t2')'
> say
>
> * get the x- and y-limits for virtual page
>
> 'q gxinfo'
> rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu = subwrd(rec3,6)
> rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu = subwrd(rec4,6)
> xptime = (xpu-xpl)/2 ; yptime = ypu-0.4
>
> * print plots into an output file?
>
> do_print = true
>
> * plot output
>
> if ( do_print = true )
> gx_file = 'maxwind2.gx'
> ps_file = 'maxwind2.ps <http://maxwind2.ps>'
> 'enable print 'gx_file
> endif
>
> * loop over time
>
> tt = t1
> while ( tt <= t2 )
>
> * use "wind1_nc_file" to define dimensions
>
> 'set dfile 1'
> 'set t 'tt
> 'q time' ; ttime = subwrd(result,3)
> say 'tt = 'tt' ('ttime')'
>
> * define wind magnitude variable for this time level
>
> 'set x 1 'nx
> 'set y 1 'ny
> 'set z 1 'nz
> 'define wind = mag('wind1_var'.1,'wind2_var'.2)'
>
> * find max wind and its location
>
> 'set x 1'
> 'set y 1'
> 'set z 1'
> 'set gxout print'
> 'd max(max(max(wind,x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> max_wind = subwrd(rec,1)
> say ' max wind = 'max_wind
> 'd maxloc(max(max(wind,y='y1',y='y2'),z='z1',z='z2'),x='x1',x='x2')'
> rec = sublin(result,(z2-z1+2)*(x2-x1+1)+3)
> xc = subwrd(rec,1)
> 'd maxloc(max(max(wind,x='x1',x='x2'),z='z1',z='z2'),y='y1',y='y2')'
> rec = sublin(result,(z2-z1+2)*(y2-y1+1)+3)
> yc = subwrd(rec,1)
> 'd maxloc(max(max(wind,x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> zc = subwrd(rec,1)
> say ' (xc,yc,zc) = ('xc','yc','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
>
> * find "world" coordinates of (xc,yc) and convert "world" coordinates
> * to "xy" coordinates for plotting track of hurricane center...
>
> 'set x 'xc
> lonc = subwrd(result,4)
> 'set y 'yc
> latc = subwrd(result,4)
> 'set z 'zc
> levc = subwrd(result,4)
> say ' (lonc,latc,levc) = ('lonc','latc','levc')'
> 'd wind(lon='lonc',lat='latc',lev='levc')'
> rec = sublin(result,2)
> wind_lnc_ltc_lvc = subwrd(rec,1)
> say ' wind(lon='lonc',lat='latc',lev='levc') = 'wind_lnc_ltc_lvc
>
> * plot the wind at the level of the max and place a mark at the max
>
> 'c'
> 'set vpage 'xpl' 'xpu' 'ypl' 'ypu
> 'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
> 'set x 'x1' 'x2 ; lon1 = subwrd(result,4) ; lon2 = subwrd(result,5)
> 'set y 'y1' 'y2 ; lat1 = subwrd(result,4) ; lat2 = subwrd(result,5)
> 'set z 'zc
> 'set gxout shaded'
> 'set grads off'
> 'd 'wind
> 'run cbar 0.7'
> 'set strsiz 0.16'
> 'set string 1 c 6'
> 'draw string 'xptime' 'yptime' 'ttime
> 'draw title Wind Magnitude ('wind_units') \
> max='wind_lnc_ltc_lvc' at: \ (x,y,z)=('xc','yc','zc'),
> (lon,lat,lev)=('lonc'E,'latc'N,'levc''lev_units')'
> 'draw xlab plot: 'lon1'E to 'lon2'E\'
> 'draw ylab \ plot: 'lat1'N to 'lat2'N'
> 'q w2xy 'lonc' 'latc
> xpos = subwrd(result,3)
> ypos = subwrd(result,6)
> say ' (xpos,ypos) = ('xpos','ypos')'
> 'set line 1 1 6'
> 'draw mark 6 'xpos' 'ypos' 0.2'
> if ( do_print = true )
> 'print'
> endif
> 'c events'
> 'undefine wind'
> *** pause()
> tt=tt+1
> endwhile
> say '...finished.'
> 'close 2' ; 'close 1'
> if ( do_print = true )
> say
> 'disable print'
> '!gxps -c -i 'gx_file' -o 'ps_file
> '!ps2pdf 'ps_file
> '!rm 'gx_file' 'ps_file
> endif
> return
>
> function pause()
>
> * Pauses till ENTER key is pressed
>
> say
> say '*********************************'
> say 'Press "ENTER" key to continue...'
> say '*********************************'
> pull answer
> 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."
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: maxwind2.gs
Url: http://gradsusr.org/pipermail/gradsusr/attachments/20110504/52cacc2e/attachment-0003.pl
-------------- next part --------------
A non-text attachment was scrubbed...
Name: maxwind2.pdf
Type: application/pdf
Size: 414727 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20110504/52cacc2e/attachment-0003.pdf
More information about the gradsusr
mailing list