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 drad=500 "open "ctl hgtvar="hgtprs" slpvar="prmslmsl/100" prsvar="lev" "set display color white" "clear" "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) levs = "600 550 500 450 400 350 300" "set time "starttime " " endtime "q dims" rec=sublin(result,5) startt=subwrd(rec,11) endt=subwrd(rec,13) tt=startt ii=0 "set mpdset hires" "set gxout fwrite" "set fwrite vtu.grads" 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) if (xlon < 0) xlon=xlon+360 endif xlat=subwrd(posns,2) pmin=subwrd(posns,3) if (xlon > -900 & xlat > -900 & xlon != "" ) "set gxout stat" Sumx=0 Sumy=0 Sumx2=0 Sumxy=0 zz=1 "set gxout contour" "set lat "xlat-30" " xlat+30 "set lon "xlon-30 " "xlon+30 "define dx=(lon-"xlon")*111*cos(lat*3.1415/180)" "define dy=(lat-"xlat")*111" "define dist=pow(dx*dx+dy*dy,0.5)" while (subwrd(levs,zz) != "") plev=subwrd(levs,zz) "set z 1" "define hlev="pinterp(hgtvar,prsvar,plev) "define grid=maskout(hlev,"drad"-dist)" "clear" "set gxout grfill" "set cint 10" "set grads off" "d grid" "draw title "name"("year") "time"\"plev "hPa" if (tt > startt) "q w2xy "xlon " " xlat xloc1=subwrd(result,3) yloc1=subwrd(result,6) "q w2xy "xlonold " " xlatold xloc2=subwrd(result,3) yloc2=subwrd(result,6) "draw line "xloc1 " " yloc1 " " xloc2 " " yloc2 endif "set gxout stat" "d grid" rec=sublin(result,8) minval=subwrd(rec,4) maxval=subwrd(rec,5) grad.zz=maxval-minval plev.zz=math_log(plev) Sumx=Sumx+plev.zz Sumy=Sumy+grad.zz Sumx2=Sumx2+(plev.zz)*(plev.zz) Sumxy=Sumxy+(plev.zz)*(grad.zz) say plev " " grad.zz "draw xlab DZ="substr(maxval-minval,1,6) zz=zz+1 endwhile n=zz-1 maxp=subwrd(levs,1) minp=subwrd(levs,n) slope=(Sumx*Sumy-n*Sumxy)/(Sumx*Sumx-n*Sumx2) if (slope > 2000 | slope < -2000) slope=9.999e20 endif say time " " xlon " " xlat " " slope " " stype ii=ii+1 type.ii=stype "set gxout fwrite" if (ii = 1) startctl = time endif "set x 1" "set y 1" "set z 1" "set t 1" "d "slope endif "clear" xlatold=xlat xlonold=xlon tt=tt+1 endwhile "disable fwrite" "!echo 'dset ^vtu.grads' > vtu.ctl" "!echo 'undef 9.999E+20' >> vtu.ctl" "!echo 'options "bytes"' >> vtu.ctl" "!echo 'title EP Prof File' >> vtu.ctl" "!echo 'ydef 1 linear -90.000000 2.5' >> vtu.ctl" "!echo 'xdef 1 linear 0.000000 2.500000' >> vtu.ctl" "!echo 'tdef "ii" linear "startctl" 6hr' >> vtu.ctl" "!echo 'zdef 1 levels 1000' >> vtu.ctl" "!echo 'vars 1' >> vtu.ctl" "!echo 'vtu 1 99 slope of lin reg' >> vtu.ctl" "!echo 'ENDVARS' >> vtu.ctl" "quit" 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) ************************************************************************* 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 1 "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)