************************************************************************
*
*  To use:  "grads -p"
*
*  Script calculates and plots zonally averaged cross-sections of
*  height interpolated "radon222" from the GrADS files specified below
*  (with an option to overlay the level of minimum temperature).
*
*  The sample standard deviation is also calculated and plotted
*  (contour plot overlay).
*
*  From statistics (reference: "Statistical Concepts and Methods"
*  by Bhattacharyya and Johnson (1977)), the sample variance s**2
*  for a set of "N" observations of a variable "xi" (x1,...,xN)
*  is defined as:
*
*    s**2 = [(x1-xbar)**2 + (x2-xbar)**2 + (x3-xbar)**2 +....
*          + (xN-xbar)**2 ]/(N-1) ,  (from page 34)
*
*  where the sample mean = xbar = (x1+x2+x3+...+xN)/N , (from page 27)
*
*  Note, the sample variance sums over "N" observations, and divides the
*  sum by "N-1"...
*
*  The sample standard deviation is defined as the square root of the
*  sample variance (from page 35):
*
*    s = sqrt(s**2)
*
************************************************************************
function radon (args)

'reinit'

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

rc = gsfallow("on")

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

do_pause = 'no'

************************************************************************
*  GrADS meta file output?
************************************************************************

do_print = 'yes'

************************************************************************
*  Information for GrADS files
************************************************************************

dir.1 = '/net2/cjs/ljd/lima/ch3i1/history/atmos_month'
dir.2 = '/net2/cjs/lwh/fms/lima/m45_am2p14_ch3i/history/atmos_month'
sdir_months = 'Jul'
mod_years = '1983 1998'
ctl_file.1 = 'radon222.o_zl.ctl'
ctl_file.2 = 'Radon222.o_zl.ctl'
ctl_file_temp = 'temp.o_zl.ctl'

************************************************************************
*  Information for the plots
************************************************************************

plot_min_max = 'no'
plot_sd_min_max = 'no'

var_name = 'radon222'
*
*  Input Units:  10e-21 vmr
*
*  Note:  10e-21 vmr = 1.53 pCi/SCM,
*              1 pCi = 37.8 mBq
*
var_scale = 1
var_units = '(10`a-21`n volume mixing ratio)'
*
*  Conversion factor from 10e-21 vmr to mBq/SCM:
*
*   ([37.8 mBq]/[1 pCi])*([1.53 pCi/SCM]/[10e-21 vmr])
*         = 57.834 mBq/SCM /(10e-21 vmr)
*
*var_scale = 57.834
*var_units = '(mBq/SCM)'
*
*  Height levels (scale GrADS ctl file levels from "m" to "km")
*
m_to_km = 1./1000.
zlevs = '0 18'
ylabel_int = 2
ylabel = 'Height (km)'
*
*  Define the contour levels and contour colors for the variable...
*
define_clevcols('vmr')
*
*  Define the contour levels for the sample standard deviation...
*
sd_level = '0.2 0.5 0.8 1.1'
*sd_style = ' 5   3   2   1'
sd_style = ' 1   1   1   1'
*sd_label = 'off'
*sd_label = 'on'
sd_label = 'manual' ; csiz = 0.08 ; cformat = '%3.1f'
*
*  Plot level of minimum temperature?
*
plot_min_temp = 'yes'
var_name_temp = 'temp'

************************************************************************
*  Plot variables' page limits
************************************************************************

pvar_vpage.1 = '0 8.5 5.0 9.5'
pvar_vpage.2 = '0 8.5 0.5 5.0'
pvar_vpage.1 = '0 8.5 4.8 9.3'
pvar_vpage.2 = '0 8.5 0.3 4.8'
pvar_parea   = '1 7.0 0.5 4.0'

xl_pvar_parea = subwrd(pvar_parea,1)
yt_pvar_parea = subwrd(pvar_parea,4)

************************************************************************
*  Title information for the plots
************************************************************************

mod_title = 'FMS lima 1983-1998'
var_title = 'Radon-222 'var_units
avg_title = 'July'

figure_panel.1 = '(a)' ; model.1 = 'AM2-D'
figure_panel.2 = '(b)' ; model.2 = 'AM2'

************************************************************************
*  GrADS metafile for output
************************************************************************

if( do_print = 'yes' )
  'enable print fig1ab_sd-manual.gx'
endif

