*This will be an M x N multi panel image. To make it the dimensions you want, set M_ and N_.  The other odd subtractions are just due to how 
*much empty space you want around the outsides of panel (towards the edge of the page) and space between panels.

'reinit'
'set display color white' ; 'c'
'set font 1'

*Set dimensions of the plotting matrix
**********************************************
*M_ = number of rows (vertical dimension)
M_ = 2
*N_ = number of columns (horizontal dimension)
N_ = 2
**********************************************
*Set page size and/or aspect ratio
************************************
x_page_dim = 11
y_page_dim = 8.5
aspect_ratio = x_page_dim/y_page_dim
xpixels = 1280
ypixels = xpixels/aspect_ratio
************************************
*Set spacing parameters (in real page inches)
***************************************************************
l_margin = 0.5
r_margin = 0.2
b_margin = 0.5
*Top margin includes space for a title
t_margin = 0.1
h_pad = 0.6
v_pad = 0.33
*These two serve as markers for where to start the virtual pages (so where the lower left corner of the images will be)
*These two control the actual width of each panel
panel_height = (y_page_dim-t_margin-b_margin-(M_-1)*v_pad)/M_
panel_width = (x_page_dim-l_margin-r_margin-(N_-1)*h_pad)/N_
title_y = (y_page_dim-0.5*t_margin)
****************************************************************

* Set removed_ens to anything other than 1, 2, or 3 to eliminate plotting sub-ensemble spread altogether
removed_ens = 3
thermo = 1
'open /work/jdduda/analysis/FLSM/WRFPRS_FLSM.ctl'
'open /work/jdduda/analysis/MLSM/WRFPRS_MLSM.ctl'
'open /work/jdduda/analysis/LSMO/WRFPRS_LSMO.ctl'
ens_size = 10
count_size = 9
if (0)
a = 1
while (a <= ens_size)
 'open /work/jdduda/analysis/FLSM/m'a'/wrfout_scalar.ctl'
 a = a + 1
endwhile
a = 1
while (a <= ens_size)
 'open /work/jdduda/analysis/MLSM/m'a'/wrfout_scalar.ctl'
 a = a + 1
endwhile
a = 1
while (a <= ens_size)
 'open /work/jdduda/analysis/LSMO/m'a'/wrfout_scalar.ctl'
 a = a + 1
endwhile
endif

pres_lev = 700
name.1 = "FLSM"
name.2 = "MLSM"
name.3 = "LSMO"

'set grads off'
'set tlsupp year'
'q dims'
xline = sublin(result,2) ; yline = sublin(result,3)
lon1 = subwrd(xline,6) ; lon2 = subwrd(xline,8)
lat1 = subwrd(yline,6) ; lat2 = subwrd(yline,8)
'set xlopts 1 1 0.09' ; 'set ylopts 1 1 0.09'
'set x 1' ; 'set y 1' ; 'set e 1' ; 'set t 2'
*'define gridsize = asum(const(tmp2m,1),lon='lon1',lon='lon2',lat='lat1',lat='lat2')'
'set t 2 37'

color.1 = 3
color.2 = 2
*If using HRLDAS, use color 5
*color.3 = 5
*If using LSMO, use color 9
color.3 = 9

line.1 = 1
line.2 = 2
line.3 = 1
line.4 = 6
line.5 = 1
line.6 = 1
line.7 = 3
line.8 = 1
line.9 = 1
line.10 = 1
mark.1 = 0
mark.2 = 0
mark.3 = 5
mark.4 = 0
mark.5 = 1
mark.6 = 6
mark.7 = 0
mark.8 = 9
mark.9 = 12
mark.10 = 3

a = M_-1
b = 0
while (a >= 0)
*a is the row
 while (b < N_)
