function main(args)

ensemble_size = subwrd(args,1)
first_fhr = subwrd(args,2)
fhr = subwrd(args,3)

'set lat 30 50'
'set lon 252.5 280'
'set lev 1000'
*Equally weighted
'define shear06sum1 = 0'
'define shear01sum1 = 0'
'define SRH01sum1 = 0'
'define SRH03sum1 = 0'
'define shear06sumsq = 0'
'define shear01sumsq = 0'
'define SRH01sumsq = 0'
'define SRH03sumsq = 0'
a = 1
total_weight = 0
while (a <= ensemble_size)
'!date'
* shear(0,6000,a,30,50,252.5,280)
 zinterp(ugrdprs,hgtprs,6000,a)
 'u6km = interp'
 zinterp(vgrdprs,hgtprs,6000,a)
 'v6km = interp'
* The 1.94384 factor is to convert from m/s to kts
 'shear06 = 1.94384*mag(u6km-ugrd10m,v6km-vgrd10m)'
 'shear06sum1 = shear06sum1 + shear06'
* 'shear06sum2 = shear06sum2 + ('ensemble_size'-'a')*shear06'
 'shear06sumsq = shear06sumsq + pow(shear06,2)'
* shear(0,1000,a,30,50,252.5,280)
 zinterp(ugrdprs,hgtprs,1000,a)
 'u1km = interp'
 zinterp(vgrdprs,hgtprs,1000,a)
 'v1km = interp'
 'shear01 = 1.94384*mag(u1km-ugrd10m,v1km-vgrd10m)'
 'shear01sum1 = shear01sum1 + shear01'
 'shear01sumsq = shear01sumsq + pow(shear01,2)'
 SRH(0,1000,a,30,50,252.5,280)
 'SRH01sum1 = SRH01sum1 + SRH'
 'SRH01sumsq = SRH01sumsq + pow(SRH,2)'
* SRH(0,3000,a,30,50,252.5,280)
 'SRH03sum1 = SRH03sum1 + HLCY3000_0m.'a
 'SRH03sumsq = SRH03sumsq + pow(HLCY3000_0m.'a',2)'
 total_weight = ensemble_size-a + total_weight
 a = a + 1
endwhile


function SRH(lev_bot,lev_top,member,lat1,lat2,lon1,lon2)

*lev_bot and lev_top are in m
*Dimension environment needs to be set before the function is called

'q ctlinfo'
a = 1
while (a <= 15)
 testline = sublin(result,a)
 testwrd = subwrd(testline,1)
 if (testwrd = 'zdef')
  zmax = subwrd(testline,2)
  break
 endif
 a = a + 1
endwhile
'q dims'
xline = sublin(result,2) ; yline = sublin(result,3)
xmin = subwrd(xline,11) ; xmax = subwrd(xline,13)
ymin = subwrd(yline,11) ; ymax = subwrd(yline,13)
'define SRH = vgrdprs.'member'*0'
'define cu = ave(ugrdprs.'member',lev=1000,lev=500)'
'define cv = ave(vgrdprs.'member',lev=1000,lev=500)'

'set z 1'
'set gxout print'
xx = xmin
yy = ymin
while (xx <= xmax)
 yy = ymin
 n = 1
 while (yy <= ymax)
  'set x 'xx
  'set y 'yy
  zz = 1
  test = 0
  while (zz <= zmax)
   'set z 'zz
   'd hgtprs.'member
   line = sublin(result,2)
   height = subwrd(line,1)
   if (height > lev_bot & test = 0)
    if (zz = 1)
     z0 = 1
    else
     z0 = zz - 1
    endif
    test = 1
   endif       
   if (height > lev_top)
    z1 = zz - 1
    break
   endif
   zz = zz + 1
  endwhile 

  sum.1 = 0; sum.2 = 0; sum.3 = 0; sum.4 = 0
  zz = z0
  i = 1
  while (zz <= z1)
   'set z 'zz
   'd hgtprs.'member
   line = sublin(result,2)
   height.zz = subwrd(line,1)
* 4 terms - term 1 ****************
   if (z0 = 1)
    'd -( (vgrdprs.'member'(z+1) - vgrdprs.'member') / (hgtprs.'member'(z+1)-hgtprs.'member') )*ugrdprs.'member
   else
    'd -( (vgrdprs.'member'(z+1) - vgrdprs.'member'(z-1)) / (hgtprs.'member'(z+1)-hgtprs.'member'(z-1)) )*ugrdprs.'member
   endif              
   line = sublin(result,2)
   data.1.i = subwrd(line,1)
   if (i > 1)           
    ii = i - 1 ; zzz = zz - 1
    sum.1 = sum.1 + (data.1.i + data.1.ii)/2*(height.zz - height.zzz)
   endif
