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

'reinit'

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

batch_mode = 'no'

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

do_pause = 'yes'

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

do_print = 'yes'

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

*-----------------------------------------------------------------------
*  Physical constants and conversion factors
*-----------------------------------------------------------------------

tfreeze = 273.16
ms_to_knots = 1.0/0.51

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

ctl_file = '/archive/cjs/zetac/old/KWAJEX/old_history/SKEWT.06:30z24Jul1999_2.ctl'

*-----------------------------------------------------------------------
*  Specify names for pressure, temperature, dewpoint, wind speed, and
*  wind direction (include any "region" suffix here)
*-----------------------------------------------------------------------

region_title = 'Domain Average'
pres_var = pp
temp_var = ta
dewp_var = td
um_var = um
vm_var = vm

*region_title = 'Deep Convective Region Average'
*pres_var = pp_dcreg
*temp_var = ta_dcreg
*dewp_var = td_dcreg
*um_var = um_dcreg
*vm_var = vm_dcreg

*-----------------------------------------------------------------------
*  Specify the pressure level for the top of the sounding (hPa)
*  and time level number(s) for the sounding...
*-----------------------------------------------------------------------

ptop = 100
tlev = '1 96 192'

*-----------------------------------------------------------------------
*  Allow external GrADS functions
*-----------------------------------------------------------------------

rc = gsfallow("on")

*-----------------------------------------------------------------------
*  Open GrADS ctl file...
*-----------------------------------------------------------------------

'open 'ctl_file

*-----------------------------------------------------------------------
*  Define GrADS dataset file information...
*-----------------------------------------------------------------------

'q file'
rec5 = sublin(result,5)
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt = subwrd(rec5,12)

'q ctlinfo'
say result
dset_rec = sublin(result,1) ; ieee_file = subwrd(dset_rec,2)
title_rec = sublin(result,2)
undef_rec = sublin(result,3) ; undef_value = subwrd(undef_rec,2)
*
*  search for the "tdef" record to define the GrADS "dt"...
*
nrec = 4
rec = sublin(result,nrec) ; rec_word1 = subwrd(rec,1)
while ( rec_word1 != tdef & rec != '' )
  nrec = nrec+1
  rec = sublin(result,nrec) ; rec_word1 = subwrd(rec,1)
endwhile
if( rec_word1 != tdef )
  say
  say 'did not find a record word = tdef...'
  say '...nrec 'nrec
  say 'Record: 'rec
  say 'Record word 1: 'rec_word1
  say '...exiting script'
  say
  exit
else
  say
  tdef_rec = rec
  say '"tdef" record (nrec='nrec') is:'
  say tdef_rec
  dt = subwrd(tdef_rec,5)
  say 'GrADS "dt" is: 'dt
  say
endif
*
*  search for the "pres_var" record to define the physical units...
*  (start at the "vars" record following the "tdef" record...)
*
nrec = nrec+1
rec = sublin(result,nrec); nvars = subwrd(rec,2)
n=1
while ( rec_word1 != pres_var & n <= nvars )
  nrec = nrec+1
  rec = sublin(result,nrec) ; rec_word1 = subwrd(rec,1)
  n=n+1
endwhile
if( rec_word1 != pres_var )
  say
  say 'did not find a record word = 'pres_var'...'
  say '...nrec 'nrec
  say 'Record: 'rec
  say 'Record word 1: 'rec_word1
  say '...exiting script'
  say
  exit
else
  say
  pres_var_rec = rec
  say '"'pres_var'" record (nrec='nrec') is:'
  say pres_var_rec
  pres_var_units = subwrd(pres_var_rec,5)
  say '"'pres_var'" units: 'pres_var_units
*
*  set "pscale" to convert pressure values to (hPa)...
*
  if( pres_var_units = '(pa)' | pres_var_units = '(Pa)' )
    pscale = 1./100.
  else
    pscale = 1.
  endif
  say '"'pres_var'" scale: 'pscale
  say
endif

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

n=1 ; tt = subwrd(tlev,n)
while ( tt != '' )
*
*  reset GrADS for this sounding...
*
  'reset'
