* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* bathymap.gs
*
* This script overlays a land mask (bathymetry) in plots showing vertical cross
* section of various fields in the ocean. The bathymetry data required for this 
* purpose is obtained by the program from a bathymetry data (in NetCDF format)
* specified by the user in the input argument. 
*
* Known bugs/difficulties with the code. 
*  1) Upon calling this function, it have to close all input files opened in 
*     GrADS. This step is required since the space time dimensions of bathymetry
*     file used by this function may not agree with that of other files that
*     are opened in GrADS when this function is called, leading to errors. 
*  2) Information on  
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
* Author : Martin V. Mathew, vmartinmathew@gmail.com
* Created on 17 October 2012 
* Modified on 23 December 2012                           
*LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
function main(args)

* There are defaults for land mask color
* User must specify the input bathymetry file and its display variable name.

  if (args='')
    say 'Usage: bathymap.gs ifile dispvar <flip levels flag> <shading color> '
    say '        ifile  ==> Filename of bathymetry data in NetCDF format'
    say '        dispvar==> Variable name of bathyetry data in the given i/p file '
    say '        <flip levels flag> optional==> By default bathymetry data is '
    say '        plotted assuming negative y-axis value for depths below sea level' 
    say '        and positive values for those above sea level. But some datasets' 
    say '        have positive depth level values for depths below sea level. in '
    say '        such cases assigning "-1" to this argument will help bathymap.gs ' 
    say '        to plot the bathymetry data properly. By default this flag takes ' 
    say '        a value "1". '
    say '        <shading color>  optional==> By default a gray shading is '
    say '        given to the land mask. This optional argument allows one to '
    say '        set any other color of is preference. Use one of the 16 GrADS '
    say '        default colors.' 
    return
  else
    bathyfile  = subwrd(args,1)
    bathyvar   = subwrd(args,2)
    invrtlvl   = subwrd(args,3)
    shdeclr    = subwrd(args,4)
    if (invrtlvl = '')   ; invrtlvl = 1  ; endif
    if (shdeclr = '')   ; shdeclr = 15  ; endif
  endif

*----------------------------------------------
* Obtaining current dimension environment
*----------------------------------------------
'q dim'

lonline= sublin(result,2)
latline= sublin(result,3)
levline= sublin(result,4)
if (subwrd(result,8)=subwrd(result,21))

say 'Can not plot bathymetry data because of error in data dimension.'
say 'Please make sure that the dimension enironment is not modified after '
say 'plotting the required vertical cross section of the field.' 
say 'Only one of the x or y dimension is allowed to vary. Current dimension'
say ' environment is:'
say lonline 
say latline
return

else;* if (subwrd(result,8)=subwrd(result,21))
if (subwrd(result,8)='varying')
lonMnf=subwrd(lonline,6)
lonMxf=subwrd(lonline,8)
latMnf=subwrd(latline,6)
latMxf=latMnf
say 'Plotting bathymetry between longitudes 'lonMnf' to 'lonMxf
 
else;*if (subwrd(result,8)='varying')
lonMnf=subwrd(lonline,6)
lonMxf=lonMnf
latMnf=subwrd(latline,6)
latMxf=subwrd(latline,8)
say 'Plotting bathymetry between latitudes 'latMnf' to 'latMxf 
endif ;* if (subwrd(result,8)='varying')

endif;* if (subwrd(result,8)=subwrd(result,21))

if (subwrd(levline,3)='varying')
levMnf=subwrd(levline,6)
levMxf=subwrd(levline,8)

if levMnf>=0
say 'Vertical levels in the input data are all positive values. If the bathymetry '  
say 'data is plotted wrongly, then consider using optional -1 flag while calling '
say 'bathymap.gs to correct it.'  
endif;* if levMnf>=0

endif;*if (subwrd(levline,3)='varying')

*------------------------------------------------------------------------------
* Closing all currently opened files. This step is required since the space time
* dimensions of bathymetry file used by this function may not agree with that of
* other files that are currently open, leading to errors. 
*------------------------------------------------------------------------------
* Obtaining number of file currently opened in GraADS 
tstflg=1
incr=1
while tstflg=1
'q file 'incr
if (subwrd(result,3) ='Error:')
tstflg=0
else ;* if (subwrd(result,3) ='Error:')
incr=incr+1
endif;* if (subwrd(result,3) ='Error:')
endwhile ;* while tstflg=1
flno=incr-1

*Closing the files one by one
while flno>0
close' 'flno
flno=flno-1
endwhile ;*while flno>0 

*----------------------------------------------
* Opening the input bathymetry file
*----------------------------------------------
'sdfopen 'bathyfile

*----------------------------------------------
*  Obtaining file number of the bathymetry file
*----------------------------------------------
anslin5=sublin(result,3)
flno =subwrd(anslin5,8)

*---------------------------------------------------------------------------
* Setting latitude and longitude coordinates
*---------------------------------------------------------------------------
'set lat 'latMnf' 'latMxf
'set lon 'lonMnf' 'lonMxf
'set z 1'
'set t 1'

*---------------------------------------------------------------------------
*  Generating array containing bathymetry data along required cross section
*  and corresponding geographical coordinates.
*---------------------------------------------------------------------------
'set gxout print'
'd 'bathyvar'.'flno
maxpts=subwrd(result,4)
strtVls=10;* data for first point is the 10th word of the result
cnt=1

while(cnt<=maxpts)
bdtaary.2.cnt=subwrd(result,(strtVls+cnt-1))*invrtlvl
cnt=cnt+1
endwhile;*while(cnt<maxpts)

'set gxout stat'
'd 'bathyvar'.'flno
strtVls=16;* data for first point is the 16th word of the result
cnt=1

while(cnt<=maxpts)
bdtaary.1.cnt=subwrd(result,(strtVls+cnt-1))
cnt=cnt+1
endwhile;*while(cnt<maxpts)


*---------------------------------------------------------------------------
* Drawing the bathymetry land mask data
*---------------------------------------------------------------------------
'set gxout shaded'
'q gxinfo'
btmy=subwrd(result,20)
MaxTpY=subwrd(result,22)

*say 'Qrying location 'bdtaary.1.1' 'bdtaary.2.1
'q w2xy 'bdtaary.1.1' 'bdtaary.2.1
prevX = subwrd(result,3)
prevY = subwrd(result,6)

*truncating y values of bathymetry data that exceed the height of figure window
if prevY>MaxTpY; prevY=MaxTpY; endif

cnt=2 ;* suspect

'set line 'shdeclr
while(cnt<=maxpts)

if (cnt!=2)
prevX =crntX 
prevY =crntY
endif ;* if (cnt!=2)

*say 'Qrying location 'bdtaary.1.cnt' 'bdtaary.2.cnt
'q w2xy 'bdtaary.1.cnt' 'bdtaary.2.cnt
crntX = subwrd(result,3)
crntY = subwrd(result,6)

*truncating y values of bathymetry data that exceed the height of figure window
if crntY>MaxTpY; crntY=MaxTpY; endif

if ((prevY>btmy)&(crntY>btmy))
'draw polyf 'prevX' 'prevY' 'crntX' 'crntY' 'crntX' 'btmy' 'prevX' 'btmy
*say 'polyf x1 y1 x2 y2 x2 y3 x1 y3'prevX' 'prevY' 'crntX' 'crntY' 'crntX' 'btmy' 'prevX' 'btmy
 
endif ;*if ((prevY>btmy)&(crntY>btmy))

cnt=cnt+1
endwhile;*while(cnt<maxpts)

close' 'flno;* Closing the bathymetry file
return