how to make maximum wind plot along with level showing maximum wind for each grid
sushant puranik
sushantpuranik at GMAIL.COM
Wed Nov 25 03:25:43 EST 2009
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>wrote:
> Hi Phillip,
>
> Please find attached a script "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 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 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/<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 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/<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 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/<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 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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20091125/5a4ad62e/attachment.html
More information about the gradsusr
mailing list