*-----------------------------------------------------------------------
*
*  To use:  "grads -p"
*
*  Script calculates differences for all variables for two GrADS files,
*  at time levels t = (1, tsiz/2, tsiz), where "tsiz" is the number of
*  time levels in the GrADS dataset, as reported in the GrADS ctl file
*  TDEF record...
*
*-----------------------------------------------------------------------
'reinit'

*-----------------------------------------------------------------------
*  allow external GrADS functions...
*-----------------------------------------------------------------------

rc = gsfallow("on")

*-----------------------------------------------------------------------
*  pause between plots?
*-----------------------------------------------------------------------

do_pause = 'yes'

*-----------------------------------------------------------------------
*  check diffs?
*-----------------------------------------------------------------------

do_check_diffs = 'yes'
xinc = 1
yinc = 1
zinc = 31

*-----------------------------------------------------------------------
*  print to GrADS metafile?
*-----------------------------------------------------------------------

do_print = 'no'

*
*  define plot variables' page limits:
*
pv_vpage.1 = '0 8.5  6.75 10.0'
pv_vpage.2 = '0 8.5  3.5   6.75'
pv_vpage.3 = '0 8.5  0.25  3.5'
pv_parea   = '1 7.0  0.5   2.75'

xl = subwrd(pv_parea,1) ; xr = subwrd(pv_parea,2)
xc = xl + (xr-xl)/2
yb = subwrd(pv_parea,3) ; yt = subwrd(pv_parea,4)
yc = yb + (yt-yb)/2

*-----------------------------------------------------------------------
*  define the two directories for comparison
*-----------------------------------------------------------------------

dir1 = '/archive/cjs/zetac/KWAJEX/history/test-do_nccat'
dir2 = '/archive/cjs/zetac/KWAJEX/history/test-do_append'
*dir2 = '/archive/cjs/zetac/KWAJEX/history/test-do_templates'

*-----------------------------------------------------------------------
*  define the "describe dataset" names for each of the GrADS ctl files
*  (includes the "date" string, prepended with any "describe" string as
*  appropriate here):
*-----------------------------------------------------------------------

ctl_dd1 = '06:30z24Jul1999_2'
ctl_dd2 = '06:30z24Jul1999_2'
*ctl_dd2 = '06:30z24jul1999_2'

*-----------------------------------------------------------------------
*  set processor "sh" script to define GrADS ctl root name(s)
*  (to define diagnostic variables)
*-----------------------------------------------------------------------

sh_script = ncl_3d

if( sh_script = ncl_2d )
*
*  Non-radiation...
*
***  ctl_root = 'psurf tsurf rain t_surf q_surf shflx evap qc_int qi_int t_ref drag_heat drag_mom'
*
*  Radiation...
*
  ctl_root = 'swup_toa swdn_toa olr swup_sfc swdn_sfc lwup_sfc lwdn_sfc'
endif

if( sh_script = ncl_3d )
*
*  Non-radiation...
*
***  ctl_root = 'uu_ext vv_ext ex_ext th_ext qv_ext ttnd_ext qtnd_ext uu vv th qv ww pp ta rh qc qr qi qs qg qcond cloud vfall tdt_micro qdt_micro tke_iso tke_diff'
*
*  Radiation...
*
  ctl_root = 'tdt_sw tdt_lw'
endif

if( sh_script = uflux_3d )
  ctl_root = 'UFLUX'
endif

if( sh_script = vflux_3d )
  ctl_root = 'VFLUX'
endif

if( sh_script = massflux_3d )
  ctl_root = 'MASSFLUX'
endif

if( sh_script = QQterms_3d )
  ctl_root = 'QQ'
endif

if( sh_script = QTterms_3d )
  ctl_root = 'QT'
endif

if( sh_script = cldfr_3d )
  ctl_root = 'CLDFR'
endif

if( sh_script = regrid )
  ctl_root = 'REGRID'
endif

