Friday, May 17, 2013

Script:; Defines potential and equlivalent potential temperature given temperature and moisture inputs

This script, is really more of a function, the best way to use this would be to copy the contents of this script into the bottom of the script you wish to use it in.  This function reads in temperature and moisture information and returns a good approximation of both potential and equivalent potential temperature.

This function asks for up to four inputs, though only temperature and moisture inputs are required.
It is capable of calculating  equivalent potential temperature given four different moisture input variables:
  1. Relative Humidity (-rh)
  2. Mixing Ratio (-r)
  3. Vapor Pressure (-e)
  4. Specific Humidity (-q)
The default is relative humidity.

Example Usage: theta(tmpprs,rhprs)
                           -This example would define PT and EPT using relative humidity as the input.

In addition to calculating potential and equivalent potential temperature at given pressure levels, this function also allows you to specify surface pressure as the pressure input.  This adjusts for the fact that it would simply calculate everything at fixed pressure levels.

Example Usage: theta(tmp2m,spfhum2m,pressfc,'-q')
                           -This returns the surface PT and EPT calculated using specific humidity as the input.

Notes on Usage: As always, all input variables must be in MKS units.  If you wish to use a method other than relative humidity, but do not wish to use a specific variable for surface pressure, just type 'lev' or '-' in place of where surface pressure would normally go (e.g., theta(tmp,mixr,'-','-r'))  Be sure to put '' around any variables with a "-" in them, or GrADS will try to subtract.

I have not tested this function with all of the above variables,but this function does work well with both relative humidity and specific humidity.  As always, report any bugs here! 

*Minor bug fixed 09/2013



  1. Hi, thanks for the script. I have temperature and relative humidity data (1000 mb to 100 mb pressure levels). With only these two inputs may I expect your function to return me EPT values ?


  2. Please comment on the following script ... IS IT OK ?


    'open TMP.ctl'
    'open RH.ctl'

    'set fwrite THETAE.grd'
    'set undef -999.9'

    say "tt="tt

    say "zz="zz

    'set x 1 144'
    'set y 1 73'
    'set z 'zz
    'set t 'tt
    'q dim'
    say 'press=' press
    'define es=6.1173*exp(((2.501*pow(10,6))/461.50)*(1/273.16 - 1/tmp.1))'
    'define ws=621.97*(es/(lev-es))'
    'define w=(rh.2*ws)/(100*1000)'
    'define thetae=(tmp+((2.501*pow(10,6))/1004)*w)*pow((1000/'press'),287/1004)'

    'set gxout fwrite'
    'd thetae'



    'disable fwrite'


    1. It looks like you are trying to write out a GrADS binary file, and in terms of how you are writing out your definitions of thetaE it looks pretty darn good to me. I don't see any glaring errors, or anything wrong with what you are doing, I also don't have access to GrADS right now, so I cannot 100% confirm that this is good, but it looks pretty correct to me!

  3. Thanks "Anonymous" for your nice comments ... BTW i plotted the output and the theta-e profile looks OK ... is there any way that I can attach the plot here?