use of the script isen.gs

[rhart at FSU.EDU]

Eduardo

The version I sent you worked for me

What version of grads and can you give me a url to your data

Bob
Sent via BlackBerry from T-Mobile

-----Original Message-----
From: Eduardo Agosta Scarel <eduardo.agosta at GMAIL.COM>

Date: Sat, 29 Mar 2008 11:36:38
To:GRADSUSR at LIST.CINECA.IT
Subject: Re: use of the script isen.gs


Dear Robert,   I thank you for your suggestion.   I did it the way shown below, but it did not work either! Any other suggestion?   Thanks in advance. Eduardo. __   'reinit' 'sdfopen d:\ncep-ncar\uwnd.mon.mean.nc''sdfopen d:\ncep-ncar\air.mon.mean.nc' 'set lon -180 0' 'set lat 0 90' 'set lev 1050 10' 'set t 1' 'define uisen='isen(uwnd.1,air.2,lev,320) "set z 1"'set t 1' "d uisen" ___     Â
  2008/3/28, Robert Hart <rhart at fsu.edu <mailto:rhart at fsu.edu> >: Eduardo,
The problem is possibly that you are defining  PP, u, and v only at the lowest level, since the initial dimension environment opens the files with z=1, and you immediately definethe new variables in that environment.
Change your script to this, and it should work:
'sdfopen d:\ncep-ncar\uwnd.mon.mean.nc' 'sdfopen d:\air.mon.mean.nc'
"define uisen="isen(uwnd.1,air.2,lev,320)"d uisen"
Bob HartFSU Meteorology

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