* demo.gs v1.1                a sample grads script
*
* generally grads needs 3 environment variables to be set up
*   (1) path, (2) gaddir, (3) gascrp, and location of control files.
*
* However for demo purposes, assume the 3 variables are not set.
* Works by cd-ing to directory with (1) grads, (2) grad scripts
*   (3) grads data files and running grads from that directory,
*   and (4) the location of the control files.
*
*
disk='/cdrom'
disk='e:'
day='14mar93'
mon='mar'
year='93'

count=1

say ' '
say ' '
say ' '
say ' '
say ' '
say ' '
say ' '
say ' '
say '              NCEP/NCAR Reanalysis Demo for BAMS CD-ROM v1.1'
say ' '
say 'This CD-ROM contains a preliminary climatology (1982-1994), selected'
say '  monthly means (1982-1994) and selected daily means (1993) from'
say '  the NCEP/NCAR Reanalysis.  The data can be displayed on a PC'
say '  using PC-GrADS (included) and on many workstations (workstation'
say '  versions of GrADS are available on the internet).  The data are'
say '  stored in GRIB, a WMO standard, and can be easily read using a'
say '  good GRIB decoder (wgrib is included).  Before using the data for'
say '  serious work, read the documentation as some of the data were '
say '  processed differently than in the more complete Reanalysis data sets.'
say ' '
say 'Acknowledgments:'
say '  Brian Doty (COLA) for writing GrADS (c) 1988-1995 Brian Doty.'
say '  This demo was written in the GrADS scripting language and only'
say "  shows a fraction of GrADS' capabilities."
say ' '
say '  The NCEP/NCAR Reanalysis Project has been supported by the NOAA'
say '  Office of Global Programs, the National Weather Service and the'
say '  National Science Foundation.'
say ' '
say ' '
prompt '                                                     Enter return: '
pull c

while(1)
   'reinit'
   say ' '
   say ' '
   say ' '
   say ' '
   say ' '
   say ' '
   say ' '
   say ' '
   say ' '
   'clear'
   say '                   NMC/NCAR Reanalysis Demo '
   say ' '
   say '1)   mean sea-level pressure (mslp, mb) ' day ' daily ave'
   say '2)   500 mb omega (Pa/s), 850 mb T (C) ' day
   say '3)   2m T (C), 10m wind (m/s)' day
   say '4)   850 mb geo. height (m) ' day '  (departure from ' mon ' normal)'
   say '5)   200 mb geo. height (m), mslp (mb) ' day ' +- 4 days'
   say '6)   model precip. and obs. outgoing long wave radiation (OLR) ' mon year
   say '7)   model OLR - obs OLR (W/m**2) ' mon year
   say '8)   sea-surface temperature (C) ' mon year
   say '9)   zonal wind (m/s) ' mon year ' (slow)'
   say '10)  200 mb zonal wind (m/s) on equator 19' year
   say '11)  850 mb temperature (C) at 40N 80W'
   say '12)  500 mb height (m) ' day ' to +10 days'
   say '13)  potential vorticity (m**2/s/kg) on 330K surface ' day ' -4..3 days'
   say '14)  200 mb U anomaly (165W-110W,5S-5N) 5-mo running mean (slow)'
   say 'q)   Quit'
   say ' '
   say 'd)   interactive display/plot daily fields'
   say 'm)   interactive display/plot monthly fields'
