**
** Syntax: arbitrary_cross_section.gs [input_filename] [outputpath/filename] [minlat] [minlon] [maxlat] [maxlon] [begin_lev] [end_lev] [parameter] [units] [contours/shaded]
**
** ex: run arbitrary_cross_section.gs /raid2/fnmoc_models/fnmoc_model_058078240/g2006/g2006120500_058078240.000.ctl cross_section.gif 29 -81 31 -74 300 0 wtmpdsl NM sh
**
** Grads-script that takes the input_file/lat/lon/lev and draws
** a cross section from point 1 to point 2 with a varying level.
** The latitude 1 and longitude 1 must always be less than the
** latitude 2 and longitude 2.  Depending on which model is being
** used the orientataion of the level varies. The argument
** contour/shading refers to whether you want countour lines or
** a shaded output.  The default is contours, but if the user
** enters sh for argument 10 then a you will get a shaded output.
** The x-axis is always the distance from point 1 to point 2 in
** the specified units and the y-axis is the vertical. Units
** outputed can be Kilometer which can be passed in as KM, feet which can be
** passed in as EN, meters which can be passed in as M, or nautical
** miles which can be passed in as NM.
**
**  Author: Heiner Körnich
**  Created: unknown
**
**  Author: David Ramsaur
**  Update: 6 December 2006
**      -Description, Arguments, and Input Path
**
function main(args)

_inputpath=subwrd(args,1)

'open '_inputpath

* 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(args,7)
_lev2=subwrd(args,8)
if (_lev1>_lev2)
opti=',-u'
else
opti=''
endif

* set points
_outputfilename=subwrd(args,2)
_lat1=subwrd(args,3)
_lon1=subwrd(args,4)
_lat2=subwrd(args,5)
_lon2=subwrd(args,6)
_variable=subwrd(args,9)
_units=subwrd(args,10)
_type=subwrd(args,11)



* 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)

if (_units=EN)
    distance = distance*3280.8399' feet'
endif
if (_units=NM)
    distance = distance*0.5399' nautical miles'
endif
if (_units=M)
    distance = distance*1000' meters'
endif
if (_units=KM)
    distance = distance' kilometers'
endif

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

'printim '_outputfilename


return

'reninit'