*-----------------------------------------------------------------------
*  Verify height levels...
*-----------------------------------------------------------------------

l=1
zlv.l = subwrd(zlevs,l)
if( zlv.l = '' )
  say
  say 'No height levels specified in variable "zlevs"...'
  say 'Exiting script'
  say
  exit
else
  while ( zlv.l != '' )
    l=l+1
    zlv.l = subwrd(zlevs,l)
  endwhile
  nl = l-1
  if( nl > 2 )
    say
    say 'More than two height levels ("nl" = 'nl') specified in variable "zlevs"...'
    say 'Exiting script'
    say
    exit
  else
    if( zlv.2 <= zlv.1 )
      say
      say 'Height level "2" (= 'zlv.2') not greater than height level "1" (= 'zlv.1')...'
      say 'Exiting script'
      say
      exit
    endif
  endif
endif

*-----------------------------------------------------------------------
say
say 'Calculate and plot zonally averaged cross-sections for July average...'
say
*-----------------------------------------------------------------------
*
*  local year variables...
*
y1 = subwrd(mod_years,1)
y2 = subwrd(mod_years,2)
nyears = y2-y1+1
*
*  local month variable...
*
m1 = subwrd(sdir_months,1)
*
*  Start a new frame, and plot a title at top of page...
*
ptitle_top(var_title,avg_title,pvar_parea)
*
*  Loop over models...
*
m=1
while ( m <= 2 )

  say
  say '------------'
  say 'Model: 'model.m
  say '------------'
  say

  if( plot_min_temp = 'yes' )
    say
    say 'Find and save level of minimum temperature for the zonally averaged'
    say 'July average temperature...'
    say
    say '...start temperature time average sum at year = 'y1
    say
*
*  Make local ctl file with scaled vertical levels...
*
    ctl_file = ''dir.m'/'m1'/'y1'/'ctl_file_temp
    ctl_file_scale = ''var_name_temp'.'m1'.'y1'.scale.ctl'
    scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

    say 'open GrADS file for 'm1', 'y1
    say '...'ctl_file_scale
    'open 'ctl_file_scale
    say 'GrADS file ready'

    'set dfile 1'
    'query file'
    rec5 = sublin(result,5)
    nxt = subwrd(rec5,3) ; nyt = subwrd(rec5,6) ; nzt = subwrd(rec5,9)

    'set x 1 'nxt
    'set y 1 'nyt
    'set z 1 'nzt
    'set t 1'
    'define tsum = 'var_name_temp

    'close 1'
    '!rm 'ctl_file_scale

    say
    say '...continue with running sum of July average temperature for remaining years'
    say '   (number of years = 'nyears')'
    say

    y=y1+1
    while ( y <= y2 )

      say
      say '*** 'y' ***'
      say

      ctl_file = ''dir.m'/'m1'/'y'/'ctl_file_temp
      ctl_file_scale = ''var_name_temp'.'m1'.'y'.scale.ctl'
      scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

      say 'open GrADS file for 'm1', 'y
      say '...'ctl_file_scale
      'open 'ctl_file_scale
      say 'GrADS file ready'

      say
      say '...increment temperature time average sum'
      say

      'set dfile 1'
      'query file'
      rec5 = sublin(result,5)
      nxt = subwrd(rec5,3) ; nyt = subwrd(rec5,6) ; nzt = subwrd(rec5,9)
      'set x 1 'nxt
      'set y 1 'nyt
      'set z 1 'nzt
      'set t 1'
      'define tsum = tsum + 'var_name_temp

      'close 1'
      '!rm 'ctl_file_scale

      y=y+1
    endwhile

    say
    say 'calculate July time average temperature from yearly sum...'
    say