*  say 't)   change time'
*  say 'g)   end up in grads'
   say ' '
   say ' '
   prompt 'Enter Choice: '
   pull c

   if (c = '')
       c = count 
       count = count + 1
       if (count = 13)
          count=1
       endif
   endif

   if (c = '1')
      'open ' disk '/daily/mslp.ctl'
      'set lat 20 70'
      'set lon 230 330'
      'set mpdset mres'
      'set time ' day
      'set gxout shaded'
      'clear'
      'set map 0'
      'd mslp*0.01'
      'set gxout contour'
      'set ccolor 0'
      'd mslp*0.01'
      'draw title mean sea level pressure (mb) ' day
      prompt 'press return for menu '
      pull x
   endif
   if (c = '2')
      'open ' disk '/daily/t850.ctl'
      'open ' disk '/daily/w500.ctl'
      'set lat 20 70'
      'set lon 230 330'
      'set mpdset mres'
      'set time ' day
      'set gxout shaded'
      'clear'
      'set map 0'
      'd w.2(lev=500)'
      'run cbar.gs'
      'set gxout contour'
      'set ccolor 0'
      'd t-273.16'
      'draw title 500 mb omega (Pa/s), 850 mb T (C) ' day
      prompt 'press return for menu '
      pull x
   endif
   if (c = '3')
      'open ' disk '/daily/t2m.ctl'
      'open ' disk '/daily/u10m.ctl'
      'open ' disk '/daily/v10m.ctl'
      'set lat 20 70'
      'set lon 230 330'
      'set mpdset mres'
      'set time ' day
      'set gxout shaded'
      'clear'
      'set map 0'
      'd t2m-273.16'
      'run cbar.gs'
      'set ccolor 0'
      'd u10m.2;v10m.3'
      'draw title 2m T (C), 10m wind (m/s) ' day "+6 hours"
      prompt 'press return for menu '
      pull x
   endif     
   if (c = '4')
      'open ' disk '/daily/z850.ctl'
      'open ' disk '/climate/prflcl.ctl'
      'set lev 850'
      'set time ' day
      'set lat 20 70'
      'set lon 230 330'
      'set mpdset mres'
      'clear'
      'set map 0'
      'set gxout shaded'
      'd z.1 - z.2(time=' mon '82)'
      'run cbar.gs'
      'set ccolor 0'
      'set gxout contour'
      'set cint 60'
      'd z.1'
      'draw title ' day ' 850 mb height (m) (shaded=dep from normal)'
      prompt 'press return for menu '
      pull x
   endif
   if (c = '5')
      'open ' disk '/daily/mslp.ctl'
      'open ' disk '/daily/z200.ctl'
      'set lat 20 70'
      'set lon 230 330'
      'clear'
      'set map 0'
      'set time ' day
      'q dims'
      line1=sublin(result,5)
      itime=subwrd(line1,9)
      itime1=itime+4
      itime=itime-4
*unix 'set dbuff on'
      'set lev 200'
      while (itime <= itime1)
          'set t ' itime
*dos      clear
          clear
          'set gxout shaded'
          'd mslp*0.01'
          'run cbar.gs'
          'set gxout contour'
          'd z.2'
          'q time'
          line1=sublin(result,1)
          time1=subwrd(line1,3)
          time=substr(time1,4,9)
          'draw title ' time ' 200 mb Z, mean sea-level pressure'
*unix     'swap'
          itime = itime + 1
          prompt 'press return for next'
          pull x
      endwhile
