<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Hi Charles,<br>
<br>
Yeah I originally had script variables in my .gs gradsscript, but
when I executed it (grads -blc scriptname.gs) it kept throwing up
errors all the time? I originally just copied and pasted all the
functions exactly as is, this is the output I get:<br>
<br>
Starting "/grads2/Linux/Versions/2.0.a9.oga.1/i686/grads -blc
test.gs " ...<br>
<br>
<br>
Grid Analysis and Display System (GrADS) Version 2.0.a9.oga.1<br>
Copyright (c) 1988-2010 by Brian Doty and the<br>
Institute for Global Environment and Society (IGES)<br>
GrADS comes with ABSOLUTELY NO WARRANTY<br>
See file COPYRIGHT for more information<br>
<br>
Config: v2.0.a9.oga.1 little-endian readline printim grib2 netcdf
hdf4-sds hdf5 opendap-grids,stn athena geotiff shapefile<br>
Issue 'q config' command for more detailed configuration information<br>
Loading User Defined Extensions table
</grads2/Linux/Versions/2.0.a9.oga.1/i686/gex/udxt> ... ok.<br>
GX Package Initialization: Size = 11 8.5<br>
Running in Batch mode<br>
No hardcopy metafile open<br>
All files closed; all defined objects released;<br>
All GrADS attributes have been reinitialized<br>
Non-numeric args to numeric operation<br>
Error occurred processing function arguments<br>
Error occurred on line 198<br>
In file test.gs<br>
Error occurred on line 62<br>
In file test.gs<br>
ga-> quit<br>
No hardcopy metafile open<br>
GX package terminated<br>
<br>
Even if I change line 62 to say
Temp500V=virtual2(t500mb+273.15,dp500mb+273.15,,500)-273.15 (or even
something simple like t500mb = TMPprs-273.14) it still shows errors:
<br>
Non-numeric args to numeric operation<br>
Error occurred on line 62<br>
<br>
<br>
Here's the code I am using now, I've highlighted the lines in bold
that it's showing the errors with:<br>
<br>
<br>
'c'<br>
'reinit'<br>
'open
/home/subdomains/forecasts/public_html/grib/grib2_cmc/2012030212/gfs.12z.ctl'<br>
'set rgb 50 0 128 0'<br>
'set xlab off'<br>
'set ylab off'<br>
'set frame off'<br>
'set grads off'<br>
'set vpage 0 11 0 8.5'<br>
'set parea 0 11 0 8.5'<br>
'set mproj scaled'<br>
'set mpdraw off'<br>
'set map 0 1 5'<br>
'set csmooth on'<br>
'set t 1'<br>
'set lon 109 159'<br>
'set lat -48 -5'<br>
<br>
'set rgb 50 245 201 255'<br>
'set rgb 51 234 144 255'<br>
'set rgb 52 209 93 236'<br>
'set rgb 53 175 62 202'<br>
'set rgb 54 137 7 168'<br>
'set rgb 55 5 83 206'<br>
'set rgb 70 8 74 179'<br>
'set rgb 56 30 110 235'<br>
'set rgb 57 60 150 245'<br>
'set rgb 58 150 210 250'<br>
'set rgb 59 180 235 255'<br>
'set rgb 60 255 46 46'<br>
'set rgb 61 201 36 36'<br>
'set rgb 62 255 255 255'<br>
'set rgb 66 254 254 254'<br>
'set rgb 99 255 0 0'<br>
<br>
'set clevs -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0'<br>
'set ccols 50 50 51 52 53 54 70 55 56 57 58 59 62'<br>
<br>
'define tsfc = TMP2m-273.14'<br>
'define dpsfc =
tsfc-((14.55+0.114*tsfc)*(1-0.01*RHprs)+pow((2.5+0.007*tsfc)*(1-0.01*RHprs),3)+(15.9+0.117*tsfc)*pow((1-0.01*RHprs),14))'<br>
'define psfc = PRESsfc/100'<br>
<br>
'set lev 1000'<br>
'define t1000mb = TMPprs-273.14'<br>
<br>
'set lev 925'<br>
'define t925mb = TMPprs-273.14'<br>
'define dp925mb =
t925mb-((14.55+0.114*t925mb)*(1-0.01*RHprs)+pow((2.5+0.007*t925mb)*(1-0.01*RHprs),3)+(15.9+0.117*t925mb)*pow((1-0.01*RHprs),14))'<br>
<br>
'set lev 850'<br>
'define t850mb = TMPprs-273.14'<br>
'define dp850mb =
t850mb-((14.55+0.114*t850mb)*(1-0.01*RHprs)+pow((2.5+0.007*t850mb)*(1-0.01*RHprs),3)+(15.9+0.117*t850mb)*pow((1-0.01*RHprs),14))'<br>
<br>
'set lev 700'<br>
'define t700mb = TMPprs-273.14'<br>
'define dp700mb =
t700mb-((14.55+0.114*t700mb)*(1-0.01*RHprs)+pow((2.5+0.007*t700mb)*(1-0.01*RHprs),3)+(15.9+0.117*t700mb)*pow((1-0.01*RHprs),14))'<br>
<br>
'set lev 500'<br>
'define t500mb = TMPprs-273.14'<br>
'define dp500mb =
t500mb-((14.55+0.114*t500mb)*(1-0.01*RHprs)+pow((2.5+0.007*t500mb)*(1-0.01*RHprs),3)+(15.9+0.117*t500mb)*pow((1-0.01*RHprs),14))'<br>
<br>
<font color="#ff0000"><b>'define
Temp500V='virtual2('t500mb+273.15','dp500mb+273.15',500)'-273.15'</b></font><br>
'define TLcl='Templcl('tsfc','dpsfc')<br>
'define PLcl='Preslcl('tsfc','dpsfc','psfc')<br>
<br>
'define PclTemp='LiftWet('TLcl','PLcl',500)<br>
'define
PclTempV='virtual2('PclTemp+273.15','PclTemp+273.15',500)'-273.15'<br>
'define LIMax=Temp500V-PclTempV'<br>
<br>
'define variable=LIMax'<br>
<br>
'set gxout shaded'<br>
'set csmooth on'<br>
'set cint 1'<br>
'd variable'<br>
<br>
'set gxout contour'<br>
'set cstyle 3'<br>
'set ccolor 0'<br>
'set cint 5'<br>
'd variable'<br>
<br>
'printim li.png x750 y645 black'<br>
'quit'<br>
<br>
<br>
<br>
<br>
function LiftWet(startt,startp,endp)<br>
<br>
*--------------------------------------------------------------------<br>
* Lift a parcel moist adiabatically from startp to endp.<br>
* Init temp is startt in C. If you wish to see the parcel's<br>
* path plotted, display should be 1. Returns temp of parcel at
endp.<br>
*--------------------------------------------------------------------<br>
<br>
temp=startt<br>
pres=startp<br>
delp=10<br>
While (pres >= endp)<br>
temp=temp-100*delp*gammaw(temp,pres-delp/2,100)<br>
pres=pres-delp<br>
EndWhile<br>
return(temp)<br>
<br>
<br>
<br>
<br>
function gammaw(tempc,pres,rh)<br>
<br>
*-----------------------------------------------------------------------<br>
* Function to calculate the moist adiabatic lapse rate (deg C/Pa)
based<br>
* on the temperature, pressure, and rh of the environment.<br>
*----------------------------------------------------------------------<br>
<br>
tempk=tempc+273.15<br>
es=satvap2(tempc)<br>
ws=mixratio(es,pres)<br>
w=rh*ws/100<br>
tempv=virtual(tempk,w)<br>
latent=latentc(tempc)<br>
<br>
A=1.0+latent*ws/(287*tempk)<br>
B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk)<br>
Density=100*pres/(287*tempv)<br>
lapse=(A/B)/(1005*Density)<br>
return(lapse)<br>
<br>
<br>
<br>
<br>
function latentc(tempc)<br>
<br>
*-----------------------------------------------------------------------<br>
* Function to return the latent heat of condensation in J/kg given<br>
* temperature in degrees Celsius.<br>
*-----------------------------------------------------------------------<br>
<br>
val=2502.2-2.43089*tempc*1000<br>
<br>
return(val)<br>
<br>
<br>
<br>
<br>
function Templcl(temp,dewp)<br>
<br>
*------------------------------------------------------<br>
* Calculate the temp at the LCL given temp & dewp in C<br>
*------------------------------------------------------<br>
<br>
tempk=temp+273.15<br>
dewpk=dewp+273.15<br>
Parta=1/(dewpk-56)<br>
Partb=log(tempk/dewpk)/800<br>
Tlcl=(1/(Parta+Partb)+56)-273.15<br>
return(Tlcl)<br>
<br>
<br>
<br>
<br>
function Preslcl(temp,dewp,pres)<br>
<br>
*-------------------------------------------------------<br>
* Calculate press of LCL given temp & dewp in C and pressure<br>
*-------------------------------------------------------<br>
<br>
Tlclk=Templcl(temp,dewp)+273.15<br>
tempk=temp+273.15<br>
theta=tempk*pow(1000/pres,0.286)<br>
plcl=1000*pow(Tlclk/theta,3.48)<br>
return(plcl)<br>
<br>
<br>
<br>
<br>
function mixratio(e,p)<br>
<br>
*------------------------------------------------------<br>
* Given vapor pressure and pressure, return mixing ratio<br>
*-------------------------------------------------------<br>
<br>
mix=0.622*e/(p-e)<br>
<br>
return(mix)<br>
<br>
<br>
<br>
<br>
function virtual2(temp,dewp,pres)<br>
<br>
*------------------------------------------------------------<br>
* Function to return virtual temperature in kelvin given temperature
in<br>
* kelvin and dewpoint in kelvin and pressure in mb<br>
*-------------------------------------------------------------<br>
<br>
if (temp > 0 & dewp > 0)<br>
<font color="#ff0000"><b> vap=satvap2(dewp-273.15)</b></font><br>
mix=mixratio(vap,pres)<br>
tempv=virtual(temp,mix)<br>
else<br>
tempv=-9999<br>
endif<br>
<br>
return (tempv)<br>
<br>
<br>
<br>
<br>
function virtual(temp,mix)<br>
<br>
*------------------------------------------------------------<br>
* Function to return virtual temperature given temperature in<br>
* kelvin and mixing ratio in g/g.<br>
*-------------------------------------------------------------<br>
<br>
tempv=temp*(1.0+0.6*mix)<br>
<br>
return (tempv)<br>
<br>
<br>
<br>
<br>
function mixratio(e,p)<br>
<br>
*------------------------------------------------------<br>
* Given vapor pressure and pressure, return mixing ratio<br>
*-------------------------------------------------------<br>
<br>
mix=0.622*e/(p-e)<br>
<br>
return(mix)<br>
<br>
<br>
<br>
<br>
function satvap2(temp)<br>
<br>
*---------------------------------------------------------------<br>
* Given temp in Celsius, returns saturation vapor pressure in mb<br>
*---------------------------------------------------------------<br>
<br>
es=6.112*exp(17.67*temp/(temp+243.5))<br>
<br>
return(es)<br>
<br>
<br>
<br>
On 3/03/2012 4:22 AM, Charles Seman wrote:
<blockquote cite="mid:4F510FF7.4010000@noaa.gov" type="cite">
<pre wrap="">Mike,
I took a look at a local version of Bob Hart's plotskew.gs script, and
found his "LiftWet" function:
**************************************************************************
function LiftWet(startt,startp,endp,display,Pmin,Pmax)
*--------------------------------------------------------------------
* Lift a parcel moist adiabatically from startp to endp.
* Init temp is startt in C. If you wish to see the parcel's
* path plotted, display should be 1. Returns temp of parcel at endp.
*--------------------------------------------------------------------
temp=startt
pres=startp
cont = 1
delp=10
While (pres >= endp & cont = 1)
If (display = 1)
xtemp=GetXLoc(temp,pres)
"q w2xy "xtemp" "pres
xloc=subwrd(result,3)
yloc=subwrd(result,6)
If (xtemp < 0 | xtemp > 100)
cont=0
Else
If (pres >= Pmin & pres < Pmax & pres < startp)
"draw line "xold" "yold" "xloc" "yloc
Endif
Endif
Endif
xold=xloc
yold=yloc
temp=temp-100*delp*gammaw(temp,pres-delp/2,100)
pres=pres-delp
EndWhile
return(temp)
**************************************************************************
What we see in this function (and I believe in much of plotskew.gs) is
the use of GrADS "script" variables, and not GrADS "defined"
variables... In plotskew.gs, GrADS "defined" variables are used for
variables which are plotted, like in this code:
------------------------------------------------------------
*-----------------------------
* Draw wind profile along side
*-----------------------------
...
zz=1
wspd=-999
cont=1
While (cont = 1 & zz < _zmax)
"set z "zz
pres=subwrd(result,4)
"define dd="sndspd
"d dd"
rec=sublin(result,8)
wspd=subwrd(rec,4)
if (math_abs(wspd) > 500 | pres > _pmax)
zz=zz+1
else
cont=0
Endif
Endwhile
------------------------------------------------------------
The above code is a nice example of how to use both GrADS "script" and
"defined" variables (and how to make "defined" variables using "script"
variables: "define dd="sndspd), and using script math functions...
<a class="moz-txt-link-freetext" href="http://grads.iges.org/grads/gadoc/mathfunctions.html">http://grads.iges.org/grads/gadoc/mathfunctions.html</a>
Could you recode your script to use "script" variables and not "defined"
variables? Looking at
'temp='startt
'pres='startp
'delp=10'
While (pres >= endp)
'temp=temp-100*delp*'gammaw('temp','pres-delp/2',100)
'pres=pres-delp'
EndWhile
one thing to note is the variable "pres" is a "defined" variable, but
used as a "script" variable in the "While()" statement...
Hope this helps,
Chuck
On 03/02/2012 08:23 AM, Mike Manning wrote:
</pre>
<blockquote type="cite">
<pre wrap="">Cool will keep this in mind as well - I had been making good progress
with everything until I hit the function for LiftWet. On plotskew.gs
(<a class="moz-txt-link-freetext" href="http://moe.met.fsu.edu/~rhart/software/plotskew.gs">http://moe.met.fsu.edu/~rhart/software/plotskew.gs</a>) everything there is
taken for a single point on a chart - I guess because I'm trying to use
an entire chart section, it's why I'm having dramas?
I've got all other functions working perfectly.. here's one for example:
function Templcl(temp,dewp)
*------------------------------------------------------
* Calculate the temp at the LCL given temp & dewp in C
*------------------------------------------------------
'tempk='temp'+273.15'
'dewpk='dewp'+273.15'
'Parta=1/(dewpk-56)'
'Partb=log(tempk/dewpk)/800'
'Tlcl=(1/(Parta+Partb)+56)-273.15'
return(Tlcl)
This is the function I am having trouble with, it basically gets stuck
in the loop and never exits :(
function LiftWet(startt,startp,endp)
*--------------------------------------------------------------------
* Lift a parcel moist adiabatically from startp to endp.
* Init temp is startt in C. If you wish to see the parcel's
* path plotted, display should be 1. Returns temp of parcel at endp.
*--------------------------------------------------------------------
'temp='startt
'pres='startp
'delp=10'
While (pres >= endp)
'temp=temp-100*delp*'gammaw('temp','pres-delp/2',100)
'pres=pres-delp'
EndWhile
return(temp)
I call it by 'define PclTemp='LiftWet('TLcl','PLcl',500) where TLcl and
PLcl work perfectly when I plot those 2 variables by themselves so I
know they are all good. I just have a funny feeling now that I have hit
a brick wall that I cannot get around - which probably explains why no
one else has created a Lifted Index chart. I do hope I am wrong though
Cheers, Mike
On 2/03/2012 11:24 AM, Arlindo da Silva wrote:
</pre>
<blockquote type="cite">
<pre wrap="">
On Thu, Mar 1, 2012 at 11:25 AM, Jeff Duda <<a class="moz-txt-link-abbreviated" href="mailto:jeffduda319@gmail.com">jeffduda319@gmail.com</a>
<a class="moz-txt-link-rfc2396E" href="mailto:jeffduda319@gmail.com"><mailto:jeffduda319@gmail.com></a>> wrote:
You've confused the different uses of scripting functions vs.
Grads prompt functions. You can only display a variable that you
have defined using the 'define ...' command (note the use of the
quotes signaling that it was defined as if you typed it in at the
Grads prompt). Your use of variable = satvap(...) with satvap
being a script function only is causing Grads to assume variable
is a script variable only, not one that can be displayed. The
reason you get the other error when you put the other line in
quotes is because math_exp is strictly a script version of the exp
function that you would use in the Grads command prompt. Thus, you
would want to change that line to
'es = 6.112*exp(17.67*temp/(temp+243.5))'
in order to display es or define any other variables that would
depend on it.
This is generally true but if you are using opengrads you can make the
script function satvap() available to the display command, say
ga-> d satvap(t1000)
by putting it in its own file "satvap.gsf" and loading the
corresponding UDXT as described in the gsUDF documentation
<a class="moz-txt-link-rfc2396E" href="http://opengrads.org/doc/udxt/gsudf/"><http://opengrads.org/doc/udxt/gsudf/></a>. Actually, you could in
principle implemented the lifted index as a gsUDF function.
Arlindo
Jeff Duda
On Thu, Mar 1, 2012 at 7:06 AM, Mike Manning <<a class="moz-txt-link-abbreviated" href="mailto:michael@bsch.au.com">michael@bsch.au.com</a>
<a class="moz-txt-link-rfc2396E" href="mailto:michael@bsch.au.com"><mailto:michael@bsch.au.com></a>> wrote:
Hi everyone,
I'm using data from the FNMOC model and am trying to make a
function
that calculates Lifted Index. I'm using bits of the code from the
plotskew.gs <a class="moz-txt-link-rfc2396E" href="http://plotskew.gs"><http://plotskew.gs></a> code. So far I've got the
temperature (celsius) and
dewpoint (celsius) worked out for each level. Now I'm
calculating the
parcel details.. I've pulled some of the functions and have
just been
testing to make sure I'm on the right track. If I have this
code for
example:
'c'
'reinit'
'open gfs.00z.ctl'
'set t 7'
'set lon 135 155'
'set lat -31 -10'
'set lev 1000'
't1000mb = TMPprs-273.14'
'define dp1000mb =
t1000mb-((14.55+0.114*t1000mb)*(1-0.01*RHprs)+pow((2.5+0.007*t1000mb)*(1-0.01*RHprs),3)+(15.9+0.117*t1000mb)*pow((1-0.01*RHprs),14))'
variable = satvap2(t1000mb)
'd variable'
function satvap2(temp)
*---------------------------------------------------------------
* Given temp in Celsius, returns saturation vapor pressure in mb
*---------------------------------------------------------------
es=6.112*math_exp(17.67*temp/(temp+243.5))
return(es)
'quit'
it complains about an error? If I put the code in the satvap2
function
in single quotes it then says "math_exp" is not a variable or
function"?
I'm a little lost on how to fix this one up if it's possible.
Cheers, Mike
_______________________________________________
gradsusr mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gradsusr@gradsusr.org">gradsusr@gradsusr.org</a> <a class="moz-txt-link-rfc2396E" href="mailto:gradsusr@gradsusr.org"><mailto:gradsusr@gradsusr.org></a>
<a class="moz-txt-link-freetext" href="http://gradsusr.org/mailman/listinfo/gradsusr">http://gradsusr.org/mailman/listinfo/gradsusr</a>
--
Jeff Duda
Graduate research assistant
University of Oklahoma School of Meteorology
Center for Analysis and Prediction of Storms
_______________________________________________
gradsusr mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gradsusr@gradsusr.org">gradsusr@gradsusr.org</a> <a class="moz-txt-link-rfc2396E" href="mailto:gradsusr@gradsusr.org"><mailto:gradsusr@gradsusr.org></a>
<a class="moz-txt-link-freetext" href="http://gradsusr.org/mailman/listinfo/gradsusr">http://gradsusr.org/mailman/listinfo/gradsusr</a>
--
Arlindo da Silva
<a class="moz-txt-link-abbreviated" href="mailto:dasilva@alum.mit.edu">dasilva@alum.mit.edu</a> <a class="moz-txt-link-rfc2396E" href="mailto:dasilva@alum.mit.edu"><mailto:dasilva@alum.mit.edu></a>
_______________________________________________
gradsusr mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gradsusr@gradsusr.org">gradsusr@gradsusr.org</a>
<a class="moz-txt-link-freetext" href="http://gradsusr.org/mailman/listinfo/gradsusr">http://gradsusr.org/mailman/listinfo/gradsusr</a>
</pre>
</blockquote>
<pre wrap="">
_______________________________________________
gradsusr mailing list
<a class="moz-txt-link-abbreviated" href="mailto:gradsusr@gradsusr.org">gradsusr@gradsusr.org</a>
<a class="moz-txt-link-freetext" href="http://gradsusr.org/mailman/listinfo/gradsusr">http://gradsusr.org/mailman/listinfo/gradsusr</a>
</pre>
</blockquote>
<pre wrap="">
</pre>
</blockquote>
<br>
</body>
</html>