topography plot from WRF output

[Michael J Murphy]

The script in the previous reply was very useful, thank you Javier.

I am still having an issue of blank data around the mountains in my cross sections.  I tried searching the archives for a previous post about this but did not find anything.

Using Javier's script I have been able to get a nice plot of the terrain and my chosen meteorological variables in a cross section but the problem is the data cuts off and is blank near to the mountains.  This is obviously caused by the fact that the model output is interpolated into a finite number of pressure levels.  I have tried dramatically increasing the number of pressure levels the data is written onto and this helps to make the blank data  areas smaller but they are still quite visible.

Does anyone have any experience with this kind of problem?  I know it sounds weird but, maybe there is some way to plot the topography in a cross section with sigma levels instead of pressure levels?

 Thanks,
MJ


----- Original Message -----
From: "Javier G. Corripio" <Javier.Corripio at UIBK.AC.AT>
Date: Thursday, August 28, 2008 12:14 am
Subject: Re: topography plot from WRF output
To: GRADSUSR at LIST.CINECA.IT

> This may help:



> The script is for my own needs, so it does not check for any exceptions.   If you prefer to use hgtsfc converted to pressure instead of pressfc use the commented lines at the beginning.   The script displays the model surface height as a grey mask across a transect with pressure values in the Y-axis.

> Best,

> Javier

> *plotsfc.gs  Javier.Corripio at uibk.ac.at

> *'define earthg=9.80665'     ;*acceleration due to gravity (m s-2)
> *'define md=28.966'           ;*molecular weight of dry air
> *'define atmosr=8.31432'     ;*universal gas constant (j deg-1 kmol-1)
> *'define p0=prmslmsl(lev=1013)'
> *'define T0=tmpprs(lev=1013)'
> *'define lrate=(tmpprs-tmpprs(lev=1013))/hgtsfc(lev=1013)'
> *'define base=(T0/(T0+lrate*hgtsfc))'
> *'define exponent=((earthg*md)/(atmosr*lrate*1000))'
> *'define pz=p0*POW(base,exponent)'

> *get gx info
> 'q gxinfo'
> data=sublin(result,3)
> px0 = subwrd(data,4)
> px1 = subwrd(data,6)
> data=sublin(result,4)
> py0 = subwrd(data,4)
> py1 = subwrd(data,6)

> say px0 ' ' px1 ' ' py0 ' ' py1

> *get dim info
> 'q dim'
> data=sublin(result,4)
> lev0 = subwrd(data,6)
> lev1 = subwrd(data,8)

> say lev0 ' ' lev1

> 'set lev ' lev0

> *set plotting area
> 'set parea ' px0 ' 'px1 ' ' py0 ' ' py1

> *set Yaxis range
> 'set vrange ' lev0 ' ' lev1

> 'set gxout linefill'
> 'set lfcols 15 15'
> 'set cmark 0'
> 'd pressfc/100;const(pressfc/100,' lev0 ')'

> 'set gxout line'
> 'set parea off'
> 'set lev ' lev0 ' ' lev1



> On 28 Aug 2008, at 08:18, Michael J Murphy wrote:

> Grads users,

> I am trying to plot topography in a cross section with pressure as the vertical coordinate.  I have read a lot about this in the user archives and even found a script (orog.gs) that was written to do it. 
 
This script uses the variable psfc (press. at the surface).  Unfortunately the variables that could be used in the script that are available from my WRF-ARW model output (terrain height, psfc, slp) are only on the lowest level.  Because of this when try:

> set lat 21.5
> set lev 1000 700
> d psfc

> The result is nothing being visibly plotted (because psfc is only on the lowest level).  Now if I set the vertical axis to a single pressure value like so:

> set lat 21.5
> set lev 1000
> d psfc

> The result is an x-y graph of the surface pressure with longitude on the x-axis and pressure on the y-axis.  The result is the similar when I try to plot terrain height and slp as well.  

> Because psfc is only on the lowest level I cannot do the  var = lev-psfc/100  command.  What I really need is terrain height as a function of pressure!!
> Maybe someone has run into this problem with WRF output data and has figure this out?  Maybe I would be better off using a different program such as RIP?  Any advice will be greatly appreciated.

> Thanks,
> MJ
 


> -- 
> ***********************************************
> Javier G. Corripio
> Faculty of Geo- and Atmospheric Sciences
> University of Innsbruck
> Innrain, 52
> A-6020 Innsbruck
> Austria
> Tel.: +43 512 507 5418
> Fax:  +43 512 507 2895
> email: Javier.Corripio at uibk.ac.at
> url: http://www.uibk.ac.at/geographie/personal/corripio/
> http://web.mac.com/jgcorripio/
> http://www.meteoexploration.com





-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20080828/17fb8c29/attachment.html