[gradsusr] Lifted Index Formula
Charles Seman
charles.seman at noaa.gov
Mon Mar 5 12:24:49 EST 2012
Hi Matt,
Sorry, but I'm not able to help much with this. I don't know what
exactly is preventing the script variables version from working. One
thing to note is that with script variables, in order to use a dataset
or local GrADS-defined variable, you have to specify a single (x,y,z,t)
point, do a "d somevar", and then get the value of the "result" using
something like:
'set x 'xvalue
'set y 'yvalue
'set z 'zvalue
'set t 'tvalue
'd somevar'
somevarvalue = subwrd(result,4) <- this usually works without using
"sublin" first but depends on what "gxout" is set to
* note "somevar" could be a dataset variable or some variable you
* defined locally in GrADS
I'm assuming that the variables like "TMPprs" are from the FNMOC model
output dataset? (Don't know what's in your dataset.)
Hope this helps,
Chuck
On 03/02/2012 09:00 PM, Mike Manning wrote:
> 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
>
>
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
--
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 reflect any
official or unofficial position of the United States Federal Government,
the United States Department of Commerce, or NOAA."
More information about the gradsusr
mailing list