* To obtain (i) the largest of the maximum wind magnitudes from 12 standard pressure levels
* (related to the subtropical jet stream) over a defined area, 
* (ii) its corresponding pressure level and (iii) the grid coordinates.
*
* Done by Kumar Subramaniam (Malaysian Meteorological Department)
* and Ooi See Hai (a retired meteorologist) on 25 March 2010
*
  'reinit'
*
  'sdfopen d:\cdc\uvw\uwnd.2009.nc'
  'sdfopen d:\cdc\uvw\vwnd.2009.nc'
*
  l.1=1000
  l.2=925
  l.3=850
  l.4=700
  l.5=600
  l.6=500
  l.7=400
  l.8=300
  l.9=250
  l.10=200
  l.11=150
  l.12=100
*
  'set t 60'
*
*
  n=1
  while (n <= 2)
*
* Array of maximum wind over a specified area for the respective pressure levels
* ..............................................................................
*
  i=1
  while (i <= 12)
    'set lev 'l.i
    'wind=mag(uwnd,vwnd.2)'
    'd max(max(wind,lon=80,lon=140),lat=20,lat=40)'
    MAXING=sublin(result,11)
    val=subwrd(MAXING,4)
    mw.i=val    
    i=i+1 
  endwhile
*
*
*
    if (n = 1)
*
* Obtain largest of the maximum magnitude 
* ...................................................
* (Sorting using selection method in ascending order)
*
      j=1
      while (j <= 12)
        id=j
        k=j+1
        while (k <= 12)
          if (mw.k < mw.id)
            id=k
          endif
          k=k+1
        endwhile
*
        if(j != id)
          large= mw.j
          mw.j=mw.id
          mw.id=large
        endif
*
        j=j+1
      endwhile
*
      print
      print
      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 <= 12)
        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=80,lon=140),lat=20,lat=40)'
      line=sublin(result,11)
*
* Note: There are 9 grid points (or MAXing lines) between the latitude band.
*       Add 1 MAXLOCing and 1 Result value lines.  Total = 11
*
      ygrd=subwrd(line,4)
*
      'd maxloc(max(wind,lat=20,lat=40),lon=80,lon=140)'
      line=sublin(result,27)
      xgrd=subwrd(line,4)
*
      'set x 'xgrd
      lonval=subwrd(result,4)
      'set y 'ygrd
      latval=subwrd(result,4)
*
      print
      print '  at grid points : 'latval' deg N and 'lonval' deg E '
      print
      print
*
*
     endif
*
     n=n+1
   endwhile
*
*