in=read(desc.dat) val=sublin(in,2) starttime=getgdtg(val) in=read(desc.dat) val=sublin(in,2) endtime=getgdtg(val) in=read(desc.dat) name=sublin(in,2) in=read(desc.dat) ctl=sublin(in,2) trackfile=stormtrack "open "ctl "q config" rec=sublin(result,1) bytes=subwrd(rec,4) if (bytes = "big-endian") bytes="big_endian" else bytes="little_endian" endif year=substr(starttime,9,4) "set display color white" "clear" "set time "starttime " " endtime "q dims" rec=sublin(result,5) startt=subwrd(rec,11) endt=subwrd(rec,13) tt=startt "set gxout fwrite" "set fwrite size.grads" ii=1 while (tt <= endt) "set t "tt "q dims" rec=sublin(result,5) time=subwrd(rec,6) posns=getpos(year,name,time,trackfile) xlon=subwrd(posns,1) xlat=subwrd(posns,2) wrad=17 inc=0.25 "set lat "xlat-10" " xlat+10 "set lon "xlon-16 " "xlon+16 "define wind=mag(ugrdprs(lev=925),vgrdprs(lev=925))" "define dy=111*(lat-"xlat")" "define dx=cos(lat*3.1415/180)*111*(lon-"xlon")" "define dist=pow(dx*dx+dy*dy,0.5)" "define ysize=111*"inc"+dist*0" "define xsize=cos(lat*3.1415/180)*111*"inc"+dist*0" "define area=maskout(xsize*ysize,wind-"wrad")" "set gxout stat" "d area" rec=sublin(result,10) tot=subwrd(rec,2) size=math_sqrt(tot/3.1415) * Write out information to GrADS file "set gxout fwrite" "set x 1" "set y 1" "set z 1" "set t 1" "d "size tt=tt+1 ii=ii+1 xlatold=xlat xlonold=xlon endwhile "disable fwrite" * Create .ctl file ii=ii-1 "!echo 'dset ^size.grads' > size.ctl" "!echo 'undef 9.999E+20' >> size.ctl" "!echo 'title Front Prof File' >> size.ctl" "!echo 'options "bytes"' >> size.ctl" "!echo 'ydef 1 linear -90.000000 2.5' >> size.ctl" "!echo 'xdef 1 linear 0.000000 2.500000' >> size.ctl" "!echo 'tdef "ii" linear "starttime" 6hr' >> size.ctl" "!echo 'zdef 1 linear 1000 5' >> size.ctl" "!echo 'vars 1' >> size.ctl" "!echo 'size 1 99 dir' >> size.ctl" "!echo 'ENDVARS' >> size.ctl" "quit" function pinterp(field,pgrid,plev) *---------------------------------------------------------------------- * Bob Hart (hart@ems.psu.edu) / PSU Meteorology * 2/26/1999 * * 2/26/1999 - Fixed a bug that caused the script to crash on * certain machines. * * GrADS function to interpolate within a 3-D grid to a specified * pressure 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. * * Advantages: Easy to use, no UDFs. Disadvantages: Can take a few secs. * * Arguments: * field = name of 3-D grid to interpolate * * pgrid = name of 3-D grid holding pressure values at each gridpoint * If you are using regular pressure-level data, this should be * set to the builtin GrADS variable 'lev'. * * plev = pressure level 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 plev 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 * t = temperature in K * PP = pressure data in mb * * 1) Display a temperature field interpolated to 925mb * * "d "pinterp(t,PP,925) * * 2) Display a streamline field interpolated to 550mb * * "define u550="pinterp(u,PP,550) * "define v550="pinterp(v,PP,550) * "set gxout stream" * "d u550;v550" * * * PROBLEMS: Send email to Bob Hart (hart@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 * epsilon here will take care of rare instance where user requests * an interpolated level that already exists in the file. Necessary * to make maskout functions below work; otherwise, the * routine finds no set of levels that is valid for interpolation * Essentially, adding epsilon below 'tips' the * level selection to one side when two choices are available without * altering the output of the interpolated results. epsilon=0.0001 * Get full vertical dimensions of dataset "q file" rec=sublin(result,5) zsize=subwrd(rec,9) "set z 1 "zsize "q dims" rec=sublin(result,4) pmax=subwrd(rec,6) pmin=subwrd(rec,8) * Shift epsilon tip to other side in case user requests lowest level If (plev = pmax) epsilon=-epsilon Endif * Determine spatially varying bounding pressure levels for p-surface * pabove = pressure-value at level above ; pbelow = pressure value at level * below for each gridpoint "set z 1 "zsize-1 "define pabove=0.5*maskout("pgrid","plev+epsilon"-"pgrid")+0.5*maskout("pgrid","pgrid"(z-1)-"plev+epsilon")" "set z 2 "zsize "define pbelow=0.5*maskout("pgrid","plev+epsilon"-"pgrid"(z+1))+0.5*maskout("pgrid","pgrid"-"plev+epsilon")" * Isolate field values at bounding pressure levels * fabove = requested field value above pressure surface * fbelow = requested field value below pressure surface "define fabove=pabove*0+"field "define fbelow=pbelow*0+"field * Turn this 3-D grid of values (mostly undefined) into a 2-D press layer. * mean is used here only for simplicity. "set z 1" "define pabove=mean(pabove,z=1,z="zsize")" "define fabove=mean(fabove,z=1,z="zsize")" "define pbelow=mean(pbelow,z=1,z="zsize")" "define fbelow=mean(fbelow,z=1,z="zsize")" * Finally, interpolate linearly in log-pressure and create surface. "set z "zmin " " zmax "define interp=fbelow+log(pbelow/"plev")*(fabove-fbelow)/log(pbelow/pabove)" *say "Done. Newly defined variable interp has "plev"mb "field"-field." return(interp) function getpos(year,name,reqtime,trackfile) storm=year"/"name rc=0 lon=-9999 lat=-9999 tstamp="" dummy=read(trackfile) while (rc = 0 & tstamp!=reqtime) out=read(trackfile) rc=sublin(out,1) data=sublin(out,2) lat=subwrd(data,2) lon=subwrd(data,3) date=subwrd(data,1) tstamp=getgdtg(date) endwhile rc=close(trackfile) return(lon " " lat) function getmonth(month) months = "JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC" mo2=subwrd(months,month) return(mo2) function getgdtg(dtg) year=substr(dtg,1,4) month=substr(dtg,5,2) day=substr(dtg,7,2) hr=substr(dtg,9,2) month2=getmonth(month) string=hr"Z"day""month2""year return(string)