how to make maximum wind plot along with level showing maximum wind for each grid
Phillip Lin
naishi at METEO.MCGILL.CA
Mon Nov 30 14:11:16 EST 2009
Hi Chuck,
Recently, I pay attention to related answering posting you did. They are
quite helpful. Now, I have almost the same problem: I have a 3-D (namely,
(y,z,t) )gridded data file which contains u, v, pressure, etc.. And I
would like to find the maximum value of the wind speed (say like
mag(u,v))and its location (i.e., (y.z)) at a fixed t. So, may I also use
the function fndlvl ?
additionally please help give a explanation about fndlvl ?
Thanks!
Phillip
> Sushant,
>
> Sorry about the way the text showed up (I copied the Example 1 command
> from http://grads.iges.org/grads/gadoc/gradfuncfndlvl.html):
>
> Here's another attempt (if you set for a horizontal (x,y) domain, you
> should get a field of pressure levels):
>
> set x ... ...
> set y ... ...
> set lev 1000
> set t ...
> set gxout shaded
> d fndlvl(mag(u,v),max(mag(u,v),lev=700,lev=1000),lev=700,lev=1000)
>
> Hope this helps,
> Chuck
>
> Charles Seman wrote:
>> Sushant,
>>
>> I'm not sure, but perhaps the "fndlvl" function may help
>> http://grads.iges.org/grads/gadoc/gradfuncfndlvl.html
>>
>> I'd try this (not sure if it would work; based on Example 1 at above web
>> site):
>>
>> |d fndlvl (|mag(u,v)|, |max(mag(u,v)|,lev=700,lev=1000||),
>> lev=700,lev=1000||)|
>>
>> ||
>>
>> |
>> |Hope this helps,
>> Chuck
>>
>> 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 <mailto: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>
>>> <mailto:Charles.Seman at noaa.gov
>>> <mailto:Charles.Seman at noaa.gov>>> wrote:
>>>
>>> Hi Phillip,
>>>
>>> Please find attached a script "fig3ab.gs <http://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>
>>> <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>
>>> <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>>.
>>>
>>>
>>>
>>>
>>> ********************************************************************
>>> 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>>
>>>
>>> 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/>
>>> <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>
>>> <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>>.
>>>
>>>
>>>
>>> ********************************************************************
>>> 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>>
>>>
>>> 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/>
>>> <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>
>>> <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>>.
>>>
>>>
>>>
>>> ********************************************************************
>>> 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>>
>>>
>>> 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/>
>>> <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> <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
>>> <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."
>>>
>>>
>>>
>>>
>>> --
>>> 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."
>
> --
>
> 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."
>
More information about the gradsusr
mailing list