*
*  define the latitude, longitude and date for the sounding...
*
  say
  say 'Time level = 'tt
  say
  'set dfile 1'
  xmid = (nx+1)/2
  ymid = (ny+1)/2
  'set x 'xmid
  'set y 'ymid
  '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 the pressure levels for the sounding...
*
  'set x 1'
  'set y 1'
  'set z 1 'nz
  'set t 'tt
  'pb = aave('pres_var',x=1,x='nx',y=1,y='ny')'
  'set gxout print'
  zz=1 ; 'set z 'zz ; 'd pb' ; rec = sublin(result,2) ; pval = subwrd(rec,1)
  if( pval = undef_value )
    say
    say 'Caution:  pressure value is undefined...'
    say '   first level pval = 'pval
    say 'continuing with next sounding...'
    say
    n=n+1 ; tt = subwrd(tlev,n)
    continue
  else
*
*  scale to (hPa)...
*
    pval = pval*pscale
  endif
  np=0 ; pvals = ''
  while ( pval >= ptop )
    np=np+1
    pvals = ''pvals' 'pval
    zz=zz+1
    if( zz <= nz )
      'set z 'zz
      'd pb' ; rec = sublin(result,2) ; pval = subwrd(rec,1)
      if( pval = undef_value )
        say
        say 'Caution:  pressure value is undefined...'
        say '   level "'zz'" pval = 'pval
        say 'continuing with next sounding...'
        say
        do_next_sounding = 'yes'
        break
      else
*
*  scale to (hPa)...
*
        pval = pval*pscale
        do_next_sounding = 'no'
      endif
    else
      say
      say 'Caution:  lowest pressure in input data >= ptop = 'ptop
      say '   level "'zz'" pval = 'pval
      say 'continuing...'
      say
      break
    endif
  endwhile
  if( do_next_sounding = 'yes' )
    say '"do_next_sounding" = 'do_next_sounding
    n=n+1 ; tt = subwrd(tlev,n)
    continue
  endif
  if( np = 0 )
    if( pval = undef_value )
      say
      say 'Caution:  pressure value is undefined...'
      say '   first level pval = 'pval
      say 'continuing...'
      say
      continue
    else
      say
      say '...no pressure levels written'
      say '   first level pval = 'pval
      say 'exiting script...'
      say
      exit
    endif
  else
    say
    say 'Sounding will contain 'np' pressure levels...'
    say '(number of input data levels = 'nz')'
    say '...minimum pressure for the sounding = 'subwrd(pvals,np)
    say
    nlevs = np
  endif
*
*  Define sounding for GrADS fwrite file...
*
  'set x 1'
  'set y 1'
  'set z 1 'nlevs
  'set t 'tt
  'tcb = aave('temp_var',x=1,x='nx',y=1,y='ny')-'tfreeze
  'tdcb = aave('dewp_var',x=1,x='nx',y=1,y='ny')-'tfreeze
  'umb = aave('um_var',x=1,x='nx',y=1,y='ny')'
  'vmb = aave('vm_var',x=1,x='nx',y=1,y='ny')'
  'wskb = sqrt(umb*umb+vmb*vmb)*'ms_to_knots
*
* Code below is adapted from an email to GrADS Listserv on 1/2/2002 by
* Brian Doty
*
  'db1 = -atan2(umb,-vmb)*180/3.141592654'
  'db2 = const(maskout(db1,db1),0,-u)'
  'db3 = const(360+maskout(db1,-db2),0,-u)'
  'wdb1 = db2+db3'
*
* Code below is adapted from an email to GrADS Listserv on 9/20/2004 by:
*
* Ken Takahashi G.       e-mail: ken@atmos.washington.edu
* (206) 685-9305
* Dept of Atmospheric Sciences
* University of Washington, Box 351640
*
  'wdb = atan2(-umb,-vmb)*180/3.141592654'
  'wdb2 = const(wdb,1)*(const(maskout(wdb,wdb),0,-u)+const(maskout(360+wdb,-wdb),0,-u))'