*
*  Make local ctl file with scaled vertical levels...
*
    ctl_file = ''dir.m'/'m1'/'y1'/'ctl_file_temp
    ctl_file_scale = ''var_name_temp'.'m1'.'y1'.scale.ctl'
    scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

    say 'open GrADS file for 'm1', 'y1
    say '...'ctl_file_scale
    'open 'ctl_file_scale
    say 'GrADS file ready'

    'set dfile 1'
    'query file'
    rec5 = sublin(result,5)
    nxt = subwrd(rec5,3) ; nyt = subwrd(rec5,6) ; nzt = subwrd(rec5,9)
    'set x 1 'nxt
    'set y 1 'nyt
    'set z 1 'nzt
    'set t 1'
    'define tav = tsum/'nyears

    say
    say 'calculate zonal averages of the July time average...'
    say

    'set x 1'
    'set y 1 'nyt
    'set z 1 'nzt
    'set t 1'
    'define ztav = ave(tav,x=1,x='nxt')'

    say
    say 'plot zonal average of the July time average...'
    say

    'set vpage 'pvar_vpage.m
    'set parea 'pvar_parea
    'set grads off'
    'set gxout contour'
    'set clab off'
    'set ccolor 0'
    'set xlab off'
    'set ylab off'
    'd ztav'

    say
    say 'find and save location of minimum temperature...'
    say

    j=1
    while ( j <= nyt )
      'set y 'j
      'set gxout print'
      'd minloc(ztav,z=1,z='nzt')' ; rec3 = sublin(result,3) ; zmin = subwrd(rec3,1)
*      say result
      say '   ...j, zmin = 'j', 'zmin
      'query gr2xy 'j' 'zmin
      rec = sublin(result,1) ; xloc.j = subwrd(rec,3) ; yloc.j = subwrd(rec,6)
*      say result
      say '      xloc, yloc = 'xloc.j', 'yloc.j
      j=j+1
    endwhile

    'undefine tav' ; 'undefine ztav'
    'close 1'
    '!rm 'ctl_file_scale

  endif

* ----------------------------------------------------------------------

  say
  say 'Calculate and plot zonal and time averaged "'var_name'"...'
  say
  say '...start time average sum at year = 'y1
  say
*
*  Make local ctl file with scaled vertical levels...
*
  ctl_file = ''dir.m'/'m1'/'y1'/'ctl_file.m
  ctl_file_scale = ''var_name'.'m1'.'y1'.scale.ctl'
  scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

  say 'open GrADS file for 'm1', 'y1
  say '...'ctl_file_scale
  'open 'ctl_file_scale
  say 'GrADS file ready'

  'set dfile 1'
  'query file'
  rec5 = sublin(result,5)
  nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)
  if( plot_min_temp = 'yes' )
    if( nx != nxt )
      say '...nxt = 'nxt
      say '...nx  = 'nx
      exit
    endif
    if( ny != nyt )
      say '...nyt = 'nyt
      say '...ny  = 'ny
      exit
    endif
    if( nz != nzt )
      say '...nzt = 'nzt
      say '...nz  = 'nz
      exit
    endif
  endif

  'set x 1 'nx
  'set y 1 'ny
  'set z 1 'nz
  'set t 1'
  'define tsum = 'var_name

  'close 1'
  '!rm 'ctl_file_scale

  say
  say '...continue with running sum of July time average for remaining years'
  say '   (number of years = 'nyears')'
  say

  y=y1+1
  while ( y <= y2 )

    say
    say '*** 'y' ***'
    say

    ctl_file = ''dir.m'/'m1'/'y'/'ctl_file.m
    ctl_file_scale = ''var_name'.'m1'.'y'.scale.ctl'
    scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

    say 'open GrADS file for 'm1', 'y
    say '...'ctl_file_scale
    'open 'ctl_file_scale
    say 'GrADS file ready'

    say
    say '...increment time average sum'
    say

    'set dfile 1'
    'query file'
    rec5 = sublin(result,5)
    nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)
    'set x 1 'nx
    'set y 1 'ny
    'set z 1 'nz
    'set t 1'
    'define tsum = tsum + 'var_name

    'close 1'
    '!rm 'ctl_file_scale

    y=y+1
  endwhile

  say
  say 'calculate July time average from yearly sum...'
  say
*
*  Make local ctl file with scaled vertical levels...
*
  ctl_file = ''dir.m'/'m1'/'y1'/'ctl_file.m
  ctl_file_scale = ''var_name'.'m1'.'y1'.scale.ctl'
  scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

  say 'open GrADS file for 'm1', 'y1
  say '...'ctl_file_scale
  'open 'ctl_file_scale
  say 'GrADS file ready'

  'set dfile 1'
  'query file'
  rec5 = sublin(result,5)
  nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)
  'set x 1 'nx
  'set y 1 'ny
  'set z 1 'nz
  'set t 1'
  'define tav = 'var_scale'*tsum/'nyears

  say
  say 'calculate zonal averages of the July time average...'
  say

  'set x 1'
  'set y 1 'ny
  'set z 1 'nz
  'set t 1'
  'define ztav = ave(tav,x=1,x='nx')'

  say
  say 'plot zonal and time averaged "'var_name'"...'
  say

  'set vpage 'pvar_vpage.m
  'set parea 'pvar_parea
  'set grads off'
  'set clevs '_vcl.m1
  'set ccols '_vcc.m1
  'set gxout shaded'
  'set xlab on'
  'set ylab on'
  'set ylint 'ylabel_int
  'set x 1'
  'set y 1 'ny
  'set lev 'zlevs
  'set t 1'
  'd ztav'
  'run cbar.gs 0.7'
  'draw title \\\ 'model.m
