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

Charles Seman Charles.Seman at NOAA.GOV
Mon Nov 30 12:22:08 EST 2009


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



More information about the gradsusr mailing list