if( sh_script = skewt )
  ctl_root = 'SKEWT'
endif

if( sh_script = mphysics )
  ctl_root = 'MPHYSICS'
endif

if( sh_script = path )
  ctl_root = 'PATH'
endif

if( sh_script = qc_spectrum )
  ctl_root = 'qc_spec'
endif

if( sh_script = qi_spectrum )
  ctl_root = 'qi_spec'
endif

if( sh_script = qr_spectrum )
  ctl_root = 'qr_spec'
endif

if( sh_script = qs_spectrum )
  ctl_root = 'qs_spec'
endif

if( sh_script = qg_spectrum )
  ctl_root = 'qg_spec'
endif

if( sh_script = qt_spectrum )
  ctl_root = 'qt_spec'
endif

if( sh_script = w_spectrum )
  ctl_root = 'w_spec_GATE'
endif

if( sh_script = w_spectrum_z )
  ctl_root = 'w_spec_z_KWAJEX'
endif

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

say
say 'Comparing:'
say
say '*****************************************************************'
say
say 'Directory 1: 'dir1
say 'GrADS dataset describe string: 'ctl_dd1
say
say 'Directory 2: 'dir2
say 'GrADS dataset describe string: 'ctl_dd2
say
say ctl_root
say
say '*****************************************************************'
say

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

*-----------------------------------------------------------------------
*  define number of GrADS ctl files and dataset name variables...
*-----------------------------------------------------------------------

nctl = number( ctl_root )
n=1
while ( n <= nctl )
  grd1.n = ''subwrd(ctl_root,n)'.'ctl_dd1
*  say grd1.n
  grd2.n = ''subwrd(ctl_root,n)'.'ctl_dd2
*  say grd2.n
  n=n+1
endwhile

*-----------------------------------------------------------------------
*  dmget GrADS files
*-----------------------------------------------------------------------
*
*  Note, we use a "*" in the "ls" command to list "template" files, if
*  any exist...
*
say
say 'Define GrADS files to dmget...'
say
'!rm ieee_file1_list ieee_file2_list'
'!touch ieee_file1_list ieee_file2_list'
say
say 'GrADS ieee file(s):'
n=1
while ( n <= nctl )
  grads_ieee = ''dir1'/'grd1.n'*.ieee'
  say '...'grads_ieee
  '!ls -x 'grads_ieee' >> 'ieee_file1_list
  grads_ieee = ''dir2'/'grd2.n'*.ieee'
  say '...'grads_ieee
  '!ls -x 'grads_ieee' >> 'ieee_file2_list
  n=n+1
endwhile

'!rm ctl_file1_list ctl_file2_list'
'!touch ctl_file1_list ctl_file2_list'
say
say 'GrADS ctl file(s):'
n=1
while ( n <= nctl )
  grads_ctl1.n = ''dir1'/'grd1.n'.ctl'
  say '...'grads_ctl1.n
  '!ls -x 'grads_ctl1.n' >> 'ctl_file1_list
  grads_ctl2.n = ''dir2'/'grd2.n'.ctl'
  say '...'grads_ctl2.n
  '!ls -x 'grads_ctl2.n' >> 'ctl_file2_list
  n=n+1
endwhile

say
say 'dmget GrADS ieee files...'
'!dmget `cat ieee_file1_list ieee_file2_list`'
say '...GrADS model ieee files ready.'
say 'dmget GrADS ctl files...'
'!dmget `cat ctl_file1_list ctl_file2_list`'
say '...GrADS model ctl files ready.'
'!rm ieee_file1_list ieee_file2_list ctl_file1_list ctl_file2_list'

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

