[gradsusr] Finding the maximum value of Ri in a single grid point

[Jeffrey Duda]

The problem is that there is a difference between grids that can be
displayed and script variables.  They can't be handled the same way.  You
can't use if statements using displayed fields.  What you need to do is to
save the output of each computation in your script to a script variable
using the subwrd command.  Here is an example:

'define th1 = temp * (pow(100000 / press , (0.286)))'
th1 = subwrd(result,4)

Hopefully you've noticed this before, but when you enter a display command
at the Grads prompt, there will be a small amount of text output.  If you
display a variable with all of the dimensions in your environment fixed,
the text will include the value of the field you are displaying.  What
subwrd does is grab the fourth word in the string variable "result".
"result" is automatically assigned each time a display command is run, so
that's why you don't see it being assigned anywhere else.  NOTE: for some
data sets and for some environment settings, Grads will interpolate the
output before displaying, so the text output from the display command will
read

"Note: interpolating variable...
Result value: _____"

In that case you'll need to use the sublin command first:
line = sublin(result,2)
th1 = subwrd(line,4)
This tells Grads to nab the second line of the output from the result
variable (result includes all of the text output from the display
command).  Then you use subwrd to nab the 4th word of the string variable
"line".  When you do this, you'll have the fields you want as Grads
scripting variables.  Then you can run if statements on them.  To output
the value of a Grads script variable, use the say command:
say th1 OR
say "th1 = "th1

ON THE OTHER HAND, you can do this without using any scripting variables.
In this case use the max/min and maxloc/minloc functions.  Then your script
would look something like this:

'set lev 'botlevel' 'toplevel

 'define th1 = temp * (pow(100000 / press , (0.286)))'
'define u1 = u'
 'define v1 = v'
'define z1 = geopot/9.81'

 'define th2 = temp(z-1) * (pow(100000 / press(z-1) , (0.286)))'
'define u2 = u(z-1)'
 'define v2 = v(z-1)'
'define z2 = geopot(z-1)/9.81'

'define
ri=9.81/(0.5*th1+th2)*(th2-th1)/(pow((u2-u1),2)+pow((v2-v1),2))*(z2-z1)'

'd ri'
say result

'd min(ri,lev='botlev',lev='toplev')'
'd minloc(ri,lev='botlev',lev='toplev')'
'd max(ri,lev='botlev',lev='toplev')'
'd maxloc(ri,lev='botlev',lev='toplev')'

Carefully follow the location of the single quotes and note that there is
no while loop here.  Also note the (z-1) appended to the *2 variables to
represent the field one level below the current level.  The whole quote
thing confuses some people, so that's understandable.  Try this and see if
it works.  Ask if you have any questions.  Good luck.

Jeff Duda

On Thu, Nov 3, 2011 at 5:23 AM, Lippo Nester <liikkuvattaakse at gmail.com>wrote:

> Hi!
>
> I've been trying to create my very first GrADS-script. The script should
> compute values of the bulk Richardson number in a specific grid point
> between predefined model levels and then find the minimum value of Ri and
> probably the model level associated with it.
>
> I tried to do this using while and if loops, but for some reason it's not
> working as I expected. Maybe I can't understand handling variables in GrADS?
>
> In this script the top and the bottom model layers are set to 40 and 50.
> Ri is computed and I get a list of the Ri at different levels, and they
> seem to be ok. However the variable "minri" gets always the value of ri at
> the lowest model level (here: ri at level 50). Variable "minlevel" is
> always the highest model level (here: level 40)!
>
> When count is 40 at the first step of the while-loop, the value of ri
> should be saved into variable minri and the model level at minlev. After
> that, if ri < minri at lower levels (count > 40), values of minri and
> minlev should be varied by the second if-loop...
>
> Is there a more simple or more correct way to do this? I thought max()
> can't make it...
>
> Any help is much appreciated!
>
> Best regards, Lippo
>
> (version of GrADS 2.0.a1)
>
> 'set lat 50.00'
> 'set lon 10.00'
> botlevel=50
> toplevel=40
>
> count=toplevel
>
> while (count <= botlevel)
>
> 'set lev 'count
>  'define th1 = temp * (pow(100000 / press , (0.286)))'
> 'define u1 = u'
>  'define v1 = v'
> 'define z1 = geopot/9.81'
>
> 'set lev 'count - 1
>  'define th2 = temp * (pow(100000 / press , (0.286)))'
> 'define u2 = u'
>  'define v2 = v'
> 'define z2 = geopot/9.81'
>
> 'define
> ri=9.81/(0.5*th1+th2)*(th2-th1)/(pow((u2-u1),2)+pow((v2-v1),2))*(z2-z1)'
>
> 'd 'ri
> say result *THIS SEEMS TO WORK AS I GET 10 DIFFERENT VALUES OF RI
>
> if (count = toplevel) *THESE IF-LOOPS AREN'T OKAY, BUT WHY?
> minri = ri
>  minlevel = count
> endif
>
> if (ri < minri)
>  minri = ri
> minlevel = count
> endif
>
> count = count+1
> endwhile
>
> 'd 'minri
> say "Min Ri " result *ALWAYS THE VALUE OF RI AT LEVEL 50 (AT THE BOTTOM
> LEVEL)
> 'd 'minlevel
> say "Model level " result *ALWAYS LEVEL 40 (TOP LEVEL)
>
>
>
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
>


-- 
Jeff Duda
Iowa State University
Meteorology Graduate Student
www.meteor.iastate.edu/~jdduda <http://www.meteor.iastate.edu/%7Ejdduda>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20111103/26251ce0/attachment-0003.html