*  'draw xlab Latitude \\\'
  'draw ylab \\\ 'ylabel
  xfp = xl_pvar_parea + 0.1
  yfp = yt_pvar_parea + 0.2
  'set string 1 l 6'
  'set strsiz 0.13'
  'draw string 'xfp' 'yfp' 'figure_panel.m
  'set gxout stat'
  'd ztav'
  min_max = sublin(result,8)
  say '   'min_max
  if( plot_min_max = 'yes' )
    xmm = xl_pvar_parea + 0.1
    ymm = yt_pvar_parea - 0.1
    'set string 1 l 4'
    'set strsiz 0.08'
    'draw string 'xmm' 'ymm' Data: 'min_max' '
  endif

  if( plot_min_temp = 'yes' )
    say '   plot level of minimum temperature...'
    j=1
    while ( j <= ny )
      'set y 'j
*       say '      xloc, yloc = 'xloc.j', 'yloc.j
      'set line 1'
      'draw mark 3 'xloc.j' 'yloc.j' 0.03'
      j=j+1
    endwhile
  endif

  'undefine tsum' ; 'undefine tav'
  'close 1'
  '!rm 'ctl_file_scale

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

* ----------------------------------------------------------------------

  say
  say 'Calculate and plot standard deviation for zonally averaged "'var_name'"...'
  say
  say '...start sample variance sum at year = 'y1
  say
*
*  Make local ctl file with scaled vertical levels...
*
  ctl_file = ''dir.m'/'m1'/'y1'/'ctl_file.m
  ctl_file_scale = ''var_name'.'m1'.'y1'.scale.ctl'
  scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

  say 'open GrADS file for 'm1', 'y1
  say '...'ctl_file_scale
  'open 'ctl_file_scale
  say 'GrADS file ready'

  'set dfile 1'
  'query file'
  rec5 = sublin(result,5)
  nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)

  say
  say '...scale July time average'
  say

  'set x 1 'nx
  'set y 1 'ny
  'set z 1 'nz
  'set t 1'
  'define tav1 = 'var_scale'*'var_name

  say
  say '...calculate zonal average of July time average'
  say

  'set x 1'
  'set y 1 'ny
  'set z 1 'nz
  'set t 1'
  'define ztav1 = ave(tav1,x=1,x='nx')'

  say
  say '...start sample variance sum of zonally-averaged July time average'
  say

  'define svsum = pow((ztav1-ztav),2)'

  'undefine tav1' ; 'undefine ztav1'
  'close 1'
  '!rm 'ctl_file_scale

  say
  say '...running sum for sample variance for remaining years...'
  say '   (number of years = 'nyears')'
  say

  y=y1+1
  while ( y <= y2 )

    say
    say '*** 'y' ***'
    say

    ctl_file = ''dir.m'/'m1'/'y'/'ctl_file.m
    ctl_file_scale = ''var_name'.'m1'.'y'.scale.ctl'
    scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

    say 'open GrADS file for 'm1', 'y
    say '...'ctl_file_scale
    'open 'ctl_file_scale
    say 'GrADS file ready'

    'set dfile 1'
    'query file'
    rec5 = sublin(result,5)
    nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)

    say
    say '...scale July time average'

    'set x 1 'nx
    'set y 1 'ny
    'set z 1 'nz
    'set t 1'
    'define tav1 = 'var_scale'*'var_name

    say '...calculate zonal average of July time average'

    'set x 1'
    'set y 1 'ny
    'set z 1 'nz
    'set t 1'
    'define ztav1 = ave(tav1,x=1,x='nx')'

    say
    say '...increment sample variance sum'
    say

    'define svsum = svsum + pow((ztav1-ztav),2)'

    'undefine tav1' ; 'undefine ztav1'
    'close 1'
    '!rm 'ctl_file_scale

    y=y+1
  endwhile

  say
  say 'calculate sample standard deviation of zonally averaged "'var_name'"...'
  say
