how to make maximum wind plot along with level showing maximum wind for each grid

Charles Seman Charles.Seman at NOAA.GOV
Wed Nov 25 13:19:48 EST 2009


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."



More information about the gradsusr mailing list