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