* Grads-script for a cross section
* usage: 1) first open file
*        2) set desired vertical span
*        3) run: section(lon1,lat1,lon2,lat2,variable,type)
*                with lon1<lon2
* (type='' for contours, type='sh' for shadings)
* problem: cross section is a straight line on a lat-lon rectangular plot
*          not the shortest distance
function main(args)

* resolution in lon-direction
dlon=1

* count words in args
nargs=1
while(subwrd(args,nargs)!='')
nargs=nargs+1
endwhile
nargs=nargs-1
* check consistency of script call
if (nargs<5)
say 'section.gs plots a cross-section from the point (lon1,lat1)'
say ' to the point (lon2,lat2) for the variable var.'
say 'Usage of section.gs:'
say ' 1) open file'
say ' 2) set desired vertical span, eg. set lev 1000 100'
say ' 3) run: section lon1 lat1 lon2 lat2 var [sh]'
say '    (lon1<lon2, optional argument sh yields a shaded plot)'
return
endif

* check vertical coordinate
'q dims'
res=sublin(result,4)
lev1=subwrd(res,6)
lev2=subwrd(res,8)
if (lev1>lev2)
opti=',-u'
else
opti=''
endif

* set points
lon1=subwrd(args,1)
lat1=subwrd(args,2)
lon2=subwrd(args,3)
lat2=subwrd(args,4)
variable=subwrd(args,5)
type=subwrd(args,6)

* calculate distance between points
'set x 1'
'set y 1'
'set z 1'
'pi=3.1415927'
'p1= cos(('lon1-lon2')/180*pi)'
'p2 = cos(('lat1-lat2')/180*pi)'
'p3 = cos(('lat1+lat2')/180*pi)'

'dp= (acos(((p1*(p2+p3))+(p2-p3))/2))'
'd dp*6.371e3'
distance=subwrd(result,4)
say 'distance= 'distance
'set lev 'lev1' 'lev2


* interpolate on cross-section
'set x 1'
'set y 1'
lon = lon1
'collect 1 free'
while (lon <= lon2)
  lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)
  'collect 1 gr2stn('variable','lon','lat')'
  lon = lon + dlon
endwhile

* plot interpolated data
'set x 14 16'
*'set xaxis 'lon1' 'lon2
'set xaxis 0 'distance
'set clab on'
if(type='sh')
'set gxout shaded'
'd coll2gr(1'opti')'
*'d coll2gr(1)'
'cbar'
else
'set gxout contour'
'set clab forced'
*'set cint .5'
*'d coll2gr(1,-u)'
'd coll2gr(1'opti')'
endif


return
