function pvinterp(field,tgrid,ugrid,vgrid,pgrid,pvlev,smooth) *---------------------------------------------------------------------- * Clark Evans (evans36-at-uwm-dot-edu) 7/28/12 * Logic following related code by Bob Hart (rhart@fsu.edu) * * GrADS function to interpolate within a 3-D grid to a specified * potential vorticity (PV) level. * * Arguments: * field = name of 3-D grid to interpolate * * tgrid = name of 3-D grid holding temperature values (deg K) at each * gridpoint. * * ugrid = name of 3-D grid holding u-wind components (m s-1) at each * gridpoint. * * vgrid = name of 3-D grid holding v-wind components (m s-1) 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 built-in GrADS variable "lev". * * pvlev = PV-level (PVU; 1e-6 m K Pa-1 s-3) at which to interpolate * * smooth = number of times to run each variable through a nine-point * smoother prior to interpolation. Since the salient features * that this (and other similar) routine is designed to highlight * are predominantly larger-scale features, a small number (0-5) * is recommended for synoptic to meso-alpha scale data. * Larger numbers (10+) are recommended for finer scale data. * Set to zero for no smoothing. * * Function returns: defined grid interp holding interpolated values * * EXAMPLE -- defining potential temperature on the 1.5 PVU surface: * (temp=temperature, U=u-wind, V=v-wind data, zmax=highest data level) * "set z 1 zmax" * "define temppv="pvinterp(temp,temp,U,V,lev,1.5) * "define levpv="pvinterp(lev,temp,U,V,lev,1.5) * "define thetapv=temppv*pow(1000/levpv,0.286)" * "set z 1" * "d thetapv" * * NOTE: You must include a copy of this function in any script that calls it. * * NOTE: Areas having pvlev below the bottom or above the upper level * of available data will be undefined in the output field. No * extrapolation is performed on any data. *----------------------------------------------------------------------- *-------------------- BEGINNING OF FUNCTION ---------------------------- *----------------------------------------------------------------------- * Run everything through a smoother. i=1 while(i<=smooth) "define "field"=smth9("field")" "define "ugrid"=smth9("ugrid")" "define "vgrid"=smth9("vgrid")" "define "tgrid"=smth9("tgrid")" i=i+1 endwhile * Get initial dimensions of dataset so that exit dimensions will be * same. Thus, zmin is the lowest z level specified by the user, while * zmax is the highest z level specified by the user. "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. Thus, zsize is the total number of * z levels in the dataset. "q file" rec=sublin(result,5) zsize=subwrd(rec,9) * Determine spatially varying bounding pressure levels for surface * pvabove = PV value at level above ; pvbelow = PV value at level * below for each gridpoint "set z 1 "zsize "define dudy=cdiff("ugrid",y)/(111177*cdiff(lat,y))" "define dvdx=cdiff("vgrid",x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))" "define dt=("tgrid"(z+1)*pow(1000/"pgrid"(z+1),0.286))-("tgrid"(z-1)*pow(1000/"pgrid"(z-1),0.286))" "define dtdp=dt/(100*("pgrid"(z+1)-"pgrid"(z-1)))" * Define PV across the domain at the appropriate level and at levels from 2-end * and 1-endminus1. Then, isolate the appropriate levels for each file. "set z 1 "zsize "define pv=-9.8*(2*7.29e-5*sin(lat*3.1415/180)+((cdiff("vgrid",x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180)))-(cdiff("ugrid",y)/(111177*cdiff(lat,y)))))*((("tgrid"(z+1)*pow(1000/"pgrid"(z+1),0.286))-("tgrid"(z-1)*pow(1000/"pgrid"(z-1),0.286)))/(100*("pgrid"(z+1)-"pgrid"(z-1))))*1e6" "set z 2 "zsize "define pvm=-9.8*(2*7.29e-5*sin(lat*3.1415/180)+((cdiff("vgrid"(z-1),x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180)))-(cdiff("ugrid"(z-1),y)/(111177*cdiff(lat,y)))))*((("tgrid"*pow(1000/"pgrid",0.286))-("tgrid"(z-2)*pow(1000/"pgrid"(z-2),0.286)))/(100*("pgrid"-"pgrid"(z-2))))*1e6" "set z 1 "zsize-1 "define pvp=-9.8*(2*7.29e-5*sin(lat*3.1415/180)+((cdiff("vgrid"(z+1),x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180)))-(cdiff("ugrid"(z+1),y)/(111177*cdiff(lat,y)))))*((("tgrid"(z+2)*pow(1000/"pgrid"(z+2),0.286))-("tgrid"*pow(1000/"pgrid",0.286)))/(100*("pgrid"(z+2)-"pgrid")))*1e6" "define pvabove=0.5*maskout(pv,pv-"pvlev")+0.5*maskout(pv,"pvlev"-pvm)" "define pvbelow=0.5*maskout(pv,pvp-"pvlev")+0.5*maskout(pv,"pvlev"-pv)" * Isolate field values at bounding PV levels, * fabove = requested field value above PV surface * fbelow = requested field value below PV surface * pabove = pressure grid values above PV surface * pbelow = pressure grid values below PV surface "define fabove=pvabove*0+"field "define fbelow=pvbelow*0+"field "define pabove=pvabove*0+"pgrid "define pbelow=pvbelow*0+"pgrid "set z 1" * Turn this 3-D grid of values (mostly undefined) into a 2-D PV layer. * First, find the number of valid levels containing data. * Next, from this, find which of the levels is closest to the ground. * Then, find the mean of the entire field (code below) to convert * from a 3-D grid to a 2-D grid. Since all other layers are undefined, * this simply reduces the dimensionality from 3 to 2. * NOTE: As of 7/28/12, all but the 3D->2D is not yet implemented. "define fabove=mean(fabove,z=1,z="zsize")" "define fbelow=mean(fbelow,z=1,z="zsize")" "define pvabove=mean(pvabove,z=1,z="zsize")" "define pvbelow=mean(pvbelow,z=1,z="zsize")" * Finally, interpolate linearly in theta and create PV surface. * Linear (y=mx+b) interpolation for PV doesn't work quite as well as * theta, but given dependence on theta and pressure, it should work okay. "set z "zmin " " zmax "define slope=(fabove-fbelow)/(pvabove-pvbelow)" "define b=fbelow-slope*pvbelow" "define interp=slope*"pvlev"+b" * Variable interp now holds PV field and is returned * for use by the user. say "Done. Newly defined variable interp has "pvlev" PVU "field" field." return(interp)