nf=1
*-----------------------------------------------------------------------
while ( nf <= nctl )
*-----------------------------------------------------------------------
  say
  say '################################################################'
  say ' Plotting GrADS data for processor "'sh_script'"'
  say ' for GrADS ctl files:'
  say ' 'grads_ctl1.nf
  say ' 'grads_ctl2.nf
  say '################################################################'
  say
  'open 'grads_ctl1.nf
  'open 'grads_ctl2.nf
  say
  say ' Characteristics of GrADS datasets:'
  say
  say '************************************************'
  say
  'set dfile 1'
  'q ctlinfo'
    rec = sublin(result,3); undef1 = subwrd(rec,2)
  'q file'
    say result
    rec5 = sublin(result,5)
    xsiz1 = subwrd(rec5,3); ysiz1 = subwrd(rec5,6); zsiz1 = subwrd(rec5,9); tsiz1 = subwrd(rec5,12)
    rec6 = sublin(result,6)
    nvar1 = subwrd(rec6,5)
    nv=1
    while ( nv <= nvar1 )
      rec = sublin(result,nv+6)
      var1.nv = subwrd(rec,1)
      nv=nv+1
    endwhile
  say
  'set dfile 2'
  'q ctlinfo'
    rec = sublin(result,3); undef2 = subwrd(rec,2)
  'q file'
    say result
    rec5 = sublin(result,5)
    xsiz2 = subwrd(rec5,3); ysiz2 = subwrd(rec5,6); zsiz2= subwrd(rec5,9); tsiz2 = subwrd(rec5,12)
    rec6 = sublin(result,6)
    nvar2 = subwrd(rec6,5)
    nv=1
    while ( nv <= nvar2 )
      rec = sublin(result,nv+6)
      var2.nv = subwrd(rec,1)
      nv=nv+1
    endwhile
  say
  say '************************************************'
  say
*
*  Verify dataset info for files...
*
  if( undef2 != undef1 )
    say ' undef1 = 'undef1
    say ' undef2 = 'undef2
    say ' Exiting script.'
    exit
  else
    undef = undef1
  endif

  if( xsiz2 != xsiz1 )
    say ' xsiz1 = 'xsiz1
    say ' xsiz2 = 'xsiz2
    say ' Exiting script.'
    exit
  else
    xsiz = xsiz1
  endif
  if( ysiz2 != ysiz1 )
    say ' ysiz1 = 'ysiz1
    say ' ysiz2 = 'ysiz2
    say ' Exiting script.'
    exit
  else
    ysiz = ysiz1
  endif
  if( zsiz2 != zsiz1 )
    say ' zsiz1 = 'zsiz1
    say ' zsiz2 = 'zsiz2
    say ' Exiting script.'
    exit
  else
    zsiz = zsiz1
  endif
  if( tsiz2 != tsiz1 )
    say ' tsiz1 = 'tsiz1
    say ' tsiz2 = 'tsiz2
    say ' Exiting script.'
    exit
  else
    tsiz = tsiz1
  endif
  if( nvar2 != nvar1 )
    say ' nvar1 = 'nvar1
    say ' nvar2 = 'nvar2
    say ' Exiting script.'
    exit
  else
    nvar = nvar1
    nv=1
    while ( nv <= nvar )
      if( var2.nv != var1.nv )
        say ' nvar2 != nvar1'
        say ' Exiting script.'
        exit
      else
        var.nv = var1.nv
      endif
      nv=nv+1
    endwhile
  endif