*b is the column
  if (b = 0)
   xi.a.b = l_margin + b*panel_width 
  else
   xi.a.b = l_margin + b*(panel_width + h_pad)
  endif
  xf.a.b = xi.a.b + panel_width
  if (a = M_-1)
   yf.a.b = (y_page_dim-t_margin) - (M_-a-1)*panel_height
  else
   yf.a.b = (y_page_dim-t_margin) - (M_-a-1)*(panel_height + v_pad)
  endif
  yi.a.b = yf.a.b - panel_height
  b = b + 1
 endwhile
 b = 0
 a = a - 1
endwhile

a = M_-1
b = 0
n = 1
while (a >= 0)
*a is the row
 while (b < N_)
*b is the column
  'set parea 'xi.a.b' 'xf.a.b' 'yi.a.b' 'yf.a.b

  x_pos = xi.a.b - 0.4
  y_pos = 0.5*(yi.a.b+yf.a.b)
  'set string 1 c 6 90'
  'set strsiz 0.125'

  'set cthick 6'

  if (n = 1)
   if (thermo)
    'set vrange 0 500'
   else
    'set lev 'pres_lev
    'set vrange 0 2.25'
   endif
   ens = 1
   while (ens <= 3)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*vgrdprs'
    'define sumsq = 0*vgrdprs'
    ee = 1
    while (ee <= ens_size)
     'set e 'ee
     if (thermo)
      say "fuck "ee
      'total = total + capesfc.'ens
      'sumsq = sumsq + pow(capesfc.'ens',2)'
     else
      'total = total + vgrdprs.'ens
      'sumsq = sumsq + pow(vgrdprs.'ens',2)'
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('ens_size'-1))*(sumsq-(1/'ens_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    if (ens = removed_ens) ; 'set ccolor 1' ; endif
    'set cmark 0'
    'set x 1' ; 'set y 1'
    'd sqrt(avevar)'

* plot spread with each member taken out
    if (ens = removed_ens)
    e3 = 1
    while (e3 <= ens_size)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*vgrdprs'
    'define sumsq = 0*vgrdprs'
    ee = 1 
    while (ee <= ens_size)
     if (ee != e3)
     'set e 'ee
     if (thermo)
      'total = total + capesfc.'ens
      'sumsq = sumsq + pow(capesfc.'ens',2)'
     else   
      'total = total + vgrdprs.'ens
      'sumsq = sumsq + pow(vgrdprs.'ens',2)'
     endif
     endif            
     ee = ee + 1
    endwhile
    'define var1 = (1/('count_size'-1))*(sumsq-(1/'count_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    'set cthick 1' ; 'set ccolor 'color.ens ; 'set cstyle 'line.e3 ; 'set digsize 0.05' ;  'set cmark 'mark.e3    
    'set x 1' ; 'set y 1'
    'd sqrt(avevar)'
    e3 = e3 + 1
    endwhile
    endif
    ens = ens + 1
   endwhile
   if (thermo)
    'draw string 'x_pos' 'y_pos' SBCAPE (J/kg)'
   else
    'draw string 'x_pos' 'y_pos' 'pres_lev'-mb v-wind (m/s)'
   endif
  endif
  if (n = 2)
   if (thermo)
    'set lev 0 50'
   else
    'set lev 'pres_lev
    'set vrange 0 2.25'
   endif
   ens = 1
   while (ens <= 3)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*ugrdprs'
    'define sumsq = 0*ugrdprs'
    ee = 1
    while (ee <= ens_size)
     'set e 'ee
     if (thermo)
      'total = total + cinsfc.'ens
      'sumsq = sumsq + pow(cinsfc.'ens',2)'
     else
      'total = total + ugrdprs.'ens
      'sumsq = sumsq + pow(ugrdprs.'ens',2)'
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('ens_size'-1))*(sumsq-(1/'ens_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    if (ens = removed_ens) ; 'set ccolor 1' ; endif
    'set cmark 0'
    'set x 1' ; 'set y 1'
    'd sqrt(avevar)'

* plot spread with each member taken out
    if (ens = removed_ens)
    e3 = 1
    while (e3 <= ens_size)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*ugrdprs'
    'define sumsq = 0*ugrdprs'
    ee = 1
    while (ee <= ens_size)
     'set e 'ee
     if (ee != e3)
     if (thermo)
      'total = total + cinsfc.'ens
      'sumsq = sumsq + pow(cinsfc.'ens',2)'
     else
      'total = total + ugrdprs.'ens
      'sumsq = sumsq + pow(ugrdprs.'ens',2)'
     endif
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('count_size'-1))*(sumsq-(1/'count_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    'set cthick 1' ; 'set ccolor 'color.ens ; 'set cstyle 'line.e3 ; 'set digsize 0.05' ;  'set cmark 'mark.e3    
    'set x 1' ; 'set y 1'
    'd sqrt(avevar)'
    e3 = e3 + 1
    endwhile
    endif
    ens = ens + 1
   endwhile
   if (thermo)
    'draw string 'x_pos' 'y_pos' SBCIN (J/kg)'
   else
    'draw string 'x_pos' 'y_pos' 'pres_lev'-mb u-wind (m/s)'
   endif
  endif
  if (n = 3)
   if (thermo)
    'set vrange 0 400'
   else
    'set lev 'pres_lev
    'set vrange 0 1.25'
   endif
   ens = 1
   while (ens <= 3)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*tmpprs'
    'define sumsq = 0*tmpprs'
    ee = 1
    while (ee <= ens_size)
     'set e 'ee
     if (thermo)
      'total = total + cape90_0mb.'ens
      'sumsq = sumsq + pow(cape90_0mb.'ens',2)'
     else
      'total = total + tmpprs.'ens
      'sumsq = sumsq + pow(tmpprs.'ens',2)'
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('ens_size'-1))*(sumsq-(1/'ens_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    if (ens = removed_ens) ; 'set ccolor 1' ; endif
    'set cmark 0'
    'set x 1' ; 'set y 1'
    'd sqrt(avevar)'

* plot spread with each member taken out
    if (ens = removed_ens)
    e3 = 1
    while (e3 <= ens_size)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*tmpprs'
    'define sumsq = 0*tmpprs'
    ee = 1
    while (ee <= ens_size)
     if (ee != e3)
     'set e 'ee
     if (thermo)
      'total = total + cape90_0mb.'ens
      'sumsq = sumsq + pow(cape90_0mb.'ens',2)'
     else
      'total = total + tmpprs.'ens
      'sumsq = sumsq + pow(tmpprs.'ens',2)'
     endif
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('count_size'-1))*(sumsq-(1/'count_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    'set cthick 1' ; 'set ccolor 'color.ens ; 'set cstyle 'line.e3 ; 'set digsize 0.05' ;  'set cmark 'mark.e3    
    'set x 1' ; 'set y 1'
    'd sqrt(avevar)'
    e3 = e3 + 1
    endwhile
    endif

    ens = ens + 1
   endwhile
   if (thermo)
    'draw string 'x_pos' 'y_pos' MLCAPE (J/kg)'
   else
    'draw string 'x_pos' 'y_pos' 'pres_lev'-mb `3z`1 (K)'
   endif
  endif
  if (n = 4)
   if (thermo)
    'set vrange 0 50'
   else
    'set lev 'pres_lev
    'set vrange 0 1.25'
   endif
   ens = 1
   while (ens <= 3)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    if (!thermo)
     'set e 1 10'
     'define q = spfhprs.'ens'/(1-spfhprs.'ens')'
     'set e 1'
    endif
    'define total = 0*ugrdprs'
    'define sumsq = 0*ugrdprs'
    ee = 1
    while (ee <= ens_size)
     'set e 'ee
     fn = 10*(ens-1)+2+ee
     if (thermo)
      'total = total + cin90_0mb.'ens
      'sumsq = sumsq + pow(cin90_0mb.'ens',2)'
     else
      'total = total + q'
      'sumsq = sumsq + pow(q,2)'
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('ens_size'-1))*(sumsq-(1/'ens_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    if (ens = removed_ens) ; 'set ccolor 1' ; endif
    'set cmark 0'
    'set x 1' ; 'set y 1'
    'd 1000*sqrt(avevar)'

* plot spread with each member taken out
    if (ens = removed_ens)
    e3 = 1
    while (e3 <= ens_size)
    'set lon 'lon1' 'lon2
    'set lat 'lat1' 'lat2
    'define total = 0*vgrdprs'
    'define sumsq = 0*vgrdprs'
    ee = 1
    while (ee <= ens_size)
     if (ee != e3)
     'set e 'ee
     if (thermo)
      'total = total + cin90_0mb.'ens
      'sumsq = sumsq + pow(cin90_0mb.'ens',2)'
     else
      'define q = spfhprs.'ens'/(1-spfhprs.'ens')'
      'total = total + q'
      'sumsq = sumsq + pow(q,2)'
     endif
     endif
     ee = ee + 1
    endwhile
    'define var1 = (1/('count_size'-1))*(sumsq-(1/'count_size')*pow(total,2))'
    'define avevar = tloop(amean(var1,lon='lon1',lon='lon2',lat='lat1',lat='lat2'))'
    'set ccolor 'color.ens
    'set cthick 1' ; 'set ccolor 'color.ens ; 'set cstyle 'line.e3 ; 'set digsize 0.05' ;  'set cmark 'mark.e3
    'set x 1' ; 'set y 1'
    'd 1000*sqrt(avevar)'
    e3 = e3 + 1
    endwhile
    endif

    ens = ens + 1
   endwhile
   if (thermo)
    'draw string 'x_pos' 'y_pos' MLCIN (J/kg)'
   else
    'draw string 'x_pos' 'y_pos' 'pres_lev'-mb mixing ratio (g/kg)'
   endif
  endif

*Drawing the boxes to indicate panel title
   if (n = 1)
   x1 = xf.a.b-0.6
   x2 = xf.a.b
   y1 = yi.a.b
   y2 = yi.a.b+0.75
   'set line 0'
   'draw recf 'x1' 'y1' 'x2' 'y2
   'set line 1 1 6'
   'draw rec 'x1' 'y1' 'x2' 'y2
   'set strsiz 0.08'
   xpos = 0.5*(x1+x2)
   dy = (y2-y1)/3
   dystag = dy/2
   nn = 1
   while (nn <= 3)
    'set string 'color.nn' c 3 0'
    ypos = y2-dy*(nn-1)-dystag
    'draw string 'xpos' 'ypos' 'name.nn
    nn = nn + 1         
   endwhile
   endif

* for members removed
   if (n = 2)
   x1 = xi.a.b
   x2 = xi.a.b+1.4
   y1 = yf.a.b-1.75
   y2 = yf.a.b
   'set line 0'
   'draw recf 'x1' 'y1' 'x2' 'y2
   'set line 1 1 6'      
   'draw rec 'x1' 'y1' 'x2' 'y2
   'set strsiz 0.08'      
   xpos1 = x1 + 0.1*(x2-x1)
   xpos2 = x1 + 0.8*(2/3)*(x2-x1)
   xmid = 0.5*(xpos1+xpos2)
   xpos3 = xpos2 + 0.05
   dy = (y2-y1)/10
   dystag = dy/2         
   nn = 1         
   'set string 1 l 3 0'            
   while (nn <= 10)
    'set line 'color.3' 'line.nn' 1'
    ypos = y2-dy*(nn-1)-dystag
    'draw line 'xpos1' 'ypos' 'xpos2' 'ypos
    'draw mark 'mark.nn' 'xmid' 'ypos' 0.1'
    'draw string 'xpos3' 'ypos' wo/m'nn
    nn = nn + 1
   endwhile       
   endif

  b = b + 1
  n = n + 1
 endwhile
 b = 0
 a = a - 1
endwhile

'printim LSM_thermo_spread_wo_LSMO.png x'xpixels' y'ypixels