*
*  Make local ctl file with scaled vertical levels...
*
  ctl_file = ''dir.m'/'m1'/'y1'/'ctl_file.m
  ctl_file_scale = ''var_name'.'m1'.'y1'.scale.ctl'
  scale_zdef_zlevs(ctl_file,m_to_km,ctl_file_scale)

  say 'open GrADS file for 'm1', 'y1
  say '...'ctl_file_scale
  'open 'ctl_file_scale
  say 'GrADS file ready'

  'set dfile 1'
  'query file'
  rec5 = sublin(result,5)
  nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)
  'set x 1'
  'set y 1 'ny
  'set z 1 'nz
  'set t 1'
  'define sd = sqrt(svsum/'nyears-1')'

  say
  say 'plot sample standard deviation...'
  say

  'set vpage 'pvar_vpage.m
  'set parea 'pvar_parea
  'set grads off'
  'set clab off'
  if( sd_label = 'on' )
    'set clab on'
  endif
  n=1 ; sd_clev = subwrd(sd_level,n) ; sd_cstyle = subwrd(sd_style,n)
  while ( sd_clev != '' )
    say 'n, sd_clev = 'n', 'sd_clev
    'set gxout contour'
    'set clevs 'sd_clev
    'set cstyle 'sd_cstyle
    'set ccolor 1'
    'set xlab off'
    'set ylab off'
    'set x 1'
    'set y 1 'ny
    'set lev 'zlevs
    'set t 1'
    'd sd'
    if( sd_label = 'manual' )
      say
      say 'Contour level = 'sd_clev
      say 'Do you want to plot a contour label?'
      say 'Type y(es) or Y(es), or press "ENTER" to continue...'
      pull plot
      if( plot='y' | plot='yes' | plot='Y' | plot='Yes' )
        clabel(csiz,cformat,sd_clev)
      endif
      say
      say 'Do you want to replot panel?'
      say 'Type y(es) or Y(es), or press "ENTER" to continue...'
      pull replot
      if( replot='y' | replot='yes' | replot='Y' | replot='Yes' )
        say
        say 'plot zonal and time averaged "'var_name'"...'
        say
        'set gxout shaded'
        'set clevs '_vcl.m1
        'set ccols '_vcc.m1
        'set xlab off'
        'set ylab off'
        'set x 1'
        'set y 1 'ny
        'set lev 'zlevs
        'set t 1'
        'd ztav'
        'set gxout stat'
        'd ztav'
        min_max = sublin(result,8)
        say '   'min_max
        if( plot_min_max = 'yes' )
          xmm = xl_pvar_parea + 0.1
          ymm = yt_pvar_parea - 0.1
          'set string 1 l 4'
          'set strsiz 0.08'
          'draw string 'xmm' 'ymm' Data: 'min_max' '
        endif
        if( plot_min_temp = 'yes' )
          say
          say '   plot level of minimum temperature...'
          say
          j=1
          while ( j <= ny )
            'set y 'j
*           say '      xloc, yloc = 'xloc.j', 'yloc.j
            'set line 1'
            'draw mark 3 'xloc.j' 'yloc.j' 0.03'
            j=j+1
          endwhile
        endif
        n=0
      endif
    endif
    n=n+1
    sd_clev = subwrd(sd_level,n) ; sd_cstyle = subwrd(sd_style,n)
  endwhile
  'set gxout stat'
  'd sd'
  min_max = sublin(result,8)
  say '   'min_max
  if( plot_sd_min_max = 'yes' )
      xmm = xl_pvar_parea + 0.1
    if( plot_min_max = 'yes' )
      ymm = yt_pvar_parea - 0.25
    else
      ymm = yt_pvar_parea - 0.1
    endif
    'set string 1 l 4'
    'set strsiz 0.08'
    'draw string 'xmm' 'ymm' SSD: 'min_max
  endif

  'undefine ztav' ; 'undefine svsum'  ; 'undefine sd'
  'close 1'
  '!rm 'ctl_file_scale

  say
  say '...finished plots for: 'model.m
  say

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

  m=m+1