*
*  Plot sounding in height coordinates...
*
  'c'
  'set grads off'
  'set gxout line'
  'd tcb'
  'd tdcb'
  'draw title \\\ Sounding temperature and dewpoint'
  'draw xlab (C) \\\'
  'draw ylab \\\ Height (km)'
  if( do_pause = 'yes' )
    pause()
  endif
  'c'
  'd wskb'
  'draw title \\\ Sounding wind speed'
  'draw xlab (knots) \\\'
  'draw ylab \\\ Height (km)'
  if( do_pause = 'yes' )
    pause()
  endif
  'c'
  'd wdb1'
  'd wdb2'
  'draw title \\\ Sounding wind direction'
  'draw xlab (degrees from 0 North) \\\'
  'draw ylab \\\ Height (km)'
  l=1
  while ( l <= nlevs )
    'set z 'l
    'd wdb1-wdb2'
    rec = sublin(result,1) ; diff = subwrd(rec,4)
    if( diff != 0 )
      say '...non-zero difference = 'diff' at l = 'l
      exit
    endif
    l=l+1
  endwhile
  if( do_pause = 'yes' )
    pause()
  endif
*
*  Write results to GrADS fwrite file...
*
  '!rm skewt.ieee'
  'set gxout fwrite'
  'set fwrite -le ./skewt.ieee'
  'set z 1 'nlevs
  'wdb = wdb1'
  'd tcb'
  'd tdcb'
  'd wskb'
  'd wdb'
  'disable fwrite'
  'reset'
*
*  Begin writing a local GrADS ctl file
*
  '!rm skewt.ctl'

  res = write ('skewt.ctl','dset ^skewt.ieee')
  res = write ('skewt.ctl','title Domain Averaged Sounding')
  res = write ('skewt.ctl','options little_endian')
  res = write ('skewt.ctl',''undef_rec'')
  res = write ('skewt.ctl','xdef 1 levels 'xl'')
  res = write ('skewt.ctl','ydef 1 levels 'yl'')
  res = write ('skewt.ctl','zdef 'nlevs' levels')
  l=1
  while ( l <= nlevs )
    pval = subwrd(pvals,l)
    res = write ('skewt.ctl',' 'pval'')
    l=l+1
  endwhile
  res = write ('skewt.ctl','tdef 1 linear 'tl' 'dt'')
  res = write ('skewt.ctl','vars 4')
  res = write ('skewt.ctl',' tempb 'nlevs' 99 temperature (C)')
  res = write ('skewt.ctl',' dewpb 'nlevs' 99 dewpoint temperature (C)')
  res = write ('skewt.ctl',' wspdb 'nlevs' 99 wind speed (knots)')
  res = write ('skewt.ctl',' wdirb 'nlevs' 99 wind direction (degrees from 0 North)')
  res = write ('skewt.ctl','endvars')
  res = close ('skewt.ctl')

  status = sublin(res,1)
  if( status != 0 )
    say
    say '"skewt.ctl" status number = 'status
    say 'exiting script...'
    say
    exit
  endif
*
*  Plot sounding in pressure coordinates...
*
  'open skewt.ctl'
  'set dfile 2'
  'set x 1'
  'set y 1'
  'set z 1 'nlevs
  'set t 1'
  'q file' ; say result
  'q dims' ; say result
  'set grads off'
  'set gxout line'
  'c'
  'd tempb'
  'd dewpb'
  'draw title \\\ Sounding temperature and dewpoint'
  'draw xlab (C) \\\'
  'draw ylab \\\ Pressure (hPa)'
  if( do_pause = 'yes' )
    pause()
  endif
  'c'
  'd wspdb'
  'draw title \\\ Sounding wind speed'
  'draw xlab (knots) \\\'
  'draw ylab \\\ Pressure (hPa)'
  if( do_pause = 'yes' )
    pause()
  endif
  'c'
  'd wdirb'
  '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 'ptop
  rc = plotskew(tempb,dewpb,wspdb,wdirb)
  if( rc != 0 )
    say
    say '"plotskew" return code = 'rc
    say 'exiting script...'
    say
    exit
  endif
  'draw title \\\ 'region_title' at 'tl

  'undefine pb'
  'undefine tcb' ; 'undefine tdcb' ; 'undefine wskb' ; 'undefine wdb'
  'close 2'

  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