*
*  Further comparisons on "starting" and "final" longitude, latitude,
*  levels, and times...
*
  'set dfile 1'
  'set x 1 'xsiz
  'set y 1 'ysiz
  'set z 1 'zsiz
  'set t 1 'tsiz
  'q dims'
  rec2 = sublin(result,2) ; lon1s = subwrd(rec2,6) ; lon1f = subwrd(rec2,8)
  rec3 = sublin(result,3) ; lat1s = subwrd(rec3,6) ; lat1f = subwrd(rec3,8)
  rec4 = sublin(result,4) ; lev1s = subwrd(rec4,6) ; lev1f = subwrd(rec4,8)
  rec5 = sublin(result,5) ; time1s = subwrd(rec5,6) ; time1f = subwrd(rec5,8)
  'set dfile 2'
  'set x 1 'xsiz
  'set y 1 'ysiz
  'set z 1 'zsiz
  'set t 1 'tsiz
  'q dims'
  rec2 = sublin(result,2) ; lon2s = subwrd(rec2,6) ; lon2f = subwrd(rec2,8)
  rec3 = sublin(result,3) ; lat2s = subwrd(rec3,6) ; lat2f = subwrd(rec3,8)
  rec4 = sublin(result,4) ; lev2s = subwrd(rec4,6) ; lev2f = subwrd(rec4,8)
  rec5 = sublin(result,5) ; time2s = subwrd(rec5,6) ; time2f = subwrd(rec5,8)

  if( lon2s != lon1s )
    say ' lon1s = 'lon1s
    say ' lon2s = 'lon2s
    say ' Exiting script.'
    exit
  endif
  if( lon2f != lon1f )
    say ' lon1f = 'lon1f
    say ' lon2f = 'lon2f
    say ' Exiting script.'
    exit
  endif
  if( lat2s != lat1s )
    say ' lat1s = 'lat1s
    say ' lat2s = 'lat2s
    say ' Exiting script.'
    exit
  endif
  if( lat2f != lat1f )
    say ' lat1f = 'lat1f
    say ' lat2f = 'lat2f
    say ' Exiting script.'
    exit
  endif
  if( lev2s != lev1s )
    say ' lev1s = 'lev1s
    say ' lev2s = 'lev2s
    say ' Exiting script.'
    exit
  endif
  if( lev2f != lev1f )
    say ' lev1f = 'lev1f
    say ' lev2f = 'lev2f
    say ' Exiting script.'
    exit
  endif
  if( time2s != time1s )
    say ' time1s = 'time1s
    say ' time2s = 'time2s
    say ' Exiting script.'
    exit
  endif
  if( time2f != time1f )
    say ' time1f = 'time1f
    say ' time2f = 'time2f
    say ' Exiting script.'
    exit
  endif
*
*  define time levels for difference checks...
*
  nt = 3
  tt.1 = 1 ; tt.2 = tsiz/2 ; tt.3 = tsiz
  say
  say 'Time levels for difference checks:'
  say
  say ' tt.1 = 'tt.1', tt.2 = 'tt.2', tt.3 = 'tt.3
  say
*
*  plot t/z sections for each of the variable(s)...
*
  if( do_print = 'yes' )
    'enable print 'sh_script'.gx'
  endif
*
*  plot title page...
*
  ptitle_page(grads_ctl1.nf,grads_ctl2.nf)

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

  nv=1
*-----------------------------------------------------------------------
  while ( nv <= nvar )
*-----------------------------------------------------------------------

    var = var.nv
*
*  plot title at top of page...
*
    ptitle_top(xc,var)

    say
    say '- - - - - - - - - - - - - - - - - - -'
    say 'Plot variable "'var'" from File 1:'
    say '"'grads_ctl1.nf'"'
    say '- - - - - - - - - - - - - - - - - - -'
    say
    'set dfile 1'
    'set vpage  0 8.5  6.75 10.0'
    'set parea  1 7.0  0.5  2.75'
    'set x 1'
    'set y 1'
    'set z 1 'zsiz
    'set t 1 'tsiz
    'set gxout shaded'
    'set grads off'
    'd 'var
    'run cbar.gs 0.7'
    'draw title \\\ File 1'
*    'draw xlab Time \\\'
    'draw ylab \\\ Height (km)'
    say
    say '- - - - - - - - - - - - - - - - - - -'
    say 'Plot variable "'var'" from File 2:'
    say '"'grads_ctl2.nf'"'
    say '- - - - - - - - - - - - - - - - - - -'
    say
    'set dfile 2'
    'set vpage  0 8.5  3.5 6.75'
    'set parea  1 7.0  0.5 2.75'
    'set x 1'
    'set y 1'
    'set z 1 'zsiz
    'set t 1 'tsiz
    'set gxout shaded'
    'set grads off'
    'd 'var
    'run cbar.gs 0.7'
    'draw title \\\ File 2'
