function maxwind2(args)
  'reinit'
  wind1_nc_file = subwrd(args,1) ; wind1_var = subwrd(args,2)
  wind2_nc_file = subwrd(args,3) ; wind2_var = subwrd(args,4)
  wind_units = subwrd(args,5)
  lev_units = subwrd(args,6)
  say
  say wind1_nc_file
  say wind1_var
  say wind2_nc_file
  say wind2_var
  say 'wind units: ('wind_units')'
  say ' lev units: ('lev_units')'
  say

* open "wind1_nc_file" and "wind2_nc_file" and
* define variables for (x,y,z,t) dimension limits

  'sdfopen 'wind1_nc_file
  'q file 1'
  rec5 = sublin(result,5)
  nx1 = subwrd(rec5,3) ; ny1 = subwrd(rec5,6) ; nz1 = subwrd(rec5,9) ; nt1 = subwrd(rec5,12)
  say result
  say '(nx1,ny1,nz1,nt1) = ('nx1','ny1','nz1','nt1')'
  say
  'sdfopen 'wind2_nc_file
  'q file 2'
  rec5 = sublin(result,5)
  nx2 = subwrd(rec5,3) ; ny2 = subwrd(rec5,6) ; nz2 = subwrd(rec5,9) ; nt2 = subwrd(rec5,12)
  say result
  say '(nx2,ny2,nz2,nt2) = ('nx2','ny2','nz2','nt2')'
  say

* make sure (x,y,z,t) dimensions are the same for "wind1_nc_file" and "wind2_nc_file"
* and, if so, define general variable names for (x,y,z,t) dimension limits

  if( nx2 = nx1 )
    nx = nx1
  else
    say
    say 'nx1 = 'nx1
    say 'nx2 = 'nx2
    say '...are not equal'
    say
    exit
  endif
  if( ny2 = ny1 )
    ny = ny1
  else
    say
    say 'ny1 = 'ny1
    say 'ny2 = 'ny2
    say '...are not equal'
    say
    exit
  endif
  if( nz2 = nz1 )
    nz = nz1
  else
    say
    say 'nz1 = 'nz1
    say 'nz2 = 'nz2
    say '...are not equal'
    say
    exit
  endif
  if( nt2 = nt1 )
    nt = nt1
  else
    say
    say 'nt1 = 'nt1
    say 'nt2 = 'nt2
    say '...are not equal'
    say
    exit
  endif

* define (x,y,z,t) lower and upper limits for finding max wind

  x1 = 2 ; x2 = nx-1
  y1 = 2 ; y2 = ny-1
  z1 = 2 ; z2 = nz-1
  t1 = 2 ; t2 = nt-1
  say '(x1,x2) = ('x1','x2')'
  say '(y1,y2) = ('y1','y2')'
  say '(z1,z2) = ('z1','z2')'
  say '(t1,t2) = ('t1','t2')'
  say

* get the x- and y-limits for virtual page

  'q gxinfo'
  rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu = subwrd(rec3,6)
  rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu = subwrd(rec4,6)
  xptime = (xpu-xpl)/2 ; yptime = ypu-0.4

* print plots into an output file?

do_print = true

* plot output

  if ( do_print = true )
    gx_file = 'maxwind2.gx'
    ps_file = 'maxwind2.ps'
    'enable print 'gx_file
  endif

* loop over time

  tt = t1
  while ( tt <= t2 )

* use "wind1_nc_file" to define dimensions

    'set dfile 1'
    'set t 'tt
    'q time' ; ttime = subwrd(result,3)
    say 'tt = 'tt' ('ttime')'

* define wind magnitude variable for this time level

    'set x 1 'nx
    'set y 1 'ny
    'set z 1 'nz
    'define wind = mag('wind1_var'.1,'wind2_var'.2)'

* find max wind and its location

    'set x 1'
    'set y 1'
    'set z 1'
    'set gxout print'
    'd max(max(max(wind,x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
      rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
      max_wind = subwrd(rec,1)
      say ' max wind = 'max_wind
    'd maxloc(max(max(wind,y='y1',y='y2'),z='z1',z='z2'),x='x1',x='x2')'
      rec = sublin(result,(z2-z1+2)*(x2-x1+1)+3)
      xc = subwrd(rec,1)
    'd maxloc(max(max(wind,x='x1',x='x2'),z='z1',z='z2'),y='y1',y='y2')'
      rec = sublin(result,(z2-z1+2)*(y2-y1+1)+3)
      yc = subwrd(rec,1)
    'd maxloc(max(max(wind,x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
      rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
      zc = subwrd(rec,1)
      say ' (xc,yc,zc) = ('xc','yc','zc')'
    'd wind(x='xc',y='yc',z='zc')'
      rec = sublin(result,2)
      wind_xc_yc_zc = subwrd(rec,1)
      say ' wind(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc

* find "world" coordinates of (xc,yc) and convert "world" coordinates
* to "xy" coordinates for plotting track of hurricane center...

      'set x 'xc
      lonc = subwrd(result,4)
      'set y 'yc
      latc = subwrd(result,4)
      'set z 'zc
      levc = subwrd(result,4)
      say ' (lonc,latc,levc) = ('lonc','latc','levc')'
    'd wind(lon='lonc',lat='latc',lev='levc')'
      rec = sublin(result,2)
      wind_lnc_ltc_lvc = subwrd(rec,1)
      say ' wind(lon='lonc',lat='latc',lev='levc') = 'wind_lnc_ltc_lvc

* plot the wind at the level of the max and place a mark at the max

    'c'
    'set vpage 'xpl' 'xpu' 'ypl' 'ypu
    'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
    'set x 'x1' 'x2 ; lon1 = subwrd(result,4) ; lon2 = subwrd(result,5)
    'set y 'y1' 'y2 ; lat1 = subwrd(result,4) ; lat2 = subwrd(result,5)
    'set z 'zc
    'set gxout shaded'
    'set grads off'
    'd 'wind
    'run cbar 0.7'
    'set strsiz 0.16'
    'set string 1 c 6'
    'draw string 'xptime' 'yptime' 'ttime
    'draw title Wind Magnitude ('wind_units') \ max='wind_lnc_ltc_lvc' at: \ (x,y,z)=('xc','yc','zc'), (lon,lat,lev)=('lonc'E,'latc'N,'levc''lev_units')'
    'draw xlab plot: 'lon1'E to 'lon2'E\'
    'draw ylab \ plot: 'lat1'N to 'lat2'N'
    'q w2xy 'lonc' 'latc
      xpos = subwrd(result,3)
      ypos = subwrd(result,6)
      say ' (xpos,ypos) = ('xpos','ypos')'
    'set line 1 1 6'
    'draw mark 6 'xpos' 'ypos' 0.2'
    if ( do_print = true )
      'print'
    endif
    'c events'
    'undefine wind'
***    pause()
    tt=tt+1
  endwhile
  say '...finished.'
  'close 2' ; 'close 1'
  if ( do_print = true )
    say
    'disable print'
    '!gxps -c -i 'gx_file' -o 'ps_file
    '!ps2pdf 'ps_file
    '!rm 'gx_file' 'ps_file
  endif
return

function pause()

* Pauses till ENTER key is pressed 

 say
 say '*********************************'
 say 'Press "ENTER" key to continue...'
 say '*********************************'
 pull answer
return

