[gradsusr] Lifted Index Formula
Mike Manning
michael at bsch.au.com
Fri Mar 2 21:00:59 EST 2012
Hi Charles,
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:
Starting "/grads2/Linux/Versions/2.0.a9.oga.1/i686/grads -blc test.gs
" ...
Grid Analysis and Display System (GrADS) Version 2.0.a9.oga.1
Copyright (c) 1988-2010 by Brian Doty and the
Institute for Global Environment and Society (IGES)
GrADS comes with ABSOLUTELY NO WARRANTY
See file COPYRIGHT for more information
Config: v2.0.a9.oga.1 little-endian readline printim grib2 netcdf
hdf4-sds hdf5 opendap-grids,stn athena geotiff shapefile
Issue 'q config' command for more detailed configuration information
Loading User Defined Extensions table
</grads2/Linux/Versions/2.0.a9.oga.1/i686/gex/udxt> ... ok.
GX Package Initialization: Size = 11 8.5
Running in Batch mode
No hardcopy metafile open
All files closed; all defined objects released;
All GrADS attributes have been reinitialized
Non-numeric args to numeric operation
Error occurred processing function arguments
Error occurred on line 198
In file test.gs
Error occurred on line 62
In file test.gs
ga-> quit
No hardcopy metafile open
GX package terminated
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:
Non-numeric args to numeric operation
Error occurred on line 62
Here's the code I am using now, I've highlighted the lines in bold that
it's showing the errors with:
'c'
'reinit'
'open
/home/subdomains/forecasts/public_html/grib/grib2_cmc/2012030212/gfs.12z.ctl'
'set rgb 50 0 128 0'
'set xlab off'
'set ylab off'
'set frame off'
'set grads off'
'set vpage 0 11 0 8.5'
'set parea 0 11 0 8.5'
'set mproj scaled'
'set mpdraw off'
'set map 0 1 5'
'set csmooth on'
'set t 1'
'set lon 109 159'
'set lat -48 -5'
'set rgb 50 245 201 255'
'set rgb 51 234 144 255'
'set rgb 52 209 93 236'
'set rgb 53 175 62 202'
'set rgb 54 137 7 168'
'set rgb 55 5 83 206'
'set rgb 70 8 74 179'
'set rgb 56 30 110 235'
'set rgb 57 60 150 245'
'set rgb 58 150 210 250'
'set rgb 59 180 235 255'
'set rgb 60 255 46 46'
'set rgb 61 201 36 36'
'set rgb 62 255 255 255'
'set rgb 66 254 254 254'
'set rgb 99 255 0 0'
'set clevs -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0'
'set ccols 50 50 51 52 53 54 70 55 56 57 58 59 62'
'define tsfc = TMP2m-273.14'
'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))'
'define psfc = PRESsfc/100'
'set lev 1000'
'define t1000mb = TMPprs-273.14'
'set lev 925'
'define t925mb = TMPprs-273.14'
'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))'
'set lev 850'
'define t850mb = TMPprs-273.14'
'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))'
'set lev 700'
'define t700mb = TMPprs-273.14'
'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))'
'set lev 500'
'define t500mb = TMPprs-273.14'
'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))'
*'define Temp500V='virtual2('t500mb+273.15','dp500mb+273.15',500)'-273.15'*
'define TLcl='Templcl('tsfc','dpsfc')
'define PLcl='Preslcl('tsfc','dpsfc','psfc')
'define PclTemp='LiftWet('TLcl','PLcl',500)
'define PclTempV='virtual2('PclTemp+273.15','PclTemp+273.15',500)'-273.15'
'define LIMax=Temp500V-PclTempV'
'define variable=LIMax'
'set gxout shaded'
'set csmooth on'
'set cint 1'
'd variable'
'set gxout contour'
'set cstyle 3'
'set ccolor 0'
'set cint 5'
'd variable'
'printim li.png x750 y645 black'
'quit'
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)
function gammaw(tempc,pres,rh)
*-----------------------------------------------------------------------
* Function to calculate the moist adiabatic lapse rate (deg C/Pa) based
* on the temperature, pressure, and rh of the environment.
*----------------------------------------------------------------------
tempk=tempc+273.15
es=satvap2(tempc)
ws=mixratio(es,pres)
w=rh*ws/100
tempv=virtual(tempk,w)
latent=latentc(tempc)
A=1.0+latent*ws/(287*tempk)
B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk)
Density=100*pres/(287*tempv)
lapse=(A/B)/(1005*Density)
return(lapse)
function latentc(tempc)
*-----------------------------------------------------------------------
* Function to return the latent heat of condensation in J/kg given
* temperature in degrees Celsius.
*-----------------------------------------------------------------------
val=2502.2-2.43089*tempc*1000
return(val)
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)
function Preslcl(temp,dewp,pres)
*-------------------------------------------------------
* Calculate press of LCL given temp & dewp in C and pressure
*-------------------------------------------------------
Tlclk=Templcl(temp,dewp)+273.15
tempk=temp+273.15
theta=tempk*pow(1000/pres,0.286)
plcl=1000*pow(Tlclk/theta,3.48)
return(plcl)
function mixratio(e,p)
*------------------------------------------------------
* Given vapor pressure and pressure, return mixing ratio
*-------------------------------------------------------
mix=0.622*e/(p-e)
return(mix)
function virtual2(temp,dewp,pres)
*------------------------------------------------------------
* Function to return virtual temperature in kelvin given temperature in
* kelvin and dewpoint in kelvin and pressure in mb
*-------------------------------------------------------------
if (temp > 0 & dewp > 0)
* vap=satvap2(dewp-273.15)*
mix=mixratio(vap,pres)
tempv=virtual(temp,mix)
else
tempv=-9999
endif
return (tempv)
function virtual(temp,mix)
*------------------------------------------------------------
* Function to return virtual temperature given temperature in
* kelvin and mixing ratio in g/g.
*-------------------------------------------------------------
tempv=temp*(1.0+0.6*mix)
return (tempv)
function mixratio(e,p)
*------------------------------------------------------
* Given vapor pressure and pressure, return mixing ratio
*-------------------------------------------------------
mix=0.622*e/(p-e)
return(mix)
function satvap2(temp)
*---------------------------------------------------------------
* Given temp in Celsius, returns saturation vapor pressure in mb
*---------------------------------------------------------------
es=6.112*exp(17.67*temp/(temp+243.5))
return(es)
On 3/03/2012 4:22 AM, Charles Seman wrote:
> 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...
> http://grads.iges.org/grads/gadoc/mathfunctions.html
>
> 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:
>> 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
>> (http://moe.met.fsu.edu/~rhart/software/plotskew.gs) 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:
>>>
>>> On Thu, Mar 1, 2012 at 11:25 AM, Jeff Duda<jeffduda319 at gmail.com
>>> <mailto:jeffduda319 at gmail.com>> 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
>>> <http://opengrads.org/doc/udxt/gsudf/>. 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<michael at bsch.au.com
>>> <mailto:michael at bsch.au.com>> 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<http://plotskew.gs> 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
>>> gradsusr at gradsusr.org<mailto:gradsusr at gradsusr.org>
>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>
>>>
>>>
>>>
>>> --
>>> Jeff Duda
>>> Graduate research assistant
>>> University of Oklahoma School of Meteorology
>>> Center for Analysis and Prediction of Storms
>>>
>>>
>>> _______________________________________________
>>> gradsusr mailing list
>>> gradsusr at gradsusr.org<mailto:gradsusr at gradsusr.org>
>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>>
>>>
>>>
>>>
>>> --
>>> Arlindo da Silva
>>> dasilva at alum.mit.edu<mailto:dasilva at alum.mit.edu>
>>>
>>>
>>> _______________________________________________
>>> gradsusr mailing list
>>> gradsusr at gradsusr.org
>>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>>
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org
>> http://gradsusr.org/mailman/listinfo/gradsusr
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20120303/3bc82451/attachment-0003.html
More information about the gradsusr
mailing list