* Term 2 *******************************
   if (z0 = 1)
    'd ((ugrdprs.'member'(z+1)-ugrdprs.'member')/(hgtprs.'member'(z+1)-hgtprs.'member'))*vgrdprs.'member
   else
    'd ((ugrdprs.'member'(z+1)-ugrdprs.'member'(z-1))/(hgtprs.'member'(z+1)-hgtprs.'member'(z-1)))*vgrdprs.'member
   endif
   line = sublin(result,2)
   data.2.i = subwrd(line,1)
   if (i > 1)
    ii = i - 1 ; zzz = zz - 1
    sum.2 = sum.2 + (data.2.i + data.2.ii)/2*(height.zz - height.zzz)
   endif
* Term 3  *******************************
   'd cu'
   line = sublin(result,2)
   cx = subwrd(line,1)
   if (z0 = 1)
    'd ((vgrdprs.'member'(z+1)-vgrdprs.'member')/(hgtprs.'member'(z+1)-hgtprs.'member'))'
   else
    'd ((vgrdprs.'member'(z+1)-vgrdprs.'member'(z-1))/(hgtprs.'member'(z+1)-hgtprs.'member'(z-1)))'
   endif
   line = sublin(result,2)
   data.3.i = subwrd(line,1)
   if (i > 1)
    ii = i - 1 ; zzz = zz - 1
    sum.3 = sum.3 + cx*(data.3.i + data.3.ii)/2*(height.zz - height.zzz)
   endif       
* Term 4 **********************************
   'd cv'      
   line = sublin(result,2)
   cy = subwrd(line,1)
   if (z0 = 1)
    'd ((ugrdprs.'member'(z+1)-ugrdprs.'member')/(hgtprs.'member'(z+1)-hgtprs.'member'))'
   else
    'd ((ugrdprs.'member'(z+1)-ugrdprs.'member'(z-1))/(hgtprs.'member'(z+1)-hgtprs.'member'(z-1)))'
   endif
   line = sublin(result,2)
   data.4.i = subwrd(line,1)
   if (i > 1)
    ii = i - 1 ; zzz = zz - 1
    sum.4 = sum.4 - cy*(data.4.i + data.4.ii)/2*(height.zz - height.zzz)
   endif
*************************
   zz = zz + 1 ; i = i + 1
  endwhile 
* Sum
  helicity = sum.1 + sum.2 + sum.3 + sum.4  
  'set defval SRH 'xx' 'yy' 'helicity
  yy = yy + 1
 endwhile                  
 xx = xx + 1
endwhile

'set lat 'lat1' 'lat2
'set lon 'lon1' 'lon2

return

function zinterp(field,zgrid,zlev,member)
zlev = "(hgtsfc."member"(z=1)+"zlev")"
"q dims"
rec=sublin(result,4)
ztype=subwrd(rec,3)
if (ztype = "fixed")
   zmin=subwrd(rec,9)
   zmax=zmin
else
   zmin=subwrd(rec,11)
   zmax=subwrd(rec,13)
endif

* Get full z-dimensions of dataset.

"q file"
rec=sublin(result,5)
zsize=subwrd(rec,9)

* Determine spatially varying bounding height levels for height surface
* zabove = height-value at level above ; zbelow = height value at level
* below for each gridpoint

"set z 2 "zsize
"define zabove=0.5*maskout("zgrid"."member","zgrid"."member"-"zlev")+0.5*maskout("zgrid"."member","zlev"-"zgrid"."member"(z-1))"
"set z 1 "zsize-1
"define zbelow=0.5*maskout("zgrid"."member","zgrid"."member"(z+1)-"zlev")+0.5*maskout("zgrid"."member","zlev"-"zgrid"."member")"

* Isolate field values at bounding height levels
* fabove = requested field value above height surface
* fbelow = requested field value below height surface

"set z 2 "zsize
"define fabove=zabove*0+"field"."member
"set z 1 "zsize-1
"define fbelow=zbelow*0+"field"."member

* Turn this 3-D grid of values (mostly undefined) into a 2-D height layer
* mean is used here below for simplicity, since mean ignores undefined
* values.

"set z 1"
"define zabove=mean(zabove,z=2,z="zsize")"
"define fabove=mean(fabove,z=2,z="zsize")"
"define zbelow=mean(zbelow,z=1,z="zsize-1")"
"define fbelow=mean(fbelow,z=1,z="zsize-1")"

* Finally, interpolate linearly in height and create surface.

"set z "zmin " " zmax

"define slope=(fabove-fbelow)/(zabove-zbelow)"
"define b=fbelow-slope*zbelow"
"define interp=slope*"zlev"+b"

return