endwhile

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

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

if( batch_mode = 'yes' )
  quit
endif

return

************************************************************************
*
function define_clevcols(units)
*
* Define "clevs" for all months...
*
if( units = 'vmr' )
  _vcl.Jan = '10 20 30 40 50 60 70 80 90 100'
  _vcl.Feb = '10 15 20 25 30 35 40 45 50 55 60 65'
  _vcl.Mar = '5 10 15 20 25 30 35 40 45 50 55 60'
  _vcl.Apr = '4 8 12 16 20 24 28 32 36 40'
  _vcl.May = '3 6 9 12 15 18 21 24 27 30'
  _vcl.Jun = '3 6 9 12 15 18 21 24 27 30'
*  _vcl.Jul = '5 10 15 20 25 30 35 40 45 50 55'
*  _vcl.Jul = '0.0004 0.002 0.004 0.02 0.04 0.2 0.4 2 4 8 10 12 16 20 40'
  _vcl.Jul = '0.002 0.02 0.2 2 4 6 8 10 12 16 20 40'
  _vcl.Aug = '5 10 15 20 25 30 35 40 45 50 55'
  _vcl.Sep = '5 10 15 20 25 30 35 40 45 50 55 60'
  _vcl.Oct = '5 10 15 20 25 30 35 40 45 50 55 60'
  _vcl.Nov = '10 20 30 40 50 60 70'
  _vcl.Dec = '10 20 30 40 50 60 70'
endif
if( units = 'mBq-SCM' )
  _vcl.Jan = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Feb = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Mar = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Apr = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.May = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Jun = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Jul = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Aug = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Sep = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Oct = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Nov = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
  _vcl.Dec = '0.01 0.05 0.1 0.5 1 5 10 50 100 200 250 300 400 500 1000'
endif
*
*  Colors
*
*rainbow_ccols = '9  14  4  11  5  13  3  10  7  12  8  2  6'
define_colors()

*_vcc.Jul = '34  33  32  31  30  25  24  23  22  21  7  12  8  41  42  43'
_vcc.Jul = '34    32     31   30  25   21  18  7  12   8   41   42   43'
*
return
*
************************************************************************

************************************************************************
*
function ptitle_top(title,subtitle,plim)
*
* Plots a title at the top of the page
*
x1 = subwrd(plim,1)
x2 = subwrd(plim,2)
xc = x1 + (x2-x1)/2
'c'
'set vpage 0 8.5 0 11'
'set string 1 c 6'
'set strsiz 0.14'
'draw string 'xc' 10.4 'title
'draw string 'xc' 10.1 'subtitle
*
return
*
************************************************************************

************************************************************************
*
function define_dims(nfile)
*
* defines dimensions for verification purposes...
*
 say 'function define_dims...'
 'set dfile 'nfile
