diagonal cross section

[Frank Colby]

I had a lot of trouble making arbitrary cross-sections, but have had
success using a script that George Bryan at PSU put together.  Here's
one example that you should be able to modify.  It's quite easy to add
variables by making additional collections, and even to make calculations.

Frank

*------------------------------------------------------------------
*  GrADS script:  crs.gs
*
*  Script purpose:  (see below)
*
*  Author:  unknown.  Provided by Bob Hart or Sytske Kimball (I can't
*                     remember which).
*
*  Last modified:  10 July 2000  by George H. Bryan (PSU)
*
*********************************************************************
* The following lines will display an arbitrary X section
* from one specified point to another.
*
* lon1 is the westernmost longitude point
* lon2 is the easternmost longitude point
* lat1 is the latitude that corresponds to lon1
* lat2 is the latitude that corresponds to lon2
*
* The loop is used to interpolate between points in
* the arbitrary cross section.  This code will plot
* any cross section as long as you specify the points.
* My code plots cross sections of PV after I calculated
* PV on 11 pressure surfaces.  I have another script
* that plots cross sections of potential temperature, and
* the code is very similar to this, except theta is substituted
* for PV.
*
* Many thanks to Brian Doty at COLA for his help with this code.
*
*******************************************************************
****CHANGED THIS TO ENHANCE FLEXIBILITY
***you have to preset the dimensions, t,z,lat,lon
say 'variable'
pull var
say 'using which contour interval'
pull civ
*extract the data to fix the x-section
'set lat 40 43'
'set lon -70 -66'
'q dim'
rec1=sublin(result,2)
lon11=subwrd(rec1,6)
lon12=subwrd(rec1,8)
x11=subwrd(rec1,11)
x12=subwrd(rec1,13)
rec2=sublin(result,3)
lat11=subwrd(rec2,6)
lat12=subwrd(rec2,8)
y11=subwrd(rec2,11)
y12=subwrd(rec2,13)
dum=sublin(result,5)
time=subwrd(dum,6)
*say time
t1=1
t2=25
while(t1<=t2)

lon1=lon11
lon2=lon12
lat1=lat11
lat2=lat12
ztop=15
zbot=1

say lon1' 'lon2' 'lat1' 'lat2
say x11' 'x12' 'y11' 'y12

********
*'clear'
'set grads off'
'set zlog on'
'set z 'zbot' 'ztop
'set x 'x11
'set y 'y11
*'set y 'y11' 'y12
lon = lon1
lat = lat1
'set t 't1

'collect 1 free'

  while (lon<=lon2)
     lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)
     'collect 1 gr2stn('var','lon','lat')'
     lon = lon + 0.1
*    say lat
*    say lon
  endwhile
* IMPORTANT TO SET X CORRECTLY
  'set x 'x11' 'x12
  lon1=(lon1-1)*2
  lon2=(lon2-1)*2
  xdis=(lon2-lon1)*(lon2-lon1)
  ydis=(lat2-lat1)*(lat2-lat1)
  distance=sqrt( xdis + ydis)
  'set cint 'civ
  'set clab on'
  'set clopts -1 -1 0.15'
*  'set ccolor rainbow'
  'set gxout contour'
  'd coll2gr(1,-u)'
'draw title X-section from 'lat11'N to 'lat12'N, 'time' 't1
'set strsiz 0.08'
'draw xlab longitude'
'draw ylab sigma levels'
*endif

***reset the chosen x-section lat/lon for next plot
'set lat 'lat11' 'lat12
'set lon 'lon11' 'lon12

say t1
say '<enter>=next, q=quit, b=back, t<enter>n<enter>=n-th step'
pull ans
if (ans ="q"); break; endif;
if (ans = "b")
t1=t1-2
'q dim'
rec=sublin(result,5)
time=subwrd(rec,6)
VT=substr(time,1,3)
'draw title  blabalblalbalbalbla'var1' 'VT
endif
if (ans = "t"); say 'Time step =';pull num;
if (num > t2); say 'Timestep too high!'; num=t2-1; endif;
t1=num-1
endif
t1 = t1 + 1
*'c'
endwhile
************************************************************************
function sqrt(i)

*-----------------------------------------
* return cosine of i, where i is in radians
*------------------------------------------

"set gxout stat"
"d sqrt("i")"
rec=sublin(result,8)
val=subwrd(rec,4)
return(val)

************************************************************************


-------------- next part --------------
A non-text attachment was scrubbed...
Name: Frank_Colby.vcf
Type: text/x-vcard
Size: 237 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20051219/52fb8668/attachment.vcf