*unix 'set dbuff off'
   endif
   if (c = '6')
      'open ' disk '/monthly/prfl.ctl'
      'open ' disk '/observed/obsmnolr.ctl'
      'set lev 500'
      'set time 01' mon year
      'set gxout shaded'
      'clear'
      'd prate*86400'
      'run cbar.gs'
      'set gxout contour'
      'set ccolor 1'
      'd olr.2'
      'draw title ' mon year ' model precip (mm/day), obs OLR (W/m**2)'
      prompt 'press return for menu '
      pull x
   endif
   if (c = '7')
      'open ' disk '/monthly/prfl.ctl'
      'open ' disk '/observed/obsmnolr.ctl'
      'clear'
      'set map 0'
      'set gxout shaded'
      'd ulwrftoa-olr.2'
      'run cbar.gs'
      'draw title ' mon year ' model OLR - obs OLR (W/m**2)'
      prompt 'press return for menu '
      pull x
   endif
   if (c = '8')
      'open ' disk '/monthly/prfl.ctl'
      'open ' disk '/fixed/mask.ctl'
      'open ' disk '/fixed/orog.ctl'

      'set gxout shaded'
      'set time 01' mon year
      'clear'
      'set cint 5'
      'd maskout(tsfc-273.16, 0.5 - mask.2(t=1))'
      'set gxout contour'
      'set cint 5'
      'd maskout(tsfc-273.16, 0.5 - mask.2(t=1))'
      'd maskout(orog.3(t=1)*0.001,mask.2(t=1)-0.5)'

      'draw title observed SST (C) and orographic height (km) ' mon year
      prompt 'press return for menu '
      pull x
   endif
   if (c = '9')
      'open ' disk '/monthly/prfl.ctl'
      'set lev 1000 10'
      'define uu = u'
      'clear'
      'set lon 0'
      'define uz=ave(uu,lon=0,lon=360)'
      'clear'
      'set gxout shaded'
      'set cint 5'
      'd uz'     
      'set gxout contour'
      'set ccolor 0'
      'set cint 5'
      'd uz'     
      'draw title zonal wind (m/s) ' mon year
      'draw ylab Pressure (mb)'
      prompt 'press return for menu '
      pull x
   endif
   if (c = '10')
      'open ' disk '/monthly/prfl.ctl'
      'set time jan' year ' dec' year
      'set lat 0'
      'set lev 200'
      'clear'
      'set gxout shaded'
      'd u'
      'run cbar.gs'
      'set gxout contour'
      'd u'
      'draw title 200 mb zonal wind (m/s) 19' year
      prompt 'press return for menu '
      pull x
   endif
   if (c = '11') 
      'open ' disk '/monthly/prfl.ctl'
      'set time jan82 dec94'
      'set lat 40'
      'set lon -80'
      'set lev 850'
      'clear'
      'd t - 273.16'
      'draw title 850mb tmp (C) 40N 80W'
      'draw ylab temp (C)'
      'draw xlab Date'
      prompt 'press return for menu '
      pull x
   endif
   if (c = '12') 
      'open ' disk '/daily/z500.ctl'
      'set time ' day
      'q dims'
      line1=sublin(result,5)
      itime=subwrd(line1,9)
      itime1=itime+10
    
      'set t ' itime ' ' itime1
      'set lat 20 90'
      'set lon -270 90'
      'set lat 20 60'
      'set lon -150 -60'
      'set mproj nps'
      'clear'
      'draw title 500 mb geopotential height'
      'd z'
      prompt 'press return for menu '
      pull x
   endif
   if (c = '13') 
      'open ' disk '/daily/ipv330.ctl'
      'set time ' day
      'q dims'
      line1=sublin(result,5)
      itime=subwrd(line1,9)
      itime1=itime+3
      itime=itime-4
      'set lat 20 90'
      'set lon -270 90'
      'set mproj nps'
      'set gxout shaded'
      while (itime <= itime1)
          'set t ' itime
          'clear'
          'set map 0'
          'd ipv'
          'run cbar.gs'
          'q time'
          line1=sublin(result,1)
          time1=subwrd(line1,3)
          time=substr(time1,4,9)
          'draw title pot. vorticity on 330K sfc ' time
          itime=itime+1
          prompt 'press return for next plot '
          pull x
      endwhile
   endif
   if (c = '14')
      'open ' disk '/monthly/prfl.ctl'
      'open ' disk '/climate/prflcl.ctl'
      'set time jan82 dec82'
      'set lev 200'
      'set lat 0'
      'set lon 0'
      'define uc=aave(u.2,lon=-165,lon=-110,lat=-5,lat=5)'
      'set time jan82 dec94'
      'modify uc seasonal'
      'define uu=aave(u,lon=-165,lon=-110,lat=-5,lat=5)'
      'clear'
      'define us=ave(uu-uc,t-2,t+2)'
      'd us'
      'draw title 200mb U anom (165W-110W,5S-5N) 5-mo run. mean'
      prompt 'press return for menu '
      pull x
   endif
   if (c = 't' | c = 'T')
      say 'Enter new date: '
      say '  daily data only available for 1993'
      say '  monthly data only available from 1982-1994'
      say ' '
      say '  If you want to display daily fields, the year must be 93'
      say '  the date must be in the form of ddmmmyy'
      say '  where dd is the day of the month (2 digits)'
      say '  mmm is first 3 letters of the month; example jan, feb, mar ..'
      say '  yy is year (last 2 digits)'
      prompt 'enter date: (ex. 04mar93) '
      pull day
      mon=substr(day,3,3)
      year=substr(day,6,2)
      say ' '
      say 'you entered date=' day ' month=' mon ' year=19' year
      say 'If this is not correct use menu option t again'
      say ' '
   endif
   if (c = 'd' | c = 'D')
      'run daily.gs'
   endif
   if (c = 'm' | c = 'M')
      'run monthly.gs'
   endif
   if (c = 'g' | c = 'G')
      exit
   endif
   if (c = 'q' | c = 'Q')
      break
   endif
endwhile
'quit'