*
* save initial settings...
*
 'query dims'
 rec2 = sublin(result,2) ; range = subwrd(rec2,3)
 if( range = varying )
   lon_setting = ''subwrd(rec2,6)' 'subwrd(rec2,8)''
 else
   if( range = fixed )
     lon_setting = subwrd(rec2,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'lon_setting = 'lon_setting
 rec3 = sublin(result,3) ; range = subwrd(rec3,3)
 if( range = varying )
   lat_setting = ''subwrd(rec3,6)' 'subwrd(rec3,8)''
 else
   if( range = fixed )
     lat_setting = subwrd(rec3,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'lat_setting = 'lat_setting
 rec4 = sublin(result,4) ; range = subwrd(rec4,3)
 if( range = varying )
   lev_setting = ''subwrd(rec4,6)' 'subwrd(rec4,8)''
 else
   if( range = fixed )
     lev_setting = subwrd(rec4,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'lev_setting = 'lev_setting
 rec5 = sublin(result,5) ; range = subwrd(rec5,3)
 if( range = varying )
   time_setting = ''subwrd(rec5,6)' 'subwrd(rec5,8)''
 else
   if( range = fixed )
     time_setting = subwrd(rec5,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'time_setting = 'time_setting
*
* define dimensions...
*
 'query file'
 rec5 = sublin(result,5)
 _nx1 = subwrd(rec5,3) ; _ny1 = subwrd(rec5,6) ; _nz1 = subwrd(rec5,9)
 'set x 1 '_nx1
 'set y 1 '_ny1
 'set z 1 '_nz1
 'set t 1'
 'query dims'
 rec2 = sublin(result,2)
 _lonmin1 = subwrd(rec2,6) ; _lonmax1 = subwrd(rec2,8)
 rec3 = sublin(result,3)
 _latmin1 = subwrd(rec3,6) ; _latmax1 = subwrd(rec3,8)
 rec4 = sublin(result,4)
 _levmin1 = subwrd(rec4,6) ; _levmax1 = subwrd(rec4,8)
* say _nx1' '_lonmin1' '_lonmax1
* say _ny1' '_latmin1' '_latmax1
* say _nz1' '_levmin1' '_levmax1
*
* reset lon, lat, lev, time setttings to initial values...
*
 'set lon 'lon_setting
 'set lat 'lat_setting
 'set lev 'lev_setting
 'set time 'time_setting
*
 say '...leave function define_dims'
*
return
*
************************************************************************

************************************************************************
*
function check_dims(nfile)
*
* checks dimensions of current file for verification purposes...
*
 say 'function check_dims...'
 'set dfile 'nfile
*
* save initial settings...
*
 'query dims'
 rec2 = sublin(result,2) ; range = subwrd(rec2,3)
 if( range = varying )
   lon_setting = ''subwrd(rec2,6)' 'subwrd(rec2,8)''
 else
   if( range = fixed )
     lon_setting = subwrd(rec2,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'lon_setting = 'lon_setting
 rec3 = sublin(result,3) ; range = subwrd(rec3,3)
 if( range = varying )
   lat_setting = ''subwrd(rec3,6)' 'subwrd(rec3,8)''
 else
   if( range = fixed )
     lat_setting = subwrd(rec3,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'lat_setting = 'lat_setting
 rec4 = sublin(result,4) ; range = subwrd(rec4,3)
 if( range = varying )
   lev_setting = ''subwrd(rec4,6)' 'subwrd(rec4,8)''
 else
   if( range = fixed )
     lev_setting = subwrd(rec4,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'lev_setting = 'lev_setting
 rec5 = sublin(result,5) ; range = subwrd(rec5,3)
 if( range = varying )
   time_setting = ''subwrd(rec5,6)' 'subwrd(rec5,8)''
 else
   if( range = fixed )
     time_setting = subwrd(rec5,6)
   else
     say ' range = 'range
     exit
   endif
 endif
* say 'time_setting = 'time_setting
*
* check dimensions...
*
 'query file'
 rec5 = sublin(result,5)
 nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9)
 'set x 1 'nx
 'set y 1 'ny
 'set z 1 'nz
 'set t 1'
 'query dims'
 rec2 = sublin(result,2)
 lonmin = subwrd(rec2,6) ; lonmax = subwrd(rec2,8)
 rec3 = sublin(result,3)
 latmin = subwrd(rec3,6) ; latmax = subwrd(rec3,8)
 rec4 = sublin(result,4)
 levmin = subwrd(rec4,6) ; levmax = subwrd(rec4,8)
 if( nx != _nx1 )
   say '_nx1 = '_nx1', nx = 'nx
   exit
 endif
 if( lonmin != _lonmin1 )
   say '_lonmin1 = '_lonmin1', lonmin = 'lonmin
   exit
 endif
 if( lonmax != _lonmax1 )
   say '_lonmax1 = '_lonmax1', lonmax = 'lonmax
   exit
 endif
 if( ny != _ny1 )
   say '_ny1 = '_ny1', ny = 'ny
   exit
 endif
 if( latmin != _latmin1 )
   say '_latmin1 = '_latmin1', latmin = 'latmin
   exit
 endif
 if( latmax != _latmax1 )
   say '_latmax1 = '_latmax1', latmax = 'latmax
   exit
 endif
 if( nz != _nz1 )
   say '_nz1 = '_nz1', nz = 'nz
   exit
 endif
 if( levmin != _levmin1 )
   say '_levmin1 = '_levmin1', levmin = 'levmin
   exit
 endif
 if( levmax != _levmax1 )
   say '_levmax1 = '_levmax1', levmax = 'levmax
   exit
 endif
*
* reset lon, lat, lev, time setttings to initial values...
*
 'set lon 'lon_setting
 'set lat 'lat_setting
 'set lev 'lev_setting
 'set time 'time_setting
*
 say '...leave function check_dims'
*
return
*
************************************************************************
