<div>Dear Robert,</div>
<div> </div>
<div>I thank you for your suggestion.</div>
<div> </div>
<div>I did it the way shown below, but it did not work either!</div>
<div>Any other suggestion?</div>
<div> </div>
<div>Thanks in advance.</div>
<div>Eduardo.</div>
<div>__</div>
<div> </div>
<div>'reinit'</div>
<p>'sdfopen d:\ncep-ncar\uwnd.mon.mean.nc'<br>'sdfopen d:\ncep-ncar\air.mon.mean.nc'<br></p>
<p>'set lon -180 0' <br>'set lat 0 90' <br>'set lev 1050 10' <br>'set t 1' </p>
<p>'define uisen='isen(uwnd.1,air.2,lev,320)</p>
<p>"set z 1"<br>'set t 1'</p>
<p>"d uisen"<br></p>
<div>___</div>
<div> </div>
<div> </div>
<div> </div>
<div><br> </div>
<div><span class="gmail_quote">2008/3/28, Robert Hart <<a href="mailto:rhart@fsu.edu">rhart@fsu.edu</a>>:</span>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">Eduardo,<br><br>The problem is possibly that you are defining PP, u, and v only at the lowest level, since<br>
the initial dimension environment opens the files with z=1, and you immediately define<br>the new variables in that environment.<br><br>Change your script to this, and it should work:<br><br>'sdfopen d:\ncep-ncar\uwnd.mon.mean.nc'<br>
'sdfopen d:\air.mon.mean.nc'<br><br>"define uisen="isen(uwnd.1,air.2,lev,320)<br>"d uisen"<br><br>Bob Hart<br>FSU Meteorology<br><br><br>----- Original Message -----<br>From: Eduardo Agosta Scarel <<a href="mailto:eduardo.agosta@GMAIL.COM">eduardo.agosta@GMAIL.COM</a>><br>
Date: Friday, March 28, 2008 5:56 pm<br>Subject: use of the script <a href="http://isen.gs">isen.gs</a><br>To: <a href="mailto:GRADSUSR@LIST.CINECA.IT">GRADSUSR@LIST.CINECA.IT</a><br><br>> Hello everybody!<br>><br>> I am trying to use the script <a href="http://isen.gs">isen.gs</a> to interpolate a variable field<br>
> on isentropic surface but the final result is that there are no<br>> griddedvalues to interpolate. Perhaps the problem is the use of<br>> 'lev' instead of<br>> pressure values? How to solve it?<br>><br>
> What did I do wrong?<br>><br>> Please any suggestion to use the script will be much appreciated.<br>><br>> Eduardo<br>><br>> Below follows the main script used to test it.<br>><br>><br>><br>
> 'reinit'<br>><br>> 'sdfopen d:\ncep-ncar\uwnd.mon.mean.nc'<br>> 'sdfopen d:\air.mon.mean.nc'<br>><br>><br>> 'define PP ='lev<br>> 'define u= uwnd.1'<br>> 'define v=vwnd.2'<br>
> 'define t= air.3'<br>><br>> "set lev 1050 150"<br>><br>> "define uisen="isen(u,t,PP,320)<br>><br>> "set z 1"<br>> "d uisen"<br>><br>><br>><br>
> function isen(field,tgrid,pgrid,tlev)<br>> *-------------------------------------------------------------------<br>> ---<br>> * Bob Hart (<a href="mailto:hart@ems.psu.edu">hart@ems.psu.edu</a>) / PSU Meteorology<br>
> * 2/26/1999<br>> *<br>> * 2/26/99 - Fixed a bug that caused the script to crash on<br>> * certain machines.<br>> *<br>> * GrADS function to interpolate within a 3-D grid to a specified<br>> * isentropic level. Can also be used on non-pressure level data, such<br>
> * as sigma or eta-coordinate output where pressure is a function<br>> * of time and grid level. Can be used to create isentropic PV<br>> surfaces* (examples are given at end of documentation just prior to<br>> * function.)<br>
> *<br>> * Advantages: Easy to use, no UDFs. Disadvantages: Can take 5-20<br>> secs.*<br>> * Arguments:<br>> * field = name of 3-D grid to interpolate<br>> *<br>> * tgrid = name of 3-D grid holding temperature values (deg K) at<br>
> each* gridpoint.<br>> *<br>> * pgrid = name of 3-D grid holding pressure values (mb) at each<br>> gridpoint* If you are using regular pressure-level data,<br>> this should be<br>> * set to the builtin GrADS variable 'lev'.<br>
> *<br>> * tlev = theta-level (deg K) at which to interpolate<br>> *<br>> * Function returns: defined grid interp holding interpolated values<br>> *<br>> * NOTE: YOU NEED TO INCLUDE A COPY OF THIS FUNCTION IN YOUR SCRIPT<br>
> *<br>> * NOTE: Areas having tlev below bottom level or above upper level<br>> * will be undefined in output field. Extrapolation is NOT<br>> * performed!!<br>> *<br>> *-------------------------------------------------------------------<br>
> -----<br>> *<br>> * EXAMPLE FUNCTION CALLS:<br>> *<br>> * Sample variables: u = u-wind in m/s<br>> * v = v-wind in m/s<br>> * w = vertical velocity<br>> * t = temperature in K<br>
> * PP = pressure data in mb<br>> *<br>> * 1) Display vertical velocity field on 320K surface:<br>> *<br>> * "d "isen(w,t,PP,320)<br>> *<br>> * 2) Create & Display colorized streamlines on 320K surface:<br>
> *<br>> * "define u320="isen(u,t,PP,320)<br>> * "define v320="isen(v,t,PP,320)<br>> * "set z 1"<br>> * "set gxout stream"<br>> * "d u320;v320;mag(u320,v320)"<br>
> *<br>> * 3) Create & display a 320K isentropic PV surface:<br>> *<br>> * "set lev 1050 150"<br>> * "define coriol=2*7.29e-5*sin(lat*3.1415/180)"<br>> * "define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"<br>
> * "define<br>> dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"*<br>> "define dt=t(z-1)*pow(1000/PP(z-1),0.286)-<br>> t(z+1)*pow(1000/PP(z+1),0.286)"<br>> * "define dp=100*(PP(z-1)-PP(z+1))"<br>
> * "define dtdp=dt/dp"<br>> * "define part1="isen(dvdx,t,PP,320)<br>> * "define part2="isen(dudy,t,PP,320)<br>> * "define part3="isen(dtdp,t,PP,320)<br>> * "define pv320=-9.8*(coriol+part1-part2)*part3"<br>
> * "set z 1"<br>> * "d pv320"<br>> *<br>> * PROBLEMS: Send email to Bob Hart (<a href="mailto:hart@ems.psu.edu">hart@ems.psu.edu</a>)<br>> *<br>> *-------------------------------------------------------------------<br>
> ----<br>> *-------------------- BEGINNING OF FUNCTION ------------------------<br>> ----<br>> *-------------------------------------------------------------------<br>> ----<br>><br>> * Get initial dimensions of dataset so that exit dimensions will be<br>
> * same<br>><br>> "q dims"<br>> rec=sublin(result,4)<br>> ztype=subwrd(rec,3)<br>> if (ztype = "fixed")<br>> zmin=subwrd(rec,9)<br>> zmax=zmin<br>> else<br>> zmin=subwrd(rec,11)<br>
> zmax=subwrd(rec,13)<br>> endif<br>><br>> * Get full z-dimensions of dataset.<br>><br>> "q file"<br>> rec=sublin(result,5)<br>> zsize=subwrd(rec,9)<br>><br>> * Determine spatially varying bounding pressure levels for isen<br>
> surface* tabove = theta-value at level above ; tbelow = theta value<br>> at level<br>> * below for each gridpoint<br>><br>> "set z 1 "zsize<br>> "define theta="tgrid"*pow(1000/"pgrid",0.286)"<br>
> "set z 2 "zsize<br>> "define thetam="tgrid"(z-1)*pow(1000/"pgrid"(z-1),0.286)"<br>> "set z 1 "zsize-1<br>> "define thetap="tgrid"(z+1)*pow(1000/"pgrid"(z+1),0.286)"<br>
><br>> "define tabove=0.5*maskout<br>> (theta,theta-"tlev")+0.5*maskout(theta,"tlev"-thetam)"<br>> "define tbelow=0.5*maskout<br>> (theta,thetap-"tlev")+0.5*maskout(theta,"tlev"-theta)"<br>
><br>> * Isolate field values at bounding pressure levels<br>> * fabove = requested field value above isen surface<br>> * fbelow = requested field value below isen surface<br>><br>> "define fabove=tabove*0+"field<br>
> "define fbelow=tbelow*0+"field<br>><br>> "set z 1"<br>><br>> * Turn this 3-D grid of values (mostly undefined) into a 2-D isen<br>> layer<br>> * If more than one layer is valid (rare), take the mean of all the<br>
> * valid levels. Not the best way to deal with the multi-layer issue,<br>> * but works well, rarely if ever impacts output, and is quick.<br>> * Ideally, only the upper most level would be used. However, this<br>
> * is not easily done using current GrADS intrinsic functions.<br>><br>> "define fabove=mean(fabove,z=1,z="zsize")"<br>> "define fbelow=mean(fbelow,z=1,z="zsize")"<br>> "define tabove=mean(tabove,z=1,z="zsize")"<br>
> "define tbelow=mean(tbelow,z=1,z="zsize")"<br>><br>> * Finally, interpolate linearly in theta and create isen surface.<br>> * Linear interpolation in theta works b/c it scales as height,<br>
> * or log-P, from Poisson equation for pot temp.<br>><br>> "set z "zmin " " zmax<br>><br>> "define slope=(fabove-fbelow)/(tabove-tbelow)"<br>> "define b=fbelow-slope*tbelow"<br>
> "define interp=slope*"tlev"+b"<br>><br>> * variable interp now holds isentropic field and its named it returned<br>> * for use by the user.<br>><br>> say "Done. Newly defined variable interp has "tlev"K "field"-field."<br>
><br>> return(interp)<br>><br></blockquote></div><br>