*-----------------------------------------------------------
*
* To use, "grads -p"
*
*-----------------------------------------------------------
'reinit'

dir_CM2   = '/archive/ccsp/ipcc_ar4/CM2Q-d2_1861-2000-AllForc_h1/pp/ice_8xdaily_instA3/ts/3hr/5yr'
dir_PCMDI = '/archive/pcmdi/ipcc_ar4/CMOR/CM2Q-d2_1861-2000-AllForc_h1/pp/atmos/ts/3hr/5yr'

nc_CM2   = 'ice_8xdaily_instA3.1996010100-2000123123.LH.nc'
nc_PCMDI = 'hfls_A3.1996010100-2000123123.nc'

tinc=100

do_print = false

npts = 6
PCMDI_lons = ' 66.25    128.75     186.25     186.25    313.75    358.75'
PCMDI_lats = '-11.00     15.00      15.00     -13.00     39.00    -51.00'
CM2_lons   = ' 66.089  -231.2204  -173.8587  -173.8587  -46.29519  -1.225555'
CM2_lats   = '-11.0025   15.19512   15.19512  -13.04515   39.5    -50.5'

clevs = '0 50 100 150 200 250 300 350 400 450 500 550 600'

*
*  Begin writing an ASCII output file
*
*res = write ('hfls_A3-ice.txt','')
*res = write ('hfls_A3-ice.txt',' Surface Upward Latent Heat Flux')
res = write ('hfls_A3-ice.txt','')
res = write ('hfls_A3-ice.txt','     CM2                 PCMDI')
res = write ('hfls_A3-ice.txt','  Lon   Lat            Lon   Lat')
n=1
while ( n <= npts )
  xlon1 = subwrd(CM2_lons,n)
  ylat1 = subwrd(CM2_lats,n)
  xlon2 = subwrd(PCMDI_lons,n)
  ylat2 = subwrd(PCMDI_lats,n)
  if ( xlon1 < 0 )
    xlon1 = xlon1 + 360
  endif
  res = write ('hfls_A3-ice.txt',' 'xlon1' 'ylat1'      'xlon2' 'ylat2'')
  n = n+1
endwhile
res = write ('hfls_A3-ice.txt','')

if ( do_print = true )
  'enable print hfls_A3-ice.gx'
endif

*
* Plot a title page
*
ptitle('CM2 "LH" and PCMDI "hfls" 1996-2000')
if ( do_print = true )
  'print'
endif
pause()

'sdfopen 'dir_CM2'/'nc_CM2
'sdfopen 'dir_PCMDI'/'nc_PCMDI
*
*  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
'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

if ( tsiz2 != tsiz1 )
  say ' tsiz2 != tsiz1 '
  say ' Exiting script.'
  exit
else
  tsiz = tsiz1
endif

ts=1; tf=tsiz

res = write ('hfls_A3-ice.txt','')
res = write ('hfls_A3-ice.txt','-----------------------------------')
res = write ('hfls_A3-ice.txt','CM2 "LH" and PCMDI "hfls" 1996-2000')
res = write ('hfls_A3-ice.txt','-----------------------------------')
res = write ('hfls_A3-ice.txt','')

tt=ts
while ( tt <= tf )
  'c'
  say ' '
  say ' Time-level 'tt

  'set dfile 1'
  'set vpage  0  8.5  6.75 10.0'
  'set parea  1  7.0  0.5  2.75'
  'set lon -280 80'
  'set y 1 'ysiz1
  'set z 1 '
  'set t 'tt
  'query time'; rec = sublin(result,1); date = subwrd(rec,3)
  'set grads off'
  'set gxout shaded'
  'set clevs 'clevs
  'd LH'
  'draw title \\\ CM2 "LH" (W m`a-2`n) at time = 'date
  'run cbar 0.7'

  'set dfile 2'
  'set vpage  0  8.5  3.5  6.75'
  'set parea  1  7.0  0.5  2.75'
  'set lon -280 80'
  'set y 1 'ysiz2
  'set z 1 '
  'set t 'tt
  'query time'; rec = sublin(result,1); date = subwrd(rec,3)
  'set grads off'
  'set gxout shaded'
  'set clevs 'clevs
  'd hfls'
  'draw title \\\ PCMDI "hfls" (W m`a-2`n) at time = 'date
  'run cbar 0.7'

  n = 1
  yl = 2.75
  res = write ('hfls_A3-ice.txt','   Time = 'date'')
  res = write ('hfls_A3-ice.txt','     CM2       PCMDI')
  while ( n <= npts )
    xlon1 = subwrd(CM2_lons,n)
    ylat1 = subwrd(CM2_lats,n)
    'set dfile 1'
    'set lon 'xlon1
    'set lat 'ylat1
    'set z 1 '
    'set t 'tt
    'set grads off'
    'd LH.1';        rec = sublin(result,1);   cm2 = subwrd(rec,4)
    xlon2 = subwrd(PCMDI_lons,n)
    ylat2 = subwrd(PCMDI_lats,n)
    'set dfile 2'
    'set lon 'xlon2
    'set lat 'ylat2
    'set z 1 '
    'set t 'tt
    'set grads off'
    'd hfls.2';      rec = sublin(result,1); pcmdi = subwrd(rec,4)
      say ' '
      say '  "CM2 lon" = 'xlon1
      say '  "CM2 lat" = 'ylat1
      say '  "LH"   = 'cm2
      say '  "CM2 lon"   = 'xlon2
      say '  "CM2 lat"   = 'ylat2
      say '  "hfls" = 'pcmdi
    pflux(yl,xlon1,ylat1,cm2,xlon2,ylat2,pcmdi)
    res = write ('hfls_A3-ice.txt','    'cm2'   'pcmdi'')
    yl = yl - 0.4
    n = n+1
  endwhile
  res = write ('hfls_A3-ice.txt','')

  if ( do_print = true )
    'print'
  endif

***  pause()

  tt = tt+tinc
endwhile

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

say ' '
say '************************************************'
say ' '
say ' Finished with this script.'
say ' '
say '************************************************'
say ' '
*
************************************************************************
*
function pause()
*
* Pauses till a key is pressed
*
 say ' '
 say '*********************************'
 say 'Press ENTER key to continue...'
 say '*********************************'
 pull answer
*
return
*
************************************************************************

************************************************************************
*
function ptitle(title)
*
* Plots a title page
*
'c'
'set vpage   0 8.5 0 11'
'set string 1 c 6'
'set strsiz 0.20'
'draw string 4.25 6.0 'title
*
return
*
************************************************************************

************************************************************************
*
function pflux(yl,x1,y1,cm2,x2,y2,pcmdi)
*
* Plots "CM2", "PCMDI" flux info
*
'set vpage  0  8.5  0    3.5'
'set string 1 l 6'
'set strsiz 0.10'
'draw string 2 'yl' CM2   lon: 'x1' lat: 'y1' flux: 'cm2
'draw string 2 'yl-0.15' PCMDI lon: 'x2' lat: 'y2' flux: 'pcmdi
*
return
*
************************************************************************
