Sorry to disturb you again. <br>I am not able to solve my problem till yet.<br>Since i want to study the structure of a jet stream.I am trying to find level of maximum wind.<br>At present i am able to print maximum wind speed with the help of <br>
d max(mag(u,v),lev=700,lev=1000). (Since my data is arranged in reverse order)<br><br>But i need alongwith wind speed it should print level to which maximum value corresponds.<br>Because d maxloc(mag(u,v),lev=700,lev=1000) gives location of grid to which maximum values belongs. But instead of that i need pressure level.<br>
<br>how to solve this problem.<br><br>My single file contain data of one month (in netcdf format) and data is at 6 hour interval.<br><br>If you have any solution please send me, i am trying hard for solving this particular problem. <br>
<br>Thanking you<br><br><br><br><div class="gmail_quote">On Wed, Nov 25, 2009 at 11:49 PM, Charles Seman <span dir="ltr"><<a href="mailto:Charles.Seman@noaa.gov">Charles.Seman@noaa.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Sushant,<br>
<br>
Sorry for not explaining further... The script is designed to plot time<br>
and zonally-averaged mass fluxes and estimated tropopause heights for<br>
two of GFDL's models, and uses GrADS ctl files to access datasets that<br>
you do not have... The datasets that are accessed are quite large (think<br>
too large to send via email). Maybe the script could be useful to you<br>
as an example of how to do some calculations in GrADS? What is the<br>
problem you are trying to solve? Or what diagnostic are you trying to<br>
calculate? If it is to find the level of maximum wind, the code (see<br>
earlier email) to find the level of minimum temperature could be<br>
modified to find the level of maximum wind (the code was offered as an<br>
example of how to use GrADS to do this class of problems and not to<br>
solve your specific problem...).<div class="im"><br>
<br>
Hope this helps,<br>
Chuck<br>
<br>
sushant puranik wrote:<br>
</div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div class="im">
Thank you for immediate reply i am testing those scripts. but while<br>
running the script it shows error in the information of GrADS files<br>
<br>
level_dir = '/net2/cjs/ljd/lima/ch3i1/pp/atmos_level/av/monthly_16yr'<br>
don_dir = '/net2/cjs/ljd/lima/ch3i1/pp/atmosdon/av/monthly_16yr'<br>
ras_dir =<br>
'/net2/cjs/lwh/fms/lima/m45_am2p14_ch3i/pp/atmos_level/av/monthly_16yr'<br>
<br>
sdir_months = 'Aug Sep'<br>
month_days = '31 30'<br>
<br>
/home/Sushant/Desktop/level_ctl = 'mc.o_zl.ctl' ; level_var = 'mc'<br>
don_ctl1 = 'uceml_deep.o_zl.ctl' ; don_var1 = 'uceml_deep'<br>
don_ctl2 = 'umeml_deep.o_zl.ctl' ; don_var2 = 'umeml_deep'<br>
don_ctl3 = 'dmeml_deep.o_zl.ctl' ; don_var3 = 'dmeml_deep'<br>
<br>
ras_ctl = 'mc.o_zl.ctl' ; ras_var = 'mc'<br>
<br>
temp_ctl = 'temp.o_zl.ctl' ; temp_var = 'temp'<br>
<br>
So can you please send me those files.<br>
<br>
Also i am using ECMWF data in the netcdf format. So there is no need<br>
to use ctl files. It is possible to run the script without using such<br>
ctl file info.<br>
<br>
Sushant<br>
<br>
<br>
On Wed, Nov 25, 2009 at 1:01 AM, Charles Seman <<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div><div class="im">
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>>> wrote:<br>
<br>
Hi Phillip,<br>
<br></div>
Please find attached a script "<a href="http://fig3ab.gs" target="_blank">fig3ab.gs</a> <<a href="http://fig3ab.gs" target="_blank">http://fig3ab.gs</a>>" and<div class="im"><br>
called script functions<br>
"pause.gsf" and "scale_zdef_zlevs.gsf". "scale_zdef_zlevs.gsf"<br>
makes a<br>
local (temporary) ctl file with scaled vertical levels for plotting<br>
purposes: here specifically to convert the original zdef zlevs values<br>
from meters to kilometers (there may be a better way to do this; maybe<br>
using commands "set ylab off", followed by "set yaxis ..."). I<br>
hope all<br></div>
of the scripts needed for the <a href="http://fig3ab.gs" target="_blank">fig3ab.gs</a> <<a href="http://fig3ab.gs" target="_blank">http://fig3ab.gs</a>> script<div><div></div><div class="h5"><br>
are attached; if not,<br>
please let me know. You may have questions about these scripts;<br>
please<br>
contact me if you do...<br>
<br>
<br>
Hope this helps,<br>
Chuck<br>
<br>
Phillip Lin wrote:<br>
<br>
Hi Chuck,<br>
Thanks for your dedicated input to this problem we have.<br>
And could you be kind to send me a script you did before for<br>
the sane or<br>
similar question? Then, maybe it is much helpful for me to<br>
solve this<br>
problem!<br>
Thanks much!<br>
Phillip<br>
<br>
<br>
Hi Phillip,<br>
<br>
Sorry I didn't explain further: the code was meant as an<br>
example of how<br>
to find and plot extrema in a field of data. To do the<br>
max wind, the<br>
code would have to be modified to search in a wind speed<br>
field using<br>
"maxloc"... and the exact code would also depend on the<br>
domain...<br>
<br>
Hope this helps,<br>
Chuck<br>
<br>
Phillip Lin wrote:<br>
<br>
Hi Chuck,<br>
I also met this problem. And there are something<br>
unclear to me (only to<br>
me!)in your last email, because:<br>
(1) the level that has minimum temp isn't the same as<br>
the level where<br>
there is maximum zonal wind speed;<br>
(2) what does the meaning of the step (3) you said;<br>
sorry how can plot<br>
the<br>
values (xloc,j,yloc,j) on a (y,z)plot (this is the<br>
cross section in the<br>
height-latitude plane?<br>
(3)could you be kind to tell me the basic ideas to get<br>
the maximum jet<br>
stream, or... this means the contours of zonal wind<br>
speed in a specific<br>
level as 250mb or ...?<br>
(4)could you be kind to send me the script to do this<br>
plotting, and data<br>
file too for a test.<br>
I am quite confused about this! Thank you so much in<br>
advance!<br>
Phillip<br>
<br>
<br>
<br>
<br>
Sushant,<br>
<br>
Here's some code that:<br>
1) finds the height level "zmin" of the minimum<br>
temperature<br>
(represented<br>
by variable ztav = zonal average of weighted time<br>
average temperature)<br>
at each horizontal location "j" in a (y,z)<br>
physical domain<br>
2) uses GrADS function "gr2xy"<br>
(<a href="http://grads.iges.org/grads/gadoc/script.html#commands" target="_blank">http://grads.iges.org/grads/gadoc/script.html#commands</a>)<br>
to transform<br>
the grid point location (j,zmin) to plot location<br>
(xloc.j,yloc.j)<br>
3) plots the values (xloc.j,yloc.j) on a (y,z) plot<br>
<br>
say ' find and save location of minimum<br>
temperature...'<br>
<br>
j=1<br>
while ( j <= ny )<br>
'set y 'j<br>
'set gxout print'<br>
'd minloc(ztav,z=1,z='nz')' ; rec3 =<br>
sublin(result,3) ; zmin =<br>
subwrd(rec3,1)<br>
* say result<br>
say ' ...j, zmin = 'j', 'zmin<br>
'query gr2xy 'j' 'zmin<br>
rec = sublin(result,1) ; xloc.j =<br>
subwrd(rec,3) ; yloc.j =<br>
subwrd(rec,6)<br>
* say result<br>
say ' xloc, yloc = 'xloc.j', 'yloc.j<br>
j=j+1<br>
endwhile<br>
...<br>
if( plot_min_temp = 'yes' )<br>
say ' plot level of minimum temperature... '<br>
j=1<br>
while ( j <= ny )<br>
'set y 'j<br>
* say ' xloc, yloc = 'xloc.j', 'yloc.j<br>
'set line 1'<br>
'draw mark 3 'xloc.j' 'yloc.j' 0.03'<br>
j=j+1<br>
endwhile<br>
endif<br>
<br>
Please contact me if you have any questions about<br>
this code, or would<br>
like a script that does a plot.<br>
<br>
Hope this helps,<br>
Chuck<br>
<br>
<br>
sushant puranik wrote:<br>
<br>
<br>
Hello Sir<br>
I want to study the structure of a jet stream<br>
for that purpose i wish<br>
to make a plot showing maximum wind speed for<br>
each grid box alongwith<br>
the level in which the speed is maximum. Means<br>
it should print level<br>
with maximum wind speed value. I am trying to<br>
solve this problem since<br>
last 15 days but i am not able to solve it.<br>
Does anyone has any idea<br>
how i can make such plot using GrADS. I am<br>
using ECMWF interim<br>
reanalysis data for this purpose.<br>
<br>
Any suggestion is appreciated<br>
<br>
Thanking you.<br>
<br>
--<br>
Sushant Puranik<br>
Junior Research Fellow<br>
Dept. of Atmospheric & Space Sciences,<br>
University of Pune,<br>
Pune-07,<br>
India.<br>
<br>
<br>
--<br>
<br>
Please note that <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div></div>
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>> should be<div class="im"><br>
considered my NOAA<br>
email address, not <a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a><br></div>
<mailto:<a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a>>.<div class="im"><br>
<br>
********************************************************************<br>
Charles Seman<br>
<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div>
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>><div class="im"><br>
U.S. Department of Commerce / NOAA / OAR<br>
Geophysical Fluid Dynamics Laboratory<br>
voice: (609) 452-6547<br>
201 Forrestal Road<br>
fax: (609) 987-5063<br>
Princeton, NJ 08540-6649<br>
<a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/~cjs/</a><br></div>
<<a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/%7Ecjs/</a>><div class="im"><br>
********************************************************************<br>
<br>
"The contents of this message are mine personally<br>
and do not<br>
necessarily<br>
reflect any position of the Government or NOAA."<br>
<br>
<br>
<br>
--<br>
<br>
Please note that <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div>
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>> should be considered my NOAA<div class="im"><br>
email address, not <a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a><br></div>
<mailto:<a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a>>.<br>
<br>
********************************************************************<br>
Charles Seman<br>
<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a> <mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>><div class="im"><br>
U.S. Department of Commerce / NOAA / OAR<br>
Geophysical Fluid Dynamics Laboratory voice:<br>
(609) 452-6547<br>
201 Forrestal Road fax:<br>
(609) 987-5063<br>
Princeton, NJ 08540-6649<br>
<a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/~cjs/</a><br></div>
<<a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/%7Ecjs/</a>><div class="im"><br>
********************************************************************<br>
<br>
"The contents of this message are mine personally and do<br>
not necessarily<br>
reflect any position of the Government or NOAA."<br>
<br>
<br>
<br>
--<br>
<br>
Please note that <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div>
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>> should be considered my NOAA<br>
email address, not <a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a> <mailto:<a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a>>.<br>
<br>
********************************************************************<br>
Charles Seman<br>
<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a> <mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>><div class="im"><br>
U.S. Department of Commerce / NOAA / OAR<br>
Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547<br>
201 Forrestal Road fax: (609) 987-5063<br>
Princeton, NJ 08540-6649<br></div>
<a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/~cjs/</a> <<a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/%7Ecjs/</a>><div><div></div><div class="h5">
<br>
********************************************************************<br>
<br>
"The contents of this message are mine personally and do not<br>
necessarily<br>
reflect any position of the Government or NOAA."<br>
<br>
<br>
************************************************************************<br>
*<br>
* To use: "grads -p"<br>
*<br>
* Script calculates and plots zonally averaged cross-sections of the<br>
* time averaged height interpolated "mass flux" from the GrADS files<br>
* specified below (with an option to overlay the level of minimum<br>
* temperature):<br>
*<br>
************************************************************************<br>
function mass_flux (args)<br>
<br>
'reinit'<br>
<br>
************************************************************************<br>
* Allow external GrADS functions<br>
************************************************************************<br>
<br>
rc = gsfallow("on")<br>
<br>
************************************************************************<br>
* Pause between plots?<br>
************************************************************************<br>
<br>
do_pause = 'no'<br>
<br>
************************************************************************<br>
* GrADS meta file output?<br>
************************************************************************<br>
<br>
do_print = 'yes'<br>
<br>
************************************************************************<br>
* Information for GrADS files<br>
************************************************************************<br>
<br>
level_dir = '/net2/cjs/ljd/lima/ch3i1/pp/atmos_level/av/monthly_16yr'<br>
don_dir = '/net2/cjs/ljd/lima/ch3i1/pp/atmosdon/av/monthly_16yr'<br>
ras_dir =<br>
'/net2/cjs/lwh/fms/lima/m45_am2p14_ch3i/pp/atmos_level/av/monthly_16yr'<br>
<br>
sdir_months = 'Aug Sep'<br>
month_days = '31 30'<br>
<br>
level_ctl = 'mc.o_zl.ctl' ; level_var = 'mc'<br>
don_ctl1 = 'uceml_deep.o_zl.ctl' ; don_var1 = 'uceml_deep'<br>
don_ctl2 = 'umeml_deep.o_zl.ctl' ; don_var2 = 'umeml_deep'<br>
don_ctl3 = 'dmeml_deep.o_zl.ctl' ; don_var3 = 'dmeml_deep'<br>
<br>
ras_ctl = 'mc.o_zl.ctl' ; ras_var = 'mc'<br>
<br>
temp_ctl = 'temp.o_zl.ctl' ; temp_var = 'temp'<br>
<br>
************************************************************************<br>
* Information for the plots<br>
************************************************************************<br>
<br>
plot_min_max = 'no'<br>
*<br>
* Mass flux variable scale and units<br>
*<br>
var_scale = 1000<br>
var_units = '(g m`a-2`n s`a-1`n)'<br>
*<br>
* Height levels (scale GrADS ctl file levels from "m" to "km")<br>
*<br>
zscale = 1./1000.<br>
zlevs = '0 18'<br>
ylabel_int = 2<br>
ylabel = 'Height (km)'<br>
*<br>
* Define the contour levels...<br>
*<br>
*contour_levels = '1 3 5 7 9 11 13'<br>
*contour_levels = '0 2 4 6 8 10 12 14 16'<br>
contour_levels = '0 0.01 0.1 1 3 5 7 10 15'<br>
*<br>
* Plot level of minimum temperature?<br>
*<br>
plot_min_temp = 'yes'<br>
<br>
************************************************************************<br>
* Plot variables' page limits<br>
************************************************************************<br>
<br>
pvar_vpage.1 = '0 8.5 4.8 9.3'<br>
pvar_vpage.2 = '0 8.5 0.3 4.8'<br>
pvar_parea = '1 7.0 0.5 4.0'<br>
<br>
xl_pvar_parea = subwrd(pvar_parea,1)<br>
yt_pvar_parea = subwrd(pvar_parea,4)<br>
<br>
************************************************************************<br>
* Title information for the plots<br>
************************************************************************<br>
<br>
mod_title = 'FMS lima 1983-1998'<br>
var_title = 'Mass Flux 'var_units<br>
avg_title = 'Zonal average of weighted time average'<br>
<br>
figure_panel.1 = '(a)' ; model.1 = 'AM2-D'<br>
figure_panel.2 = '(b)' ; model.2 = 'AM2'<br>
*<br>
* Define names for the "sdir_months" subdirectories...<br>
*<br>
define_months()<br>
<br>
************************************************************************<br>
* GrADS metafile for output<br>
************************************************************************<br>
<br>
if( do_print = 'yes' )<br>
'enable print fig3ab.gx'<br>
endif<br>
<br>
*-----------------------------------------------------------------------<br>
* Determine number of months to average...<br>
*-----------------------------------------------------------------------<br>
<br>
n=1<br>
month = subwrd(sdir_months,n)<br>
if( month = '' )<br>
say<br>
say 'No months specified in variable "sdir_months"...'<br>
say 'Exiting script'<br>
say<br>
exit<br>
else<br>
while ( month != '' )<br>
n=n+1<br>
month = subwrd(sdir_months,n)<br>
endwhile<br>
nmonths=n-1<br>
endif<br>
<br>
first = subwrd(sdir_months,1) ; last = subwrd(sdir_months,nmonths)<br>
months = ''_name.first'-'_name.last''<br>
<br>
*-----------------------------------------------------------------------<br>
* Verify height levels...<br>
*-----------------------------------------------------------------------<br>
<br>
l=1<br>
zlv.l = subwrd(zlevs,l)<br>
if( zlv.l = '' )<br>
say<br>
say 'No pressure levels specified in variable "zlevs"...'<br>
say 'Exiting script'<br>
say<br>
exit<br>
else<br>
while ( zlv.l != '' )<br>
l=l+1<br>
zlv.l = subwrd(zlevs,l)<br>
endwhile<br>
nl = l-1<br>
if( nl > 2 )<br>
say<br>
say 'More than two height levels ("nl" = 'nl') specified in<br>
variable "zlevs"...'<br>
say 'Exiting script'<br>
say<br>
exit<br>
else<br>
if( zlv.2 <= zlv.1 )<br>
say<br>
say 'Height level "2" (= 'zlv.2') not greater than height<br>
level "1" (= 'plv.1')...'<br>
say 'Exiting script'<br>
say<br>
exit<br>
endif<br>
endif<br>
endif<br>
<br>
say<br>
say<br>
'------------------------------------------------------------------'<br>
say 'Calculate zonal averages of weighted time average "mass<br>
fluxes"...'<br>
say<br>
'------------------------------------------------------------------'<br>
say<br>
say 'Initialize time sum variables for time average...'<br>
say<br>
*<br>
* Use GrADS ctl files for the first "month" to define time sum<br>
variables...<br>
*<br>
month1 = subwrd(sdir_months,1)<br>
say<br>
say '--------------------'<br>
say 'Model: 'model.1<br>
say '--------------------'<br>
say<br>
<br>
ctl_file = ''level_dir'/'month1'/'level_ctl<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'mass_flux.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
say ' initialize: wtsum1m'<br>
'define wtsum1m = const('level_var',0,-a)'<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
if( plot_min_temp = 'yes' )<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file = ''level_dir'/'month1'/'temp_ctl<br>
ctl_file_scale = 'temp.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nxt = subwrd(rec5,3) ; nyt = subwrd(rec5,6) ; nzt = subwrd(rec5,9)<br>
if( nxt != nx )<br>
say '...nxt = 'nxt<br>
say '...nx = 'nx<br>
exit<br>
endif<br>
if( nyt != ny )<br>
say '...nyt = 'nyt<br>
say '...ny = 'ny<br>
exit<br>
endif<br>
if( nzt != nz )<br>
say '...nzt = 'nzt<br>
say '...nz = 'nz<br>
exit<br>
endif<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
say ' initialize: wtsum1t'<br>
'define wtsum1t = const('temp_var',0,-a)'<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
endif<br>
<br>
say<br>
say '--------------------'<br>
say 'Model: 'model.2<br>
say '--------------------'<br>
say<br>
<br>
ctl_file = ''ras_dir'/'month1'/'ras_ctl<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'mass_flux.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
say ' initialize: wtsum2m'<br>
'define wtsum2m = const('ras_var',0,-a)'<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
if( plot_min_temp = 'yes' )<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file = ''ras_dir'/'month1'/'temp_ctl<br>
ctl_file_scale = 'temp.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nxt = subwrd(rec5,3) ; nyt = subwrd(rec5,6) ; nzt = subwrd(rec5,9)<br>
if( nxt != nx )<br>
say '...nxt = 'nxt<br>
say '...nx = 'nx<br>
exit<br>
endif<br>
if( nyt != ny )<br>
say '...nyt = 'nyt<br>
say '...ny = 'ny<br>
exit<br>
endif<br>
if( nzt != nz )<br>
say '...nzt = 'nzt<br>
say '...nz = 'nz<br>
exit<br>
endif<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
say ' initialize: wtsum2t'<br>
'define wtsum2t = const('temp_var',0,-a)'<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
endif<br>
<br>
if( do_pause = 'yes' )<br>
pause()<br>
endif<br>
*<br>
* Calculate running time sum variables for "mass flux" variable...<br>
*<br>
n=1<br>
ndays=0<br>
*-----------------------------------------------------------------------<br>
while ( n <= nmonths )<br>
*-----------------------------------------------------------------------<br>
*<br>
* Name of the month...<br>
*<br>
month = subwrd(sdir_months,n)<br>
*<br>
* Days in the month...<br>
*<br>
days = subwrd(month_days,n)<br>
<br>
say<br>
say '------------------------'<br>
say 'Month = 'month', days = 'days<br>
say '------------------------'<br>
say<br>
say<br>
say '--------------------'<br>
say 'Model: 'model.1<br>
say '--------------------'<br>
say<br>
<br>
ctl_file.1 = ''level_dir'/'month'/'level_ctl<br>
ctl_file.2 = ''don_dir'/'month'/'don_ctl1<br>
ctl_file.3 = ''don_dir'/'month'/'don_ctl2<br>
ctl_file.4 = ''don_dir'/'month'/'don_ctl3<br>
*<br>
* Make local ctl files with scaled vertical levels...<br>
*<br>
ctl_file_scale.1 = 'mass_flux.1.scale.ctl'<br>
scale_zdef_zlevs(ctl_file.1,zscale,ctl_file_scale.1)<br>
ctl_file_scale.2 = 'mass_flux.2.scale.ctl'<br>
scale_zdef_zlevs(ctl_file.2,zscale,ctl_file_scale.2)<br>
ctl_file_scale.3 = 'mass_flux.3.scale.ctl'<br>
scale_zdef_zlevs(ctl_file.3,zscale,ctl_file_scale.3)<br>
ctl_file_scale.4 = 'mass_flux.4.scale.ctl'<br>
scale_zdef_zlevs(ctl_file.4,zscale,ctl_file_scale.4)<br>
<br>
say 'Open GrADS file 1:'<br>
say '...'ctl_file_scale.1<br>
'open 'ctl_file_scale.1<br>
say 'Open GrADS file 2:'<br>
say '...'ctl_file_scale.2<br>
'open 'ctl_file_scale.2<br>
say 'Open GrADS file 3:'<br>
say '...'ctl_file_scale.3<br>
'open 'ctl_file_scale.3<br>
say 'Open GrADS file 4:'<br>
say '...'ctl_file_scale.4<br>
'open 'ctl_file_scale.4<br>
say 'GrADS files ready'<br>
<br>
say<br>
say ' calculate weighted (using days = 'days') "mass flux" time<br>
sum'<br>
say ' ('level_var'.1+'don_var1'.2+'don_var2'.3+'don_var3'.4)'<br>
say ' for current month and add it to the time sum variable...'<br>
say<br>
'set dfile 1'<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define wtsum1m = wtsum1m +<br>
'days'*('level_var'.1+'don_var1'.2+'don_var2'.3+'don_var3'.4)'<br>
<br>
'close 4' ; 'close 3' ; 'close 2' ; 'close 1'<br>
'!rm 'ctl_file_scale.1<br>
'!rm 'ctl_file_scale.2<br>
'!rm 'ctl_file_scale.3<br>
'!rm 'ctl_file_scale.4<br>
<br>
if( plot_min_temp = 'yes' )<br>
<br>
ctl_file = ''level_dir'/'month'/'temp_ctl<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'temp.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
say 'GrADS file ready'<br>
say<br>
say ' calculate weighted (using days = 'days') temperature<br>
time sum'<br>
say ' for current month and add it to the time sum variable...'<br>
say<br>
<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define wtsum1t = wtsum1t + 'days'*'temp_var<br>
<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
endif<br>
<br>
say<br>
say '--------------------'<br>
say 'Model: 'model.2<br>
say '--------------------'<br>
say<br>
<br>
ctl_file = ''ras_dir'/'month'/'ras_ctl<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'mass_flux.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
say 'GrADS file ready'<br>
<br>
say<br>
say ' calculate weighted (using days = 'days') "mass flux" time<br>
sum'<br>
say ' ('ras_var'.1)'<br>
say ' for current month and add it to the time sum variable...'<br>
say<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define wtsum2m = wtsum2m + 'days'*'ras_var<br>
<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
if( plot_min_temp = 'yes' )<br>
<br>
ctl_file = ''ras_dir'/'month'/'temp_ctl<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'temp.scale.ctl'<br>
scale_zdef_zlevs(ctl_file,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
say 'GrADS file ready'<br>
say<br>
say ' calculate weighted (using days = 'days') temperature<br>
time sum'<br>
say ' for current month and add it to the time sum variable...'<br>
say<br>
<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define wtsum2t = wtsum2t + 'days'*'temp_var<br>
<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
endif<br>
<br>
say<br>
say '...finished month'<br>
say<br>
<br>
if( do_pause = 'yes' )<br>
pause()<br>
endif<br>
<br>
n=n+1<br>
ndays=ndays+days<br>
*-----------------------------------------------------------------------<br>
endwhile<br>
*-----------------------------------------------------------------------<br>
say<br>
say '---------------------------'<br>
say 'Finished weighted time sums'<br>
say '---------------------------'<br>
say<br>
say<br>
say<br>
'------------------------------------------------------------------'<br>
say 'Calculate zonal averages of the weighted time averages'<br>
say 'of the weighted time sums, and plot the data...'<br>
say<br>
'------------------------------------------------------------------'<br>
say<br>
say 'Use GrADS ctl files for the first "month" to define time<br>
averaged variables for plotting...'<br>
*<br>
* Model #1<br>
*<br>
ctl_file.1 = ''level_dir'/'month1'/'level_ctl<br>
ctl_file_temp.1 = ''level_dir'/'month1'/'temp_ctl<br>
*<br>
* Model #2<br>
*<br>
ctl_file.2 = ''ras_dir'/'month1'/'ras_ctl<br>
ctl_file_temp.2 = ''ras_dir'/'month1'/'temp_ctl<br>
*<br>
* Initialize min/max variables...<br>
*<br>
min = 1e6 ; max = -1e6<br>
*<br>
* Start a new frame, and plot a title at top of page...<br>
*<br>
ptitle_top(var_title,avg_title,months,pvar_parea)<br>
<br>
m=1<br>
while ( m <= 2 )<br>
<br>
say<br>
say '--------------------'<br>
say 'Model: 'model.m<br>
say '--------------------'<br>
say<br>
<br>
if( plot_min_temp = 'yes' )<br>
say<br>
say 'Find and save level of minimum temperature for the'<br>
say 'zonal and weighted time average temperature...'<br>
say<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'temp.scale.ctl'<br>
scale_zdef_zlevs(ctl_file_temp.m,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
say 'GrADS file ready'<br>
say<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)<br>
<br>
say ' calculate weighted time average (using ndays = 'ndays')...'<br>
<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define tav = wtsum'm't/'ndays<br>
<br>
say ' calculate zonal averages of the weighted time average...'<br>
<br>
'set x 1'<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define ztav = ave(tav,x=1,x='nx')'<br>
<br>
say ' plot zonal average of the weighted time average...'<br>
<br>
'set vpage 'pvar_vpage.m<br>
'set parea 'pvar_parea<br>
'set grads off'<br>
'set gxout contour'<br>
'set clab off'<br>
'set ccolor 0'<br>
'set xlab off'<br>
'set ylab off'<br>
'd ztav'<br>
<br>
say ' find and save location of minimum temperature...'<br>
<br>
j=1<br>
while ( j <= ny )<br>
'set y 'j<br>
'set gxout print'<br>
'd minloc(ztav,z=1,z='nz')' ; rec3 = sublin(result,3) ; zmin<br>
= subwrd(rec3,1)<br>
* say result<br>
say ' ...j, zmin = 'j', 'zmin<br>
'query gr2xy 'j' 'zmin<br>
rec = sublin(result,1) ; xloc.j = subwrd(rec,3) ; yloc.j =<br>
subwrd(rec,6)<br>
* say result<br>
say ' xloc, yloc = 'xloc.j', 'yloc.j<br>
j=j+1<br>
endwhile<br>
<br>
'undefine wtsum'm't' ; 'undefine tav' ; 'undefine ztav'<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
endif<br>
<br>
say<br>
say 'Calculate and plot zonally averaged "mass flux"...'<br>
say<br>
*<br>
* Make local ctl file with scaled vertical levels...<br>
*<br>
ctl_file_scale = 'mass_flux.scale.ctl'<br>
scale_zdef_zlevs(ctl_file.m,zscale,ctl_file_scale)<br>
<br>
say 'Open GrADS file:'<br>
say '...'ctl_file_scale<br>
'open 'ctl_file_scale<br>
say 'GrADS file ready'<br>
say<br>
say<br>
'query file'<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)<br>
<br>
say ' calculate weighted "mass flux" time average (using ndays<br>
= 'ndays')'<br>
<br>
'set x 1 'nx<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define tav = wtsum'm'm/'ndays<br>
<br>
say ' calculate zonal averages of the weighted time average...'<br>
<br>
'set x 1'<br>
'set y 1 'ny<br>
'set z 1 'nz<br>
'set t 1'<br>
'define ztav = 'var_scale'*ave(tav,x=1,x='nx')'<br>
<br>
say ' plot zonal averages of time averaged "mass flux"'<br>
'set vpage 'pvar_vpage.m<br>
'set parea 'pvar_parea<br>
'set grads off'<br>
'set clevs 'contour_levels<br>
'set gxout shaded'<br>
'set xlab on'<br>
'set ylab on'<br>
'set ylint 'ylabel_int<br>
'set x 1'<br>
'set y 1 'ny<br>
'set lev 'zlevs<br>
'set t 1'<br>
'd ztav'<br></div></div>
'run <a href="http://cbar.gs" target="_blank">cbar.gs</a> <<a href="http://cbar.gs" target="_blank">http://cbar.gs</a>> 0.7'<div><div></div><div class="h5"><br>
'draw title \\\ 'model.m<br>
* 'draw xlab Latitude \\\'<br>
'draw ylab \\\ 'ylabel<br>
xfp = xl_pvar_parea + 0.1<br>
yfp = yt_pvar_parea + 0.2<br>
'set string 1 l 6'<br>
'set strsiz 0.13'<br>
'draw string 'xfp' 'yfp' 'figure_panel.m<br>
'set gxout stat'<br>
'd ztav'<br>
min_max = sublin(result,8)<br>
say ' 'min_max<br>
min = subwrd(min_max,4)<br>
max = subwrd(min_max,5)<br>
if( plot_min_max = 'yes' )<br>
xmm = xl_pvar_parea + 0.1<br>
ymm = yt_pvar_parea - 0.1<br>
'set string 1 l 4'<br>
'set strsiz 0.08'<br>
'draw string 'xmm' 'ymm' 'min_max' '<br>
endif<br>
<br>
if( plot_min_temp = 'yes' )<br>
say ' plot level of minimum temperature... '<br>
j=1<br>
while ( j <= ny )<br>
'set y 'j<br>
* say ' xloc, yloc = 'xloc.j', 'yloc.j<br>
'set line 1'<br>
'draw mark 3 'xloc.j' 'yloc.j' 0.03'<br>
j=j+1<br>
endwhile<br>
endif<br>
<br>
'undefine wtsum'm'm' ; 'undefine tav' ; 'undefine ztav'<br>
'close 1'<br>
'!rm 'ctl_file_scale<br>
<br>
say<br>
say '...finished plot for: 'model.m<br>
say<br>
m=m+1<br>
endwhile<br>
<br>
say<br>
say '--------------'<br>
say 'Finished plots'<br>
say '--------------'<br>
say<br>
<br>
if( do_print = 'yes' )<br>
'print'<br>
'disable print'<br>
endif<br>
<br>
say<br>
say '--------------------------'<br>
say '"Mass flux" min/max values'<br>
say '--------------------------'<br>
say<br>
say 'min = 'min' max = 'max<br>
<br>
say<br>
say '************************************************'<br>
say ' Finished with this script.'<br>
say '************************************************'<br>
say<br>
<br>
if( batch_mode = 'yes' )<br>
quit<br>
endif<br>
<br>
return<br>
<br>
************************************************************************<br>
*<br>
function define_months()<br>
*<br>
* Define "names" for all months...<br>
*<br>
_name.Jan = 'January'<br>
_name.Feb = 'February'<br>
_name.Mar = 'March'<br>
_name.Apr = 'April'<br>
_name.May = 'May'<br>
_name.Jun = 'June'<br>
_name.Jul = 'July'<br>
_name.Aug = 'August'<br>
_name.Sep = 'September'<br>
_name.Oct = 'October'<br>
_name.Nov = 'November'<br>
_name.Dec = 'December'<br>
*<br>
return<br>
*<br>
************************************************************************<br>
<br>
************************************************************************<br>
*<br>
function ptitle_top(title,subtitle,months,plim)<br>
*<br>
* Plots a title at the top of the page<br>
*<br>
x1 = subwrd(plim,1)<br>
x2 = subwrd(plim,2)<br>
xc = x1 + (x2-x1)/2<br>
'c'<br>
'set vpage 0 8.5 0 11'<br>
'set string 1 c 6'<br>
'set strsiz 0.14'<br>
'draw string 'xc' 10.2 'title<br>
*'draw string 'xc' 10.2 'subtitle<br>
'draw string 'xc' 9.8 'months<br>
*<br>
return<br>
*<br>
************************************************************************<br>
<br>
<br>
<br>
<br>
--<br>
Sushant Puranik<br>
Junior Research Fellow<br>
Dept. of Atmospheric & Space Sciences,<br>
University of Pune,<br>
Pune-07,<br>
India.<br>
</div></div></blockquote>
<br><div><div></div><div class="h5">
--<br>
<br>
Please note that <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a> should be considered my NOAA<br>
email address, not <a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a>.<br>
<br>
********************************************************************<br>
Charles Seman <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br>
U.S. Department of Commerce / NOAA / OAR<br>
Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547<br>
201 Forrestal Road fax: (609) 987-5063<br>
Princeton, NJ 08540-6649 <a href="http://www.gfdl.noaa.gov/%7Ecjs/" target="_blank">http://www.gfdl.noaa.gov/~cjs/</a><br>
********************************************************************<br>
<br>
"The contents of this message are mine personally and do not necessarily<br>
reflect any position of the Government or NOAA."<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Sushant Puranik<br>Junior Research Fellow<br>Dept. of Atmospheric & Space Sciences,<br>University of Pune,<br>Pune-07,<br>India.<br>