**********************************************************************
*
*  To use, "grads -p"
*
*  Script calculates differences for all variables for two GCSS 5 LSF
*  datasets, at all height and time levels.
*
**********************************************************************
'reinit'
*
*  allow external GrADS functions...
*
rc = gsfallow("on")
*
*  pause between plots?
*
do_pause = 'no'
*
*  print to GrADS metafile?
*
do_print = 'no'
if( do_print = 'yes' )
  'enable print gcss5_diff.gx'
endif
*
*  files to difference...
*
nc1 = '/net/cjs/gcss/case5/data/subcase_B0/IFA_v2.1/zetac/gcss5-zetac-old.nc'
nc2 = '/net/cjs/gcss/case5/data/subcase_B0/IFA_v2.1/zetac/gcss5-zetac.nc'
*
*  variables...
*
nvars = 8
vars1  = 'sst psfc ta qv u  v  dthdt dqvdt'
vars2  = 'ts  ps   ta qv uu vv thlsf qvlsf'
say
say '-----------------------------------------------------------'
say 'Difference script for GCSS 5 LSF netCDF files...'
say '-----------------------------------------------------------'
say
say ' File 1: 'nc1
say ' File 2: 'nc2
say ' variables to difference:'
say '   File 1 variables: 'vars1
say '   File 2 variables: 'vars2
say
say '-----------------------------------------------------------'
say
'sdfopen 'nc1
'sdfopen 'nc2
say
say ' Characteristics of current Datasets:'
say
say '************************************************'
say
'query file 1'
say result
say
'query file 2'
say result
say '************************************************'
say
*
*  Get dataset size info
*
'query file 1'
rec5 = sublin(result,5)
xsiz1 = subwrd(rec5,3); ysiz1 = subwrd(rec5,6); zsiz1 = subwrd(rec5,9); tsiz1 = subwrd(rec5,12)
say
say ' xsiz1 'xsiz1
say ' ysiz1 'ysiz1
say ' zsiz1 'zsiz1
say ' tsiz1 'tsiz1
say
'query file 2'
rec5 = sublin(result,5)
xsiz2 = subwrd(rec5,3); ysiz2 = subwrd(rec5,6); zsiz2 = subwrd(rec5,9); tsiz2 = subwrd(rec5,12)
say
say ' xsiz2 'xsiz2
say ' ysiz2 'ysiz2
say ' zsiz2 'zsiz2
say ' tsiz2 'tsiz2
say
if( xsiz2 != xsiz1 )
  say ' xsiz2 != xsiz1'
  say ' Exiting script.'
  exit
else
  xsiz = xsiz1
endif
if( ysiz2 != ysiz1 )
  say ' ysiz2 != ysiz1'
  say ' Exiting script.'
  exit
else
  ysiz = ysiz1
endif
if( zsiz2 != zsiz1 )
  say ' zsiz2 != zsiz1'
  say ' Exiting script.'
  exit
else
  zsiz = zsiz1
endif
if( tsiz2 != tsiz1 )
  say ' Warning: tsiz2 != tsiz1'
  say ' tsiz2 = 'tsiz2
  say ' tsiz1 = 'tsiz1
  if( tsiz2 > tsiz1 )
    say ' Using "tsiz1"...'
    tsiz = tsiz1
  else
    say ' Using "tsiz2"...'
    tsiz = tsiz2
  endif
else
  tsiz = tsiz1
endif
say
*
* Plot a title page
*
ptitle(nc1,nc2)

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

ts=1; tf=tsiz

tt=ts
while ( tt <= tf )
  'c'
  say
  say '-------------------------------------'
  say ' Time-level 'tt
  say '-------------------------------------'
  say
  'set t 'tt
  'query time'; rec = sublin(result,1); date = subwrd(rec,3)
  nv=1
  while ( nv <= nvars )
    var1 = subwrd(vars1,nv)
    var2 = subwrd(vars2,nv)
    zz=1
    while ( zz <= zsiz )
      'c'
      say
      say '- - - - - - - - - - - - - - - - - - -'
      say 'Variable 'var1'... from file 1'
      say '- - - - - - - - - - - - - - - - - - -'
      say
      'set dfile 1'
      'set vpage  0  8.5  7.2 10.9'
      'set parea  2  6.5  0.2  3.0'
      'set x 1 'xsiz
      'set y 1 'ysiz
      'set z 'zz
      'define v1 = 'var1'.1'
      'set grads off'
      'set gxout shaded'
      'd v1'
      'draw title \\\ "'var1'.1" at (z,t) = ('zz','tt')'
      say
      say '- - - - - - - - - - - - - - - - - - -'
      say 'Variable 'var2'... from file 2'
      say '- - - - - - - - - - - - - - - - - - -'
      say
      'set dfile 2'
      'set vpage  0  8.5  3.7  7.3'
      'set parea  2  6.5  0.2  3.0'
      'set x 1 'xsiz
      'set y 1 'ysiz
      'set z 'zz
      'define v2 = 'var2'.2'
      'set grads off'
      'set gxout shaded'
      'd v2'
      'draw title \\\ "'var2'.2" at (z,t) = ('zz','tt')'
      say
      say '- - - - - - - - - - - - - - - - - - -'
      say 'Plot differences'
      say '- - - - - - - - - - - - - - - - - - -'
      say
      'set vpage  0  8.5  0.1  3.7'
      'set parea  2  6.5  0.2  3.0'
      'set x 1 'xsiz
      'set y 1 'ysiz
      'set z 'zz
      'set grads off'
      'set gxout shaded'
      'd v2-v1' ; rec = sublin(result,1); diff = subwrd(rec,5)
      'draw title \\\ "'var2'.2-'var1'.1" at (z,t) = ('zz','tt')'
      if( do_print = 'yes' )
        'print'
      endif
      if( diff != 0 )
        say ' Non-zero difference = 'diff
        say ' Exiting script.'
        exit
      endif
      if( do_pause = 'yes' )
        pause()
      endif
      zz = zz+1
    endwhile
    nv = nv+1
  endwhile
  tt = tt+1
endwhile

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

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

************************************************************************
*
function ptitle(file1,file2)
*
* Plots a title page
*
'c'
'set vpage 0 8.5 0 11'
'set string 1 c 6'
'set strsiz 0.16'
'draw string 4.25 7.0 Dataset Identifiers:'
'set strsiz 0.14'
'draw string 4.25 6.0 File 1:'
'draw string 4.25 5.5 'file1
'draw string 4.25 4.5 File 2:'
'draw string 4.25 4.0 'file2
*
return
*
************************************************************************
