* To obtain (i) the largest of the maximum wind magnitudes from 4 standard pressure levels
* (related to the low level jet) over a defined area, 
* (ii) its corresponding pressure level and (iii) the grid coordinates.
*
* Done by Ooi See Hai on 26 March 2010
*
  'reinit'
*
  'sdfopen d:\cdc\uvw\uwnd.2008.nc'
  'sdfopen d:\cdc\uvw\vwnd.2008.nc'
*
*
  l.1=1000
  l.2=925
  l.3=850
  l.4=700
*
  'set t 124'
  'q time'
  act=subwrd(result,3)
*
* Array header
* ............
*
  print
  print ' Pressure      Magnitude'
  print 'level(hPa) maximum wind(m/s) for defined area: 20S-40N & 30-100E '
  print '.......... ................. at 'act
  print
*
*
  n=1
  while (n <= 2)
*
* Array of maximum wind over a specified area for the respective pressure levels
* ..............................................................................
*
  i=1
  while (i <= 4)
    'set lev 'l.i
    'wind=mag(uwnd,vwnd.2)'
    'd max(max(wind,lon=30,lon=100),lat=-20,lat=40)'
    MAXING=sublin(result,27)
    val=subwrd(MAXING,4)
    mw.i=val
*
* Display array for later verification
*
    if (n=1)
       print '  'l.i'          'mw.i
    endif
    i=i+1 
  endwhile
*
*
*
    if (n = 1)
*
* Obtain largest of the maximum magnitude by sorting
* ...................................................
*
      j=1
      while (j <= 4)
        id=j
        k=j+1
        while (k <= 4)
          if (mw.k > mw.id)
            id=k
          endif
          k=k+1
        endwhile
*
        if(k != id)
          large= mw.k
          mw.k=mw.id
          mw.id=large
        endif
*
        j=j+1
      endwhile
*
      print
      print '  largest maximum wind magnitude = 'mw.id' m/s'
      print
*
    else
*
* the index of the corresponding pressure level (from unsorted data)
* ..................................................................
*
      p=1
      while (p <= 4)
        if (mw.p = large)
          pli=p
          print '  corresponding pressure level : 'l.pli' hPa'
        endif
        p=p+1
      endwhile
*
*  grid coordinates
*  ................
*
      'set lev 'l.pli
*
      'd maxloc(max(wind,lon=30,lon=100),lat=-20,lat=40)'
      line=sublin(result,27)
*
* Note: There are 25 grid points (or MAXing lines) between the latitude band.
*       Add 1 MAXLOCing and 1 Result value lines.  Total = 27
*
      ygrd=subwrd(line,4)
*
      'd maxloc(max(wind,lat=-20,lat=40),lon=30,lon=100)'
      line=sublin(result,31)
      xgrd=subwrd(line,4)
*
      'set x 'xgrd
      lonval=subwrd(result,4)
      'set y 'ygrd
      latval=subwrd(result,4)
*
      if (latval < 0)
        nlatval=-latval
      else
        nlatval=latval
      endif
*
      print
      if (latval < 0)
        print '  at grid points : 'nlatval' deg S and 'lonval' deg E '
      else
        print '  at grid points : 'nlatval' deg N and 'lonval' deg E '
      endif
      print
*
*
     endif
*
     n=n+1
   endwhile
*
*