************************************************************************
*
*  To use:  "grads -l"
*
*  Script plots SkewT/logP sounding(s) using Bob Hart's "plotskew.gs"
*  function
*
************************************************************************
function skewt_lsf

'reinit'

*-----------------------------------------------------------------------
*  Batch mode?
*-----------------------------------------------------------------------

batch_mode = 'no'

*-----------------------------------------------------------------------
*  Pause between plots?
*-----------------------------------------------------------------------

do_pause = 'yes'

*-----------------------------------------------------------------------
*  GrADS metafile output?
*-----------------------------------------------------------------------

do_print = 'yes'

if( do_print = 'yes' )
  'enable print skewt_lsf.gx'
endif

*-----------------------------------------------------------------------
*  physical constants and conversion factors
*-----------------------------------------------------------------------

tfreeze = 273.15
ms_to_knots = 1.0/0.51

*-----------------------------------------------------------------------
*  specify GrADS ctl file for the sounding
*-----------------------------------------------------------------------

ctl_file = '/net/cjs/gcss/case2/old_IFA/skewt/sound-gcss2-ifa.ctl'

*-----------------------------------------------------------------------
*  specify (x,y) location (grid point numbers), the bottom and top "zdef"
*  levels for the sounding (to define the pressure levels), and time level
*  number(s) for the sounding...
*-----------------------------------------------------------------------

xloc = 1
yloc = 1
zbot = 1 ; ztop = 38
tlev = '1 29'

*-----------------------------------------------------------------------
*  open GrADS file...
*-----------------------------------------------------------------------

'open 'ctl_file

*-----------------------------------------------------------------------
*  define and plot sounding(s)...
*-----------------------------------------------------------------------

n=1 ; tt = subwrd(tlev,n)
while ( tt != '' )
*
*  reset GrADS for this sounding...
*
  'reset'
*
*  allow external GrADS functions
*
  rc = gsfallow("on")
*
*  define the latitude, longitude and date for the sounding...
*
  say
  say 'Time level = 'tt
  say
  'set x 'xloc
  'set y 'yloc
  'set z 'zbot' 'ztop
  'set t 'tt
  'q dims'
  rec2 = sublin(result,2) ; xl = subwrd(rec2,6)
  rec3 = sublin(result,3) ; yl = subwrd(rec3,6)
  rec5 = sublin(result,5) ; tl = subwrd(rec5,6)
*
*  define sounding for Bob Hart's "plotskew.gs" function...
*
  'tfr = 'tfreeze
  'tc = t-tfr'
  'tdc = td-tfr'
  'ws = wspd*'ms_to_knots
  'wd = wdir'
*
*  plot sounding in pressure coordinates...
*
  'set grads off'
  'set gxout line'
  'c'
  'd tc'
  'd tdc'
  'draw title \\\ Sounding temperature and dewpoint'
  'draw xlab (C) \\\'
  'draw ylab \\\ Pressure (hPa)'
  if( do_pause = 'yes' )
    pause()
  endif
  'c'
  'd ws'
  'draw title \\\ Sounding wind speed'
  'draw xlab (knots) \\\'
  'draw ylab \\\ Pressure (hPa)'
  if( do_pause = 'yes' )
    pause()
  endif
  'c'
  'd wd'
  'draw title \\\ Sounding wind direction'
  'draw xlab (degrees from 0 North) \\\'
  'draw ylab \\\ Pressure (hPa)'
  if( do_pause = 'yes' )
    pause()
  endif
*
*  plot sounding using Bob Hart's "plotskew.gs" function...
*
  'set t 1'
  'set x 1'
  'set y 1'
  'set lev 1050 100'
  rc = plotskew(tc,tdc,ws,wd)
  if( rc != 0 )
    say
    say '"plotskew" return code = 'rc
    say 'exiting script...'
    say
    exit
  endif
  'draw title \\\ 'tl', (lon,lat) = ('xl','yl')'

  'undefine tc' ; 'undefine tdc' ; 'undefine ws' ; 'undefine wd'

  if( do_print = 'yes' )
    'print'
  endif

  if( do_pause = 'yes' )
    pause()
  endif

  n=n+1 ; tt = subwrd(tlev,n)
endwhile

if( do_print = 'yes' )
  'disable print'
endif

say
say '************************************************'
say ' Finished with this script.'
say '************************************************'
say

if( batch_mode = 'yes' )
  quit
endif

return
