*------------------------------------------------------------------
*  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 42.9 41.7'
'set lon -72.0 -70.0'
'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=36
while(t1<=t2)

lon1=lon11
lon2=lon12
lat1=lat11
lat2=lat12
ztop=20
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.2
*    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 ccolor 1'
  'set gxout contour'
  'd coll2gr(1,-u)'
'q dim'
rec=sublin(result,5)
time=subwrd(rec,6)
VT=substr(time,1,3)
'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)
 
************************************************************************