*    'draw xlab Time \\\'
    'draw ylab \\\ Height (km)'
    say
    say '- - - - - - - - - - - - - - - - - - -'
    say 'Plot difference field for: 'var
    say '- - - - - - - - - - - - - - - - - - -'
    say
    'set vpage  0 8.5  0.25 3.5'
    'set parea  1 7.0  0.5  2.75'
    'set x 1'
    'set y 1'
    'set z 1 'zsiz
    'set t 1 'tsiz
    'set gxout shaded'
    'set grads off'
    'd 'var'.2-'var'.1' ; diff_result = sublin(result,1)
    if( diff_result = 'Constant field.  Value = 0' )
*
*  plot "diff_result" string on plot...
*
      pdiff_result(xc,yc,diff_result)
    endif
    'draw title \\\ File 2 - File 1'
*    'draw xlab Time \\\'
    'draw ylab \\\ Height (km)'
    'run cbar.gs 0.7'

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

    if( do_check_diffs = 'yes' )
      say
      say '- - - - - - - - - - - - - - - - - - -'
      say 'Check difference field for: 'var
      say '- - - - - - - - - - - - - - - - - - -'
      say
      n=1
      while ( n <= nt )
        tt = tt.n
        say '----------'
        say ' tt = 'tt
        say '----------'
***        pause()
        zz=1
      while ( zz <= zsiz )
        say ' zz = 'zz
***        pause()
        yy=1
      while ( yy <= ysiz )
        xx=1
      while ( xx <= xsiz )
*        say ' yy, xx = 'yy', 'xx
        'set x 1'
        'set y 1'
        'set z 1'
        'set t 1'
        'd 'var'.1(x='xx',y='yy',z='zz',t='tt')'
        rec = sublin(result,1); var1_value = subwrd(rec,4)
        'd 'var'.2(x='xx',y='yy',z='zz',t='tt')'
        rec = sublin(result,1); var2_value = subwrd(rec,4)
        diff = var2_value - var1_value
        if( diff != 0 )
          say ' Non-zero difference = 'diff
          say '   xx = 'xx
          say '   "var1" value = 'var1_value
          say '   "var2" value = 'var2_value
*          say ' Exiting script.'
*          exit
        endif
        'c events'
        xx = xx+xinc
      endwhile
        yy = yy+yinc
      endwhile
        zz = zz+zinc
      endwhile
        n = n+1
      endwhile
    endif

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

    nv=nv+1
*-----------------------------------------------------------------------
  endwhile
*-----------------------------------------------------------------------

  'close 2' ; 'close 1'

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

  nf=nf+1
*-----------------------------------------------------------------------
endwhile
*-----------------------------------------------------------------------

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

************************************************************************
*
function ptitle_page(file1,file2)
*
* Plots a title page
*
'c'
'set vpage 0 8.5 0 11'
'set string 1 c 6'
'set strsiz 0.11'
'draw string 4.25 7.5 File 1:'
'draw string 4.25 7.0 'file1
'draw string 4.25 6.0 File 2:'
'draw string 4.25 5.5 'file2
*
return
*
************************************************************************

************************************************************************
*
function ptitle_top(xc,var)
*
* Plots a title at top of page
*
'c'
'set vpage 0 8.5 0 11'
'set string 1 c 6'
'set strsiz 0.14'
'draw string 'xc' 10.5 "'var'"'
*
return
*
************************************************************************

************************************************************************
*
function pdiff_result(xc,yc,diff_result)
*
* Plots "diff_result" string in center of plot...
*
'set string 1 c 6'
'set strsiz 0.11'
'draw string 'xc' 'yc' 'diff_result
*
return
*
************************************************************************
