require, "nc.i"
require, "cgd.i"
require, "plclabs_bab.i"
require, "style.i"

func plez (data, levs=, nlevs=)
/* DOCUMENT plez, data

   Makes a contour plot of a 2d array, the mesh is set from the
   indices of the input array.

*/
{
  /* check dimensions of data */
  dims = dimsof(data);
  ndims= dims(1);
  if (ndims != 2) {
    print, "whats this? dims:", dims;
    return;
  }

  /* Set contours and min/max values for colour fill */
  if (is_void(nlevs)) nlevs=10;
  if (is_void(levs))
    levs = set_contours (data, nlevs = nlevs);

  nx = dims(2);
  ny = dims(3);
  x = indgen(nx);
  y = indgen(ny);
  xc = x(*)(,-:1:ny);
  yc = y(*)(-:1:nx,);

  plmesh, yc, xc;
  plc,        data, yc, xc, levs=levs, marks=0;
  plc_labels, data, yc, xc, levs=levs, opaque=1;
}

func svp (vpnum)

/* DOCUMENT svp, viewport_number
   
   Set the current viewport (with plsys) and make a global variable
   (viewport_number) so that other functions can perform operations
   in system 0 and restore the current system when done
   */
{ 
  extern viewport_number;
  extern vport;

  plsys, vpnum;
  viewport_number = vpnum;
  vport = viewport();
}


func pln (nc_var, 
          cint=, levs=, colors=, nlab=, nlevs=, width=, concolor=,
	  neg=, type=, zlinear=, mask=, outline=, con_lab_height=,  
	  opaque=, cbar=, mxmn=, 
	  nofill=, nocont=, nocbar=, nocornlab=, noname=, flip=, setpal=,
	  proj=, xgrid=, ygrid=, anim=, x=, y=, z=, t= )
/* DOCUMENT pln, nc_var
   
   Makes a contour plot from a NC_variable structure, normally obtained by
   nc_get_var

   KEYWORDS: nlevs, number of contour levels
             cint, contour interval to use
             levs, actual contour levels, takes precedence over nlevs or cint
	     colors, list of color indices to use (numberof(levs)+1)
	     nlab, label every nlab contours
	        2 values: first contour to label and increment
	     width, 1 value: line width for contours, or,
	        a 2 element vector giving labelled and unlabelled widths
	     lwidth, line width for labelled contours
	     concolor, color to draw contours in ("red", black", etc)
	     mxmn != 0, label extrema with their values, 
                   > 1, label mxmn most extreme max's and mins
		   < 0, label abs(mxmn) extrema spaced through array
	     neg = style for negative contours (type= in plc call)
             type, one of "xy", "xz", "xt", "yz", "yt", "zt", the type
	           of 2d plot to make; of "x", "y", "z", "t", the type
		   of 1D plot to make.
	     flip != 0 transpose axes of plot
	     zlinear, nonzero: z axis in linear instead of logarthmic
	     mask, 2D array containing a region mask (such as an "ORO"
	           field), whose dimensions must match those dimensions of
		   *nc_var.data which are to be plotted.
		   Points where mask=0 are not plotted.
	     outline = 2d array, containing a region to outline. Often,
	           outline and mask will be the same array, to mask a region
		   and put an outline arround it (e.g. land or ocean)
	     con_lab_height = height (size) for contour labels
	     opaque = if nonnil, make contour labels opaque
	     cbar,  "adjust" argument to color_bar, 
	             moves color bar closer (-) or farther (+) from plot
		     default (void) gives cbar = .025
	     nofill, non-zero: no colour fill
	     nocbar, >0: no color bar, even if nofill=0
	             <0: horizontal color bar instead of vertical
	     setpal, nonzero: generate a new palette if nofill=0
	             >1, extend the current palette instead of replacing it.
		     < 0 reverse the order of the palette colors and use abs(setpal)
	     nocont, non-zero: no contours plotted
	     nocornlab, non-zero: corner label not plotted
	     noname, non-zero: name and unit label plotted
             proj, map projection (mollweide, polar, orthographic)
	           for polar and orthographic projections the external variables: 
		       polat and polon govern the projection center and rotation;
		       cirlat determines the outer latitude (+/-20 is a good choice), 
		       the circle is drawn and points outside are masked.
	     xgrid, ygrid=interval: add a grid a the specified interval (maps only)
	     anim, non-nil: the time dimension is animated, repeated anim times
	     x, y, z, t= single coordinate value to plot or coordinate range to average
*/
{
  local iregc, irego;
  
  /* make sure that the input variable is a proper structure */
  if (nameof(structof(nc_var)) != "NC_variable") {
    print, "input variable is not a structure of type NC_variable";
    info, nc_var;
    return;
  }

  /* check dimensions of data */
  dims = dimsof(*nc_var.data);
  ndims= dims(1);
  dims = dims(2:);
  if (ndims != 4) {
    print, "whats this? dims:", dims;
    return;
  }

  /* check inputs, set voids to 0 */
  if (is_void(zlinear)) zlinear = 0;
  if (is_void(nofill)) nofill = 0;
  if (is_void(nocbar)) nocbar = 0;
  if (is_void(setpal)) setpal = 0;
  if (is_void(nocont)) nocont = 0;
  if (is_void(noname)) noname = 0;
  if (is_void(nocornlab)) nocornlab = 0;
  if (is_void(cbar)) cbar = 0.025;
  if (is_void(mxmn)) mxmn = 0;
  if (anim==0) anim=[];

  // set the missing value for fillc
  local spv;
  if (!is_void(*nc_var.spv)) spv = (*nc_var.spv)(1);
  
  /* Determine default type, if necessary */
  if (is_void(type)) {
    if (dims(3) > 1) {
      if (dims(2) > 1) {
	type = "yz";
      } else if (dims(1) > 1) {
	type = "xz";
      } else if (dims(4) > 1) {
	type = "zt";
      } else {
	type = "z";
      }
    } else if (dims(2) > 1) {
      if (dims(1) > 1) {
	type = "xy";
      } else if (dims(4) > 1) {
	type = "yt";
      } else {
	type = "y";
      }
    } else if (dims(1) > 1) {
      if (dims(4) > 1) {
	type = "xt";
      } else {
	type = "x";
      }
    } else {
      type = "z";
    }
  }

  /* Set up the coordinate arrays from which the plot mesh is
     constructed and reduce the data to 2 dimensions. */

  /* xy section */
  if (type == "xy") {
    xl = *nc_var.x;
    yl = *nc_var.y;

    /* z averaging values */
    zlist = make_sublist(*nc_var.z, "z", dims(3), rng=z)
    if (numberof(zlist)==0) return;
    zwt = (*nc_var.zw)(zlist) / ((*nc_var.zw)(zlist)(sum));

    /* t averaging values */
    tlist = make_sublist(*nc_var.t, "t", dims(4), rng=t)
    if (numberof(tlist)==0) return;

    /* average data */
    if (is_void(anim)) {
      data = (*nc_var.data)(,,zlist,tlist)(,,,avg)(,,+) * zwt(+);
    } else {
      data = (*nc_var.data)(,,zlist,tlist)(,,+,) * zwt(+);
    }

  /* xz section */
  } else if (type == "xz") {
    xl = *nc_var.x;
    yl = *nc_var.z;
    if (!zlinear) yl = -7*log(yl/1000.);

    /* y averaging values */
    ylist = make_sublist(*nc_var.y, "y", dims(2), rng=y)
    if (numberof(ylist)==0) return;
    ywt = (*nc_var.yw)(ylist) / ((*nc_var.yw)(ylist)(sum));

    /* t averaging values */
    tlist = make_sublist(*nc_var.t, "t", dims(4), rng=t)
    if (numberof(tlist)==0) return;

    /* average data */
    if (is_void(anim)) {
      data = (*nc_var.data)(,ylist,,tlist)(,,,avg)(,+,) * ywt(+);
    } else {
      data = (*nc_var.data)(,ylist,,tlist)(,+,,) * ywt(+);
    }
  
  /* xt section */
  } else if (type == "xt") {
    xl = *nc_var.x;
    yl = *nc_var.t;

    /* y averaging values */
    ylist = make_sublist(*nc_var.y, "y", dims(2), rng=y)
    if (numberof(ylist)==0) return;
    ywt = (*nc_var.yw)(ylist) / ((*nc_var.yw)(ylist)(sum));

    /* z averaging values */
    zlist = make_sublist(*nc_var.z, "z", dims(3), rng=z)
    if (numberof(zlist)==0) return;
    zwt = (*nc_var.zw)(zlist) / ((*nc_var.zw)(zlist)(sum));

    /* average data */
    data = ((*nc_var.data)(,,zlist,)(,,+,) * zwt(+))(,ylist,)(,+,) * ywt(+);

  /* yz section */
  } else if (type == "yz") {
    xl = *nc_var.y;
    yl = *nc_var.z;
    if (!zlinear) yl = -7*log(yl/1000.);

    /* x averaging values */
    xlist = make_sublist(*nc_var.x, "x", dims(1), rng=x);
    if (numberof(xlist)==0) return;
    xwt = (*nc_var.xw)(xlist) / ((*nc_var.xw)(xlist)(sum));

    /* t averaging values */
    tlist = make_sublist(*nc_var.t, "t", dims(4), rng=t)
    if (numberof(tlist)==0) return;

    /* average data */
    if (is_void(anim)) {
      data = (*nc_var.data)(xlist,,,tlist)(,,,avg)(+,,) * xwt(+);
    } else {
      data = (*nc_var.data)(xlist,,,tlist)(+,,,) * xwt(+);
    }

  /* yt section */
  } else if (type == "yt") {
    xl = *nc_var.y;
    yl = *nc_var.t;

    /* x averaging values */
    xlist = make_sublist(*nc_var.x, "x", dims(1), rng=x);
    if (numberof(xlist)==0) return;
    xwt = (*nc_var.xw)(xlist) / ((*nc_var.xw)(xlist)(sum));

    /* z averaging values */
    zlist = make_sublist(*nc_var.z, "z", dims(3), rng=z)
    if (numberof(zlist)==0) return;
    zwt = (*nc_var.zw)(zlist) / ((*nc_var.zw)(zlist)(sum));

    /* average data */
    data = ((*nc_var.data)(,,zlist,)(,,+,)  *zwt(+))(xlist,,)(+,,) * xwt(+);

  /* zt section */
  } else if (type == "zt") {
    xl = *nc_var.t;
    yl = *nc_var.z;
    if (!zlinear) yl = -7*log(yl/1000);

    /* x averaging values */
    xlist = make_sublist(*nc_var.x, "x", dims(1), rng=x);
    if (numberof(xlist)==0) return;
    xwt = (*nc_var.xw)(xlist) / ((*nc_var.xw)(xlist)(sum));

    /* y averaging values */
    ylist = make_sublist(*nc_var.y, "y", dims(2), rng=y)
    if (numberof(ylist)==0) return;
    ywt = (*nc_var.yw)(ylist) / ((*nc_var.yw)(ylist)(sum));

    /* average and transpose data */
    data = transpose(((*nc_var.data)(,ylist,,)(,+,,) * ywt(+))(xlist,,)(+,,) * xwt(+));

  } else if (type == "x") {
    print, "type x not ready yet";
    return;
  } else if (type == "y") {
    print, "type y not ready yet";
    return;
  } else if (type == "z") {
    print, "type z not ready yet";
    return;
  } else if (type == "t") {
    print, "type t not ready yet";
    return;
  } else {
    print, "unknown type: ", type;
    return;
  }

  /* transpose 1st and 2nd dimensions if requested */
  if (!is_void(flip) && (flip != 0)) {
    tmp = xl;
    xl = yl;
    yl = tmp;
    data = transpose(data,[2,1]);
  }

  /* check that the reduced data has dimensions > 1 */
  nx = numberof(xl);
  ny = numberof(yl);
  if (nx <= 1) {
    print, type, " first dimension of plot is 1, array dimensions:", dims;
    return;
  }
  if (ny <= 1) {
    print, type, " second dimension of plot is 1, array dimensions:", dims;
    return;
  }

  /* set up preliminary transformation arrays */
  xc = xl(*)(,-:1:ny);
  yc = yl(*)(-:1:nx,);

  /* xy plot, set up mapping stuff */
  if (type == "xy") {
    dtr = pi/180;
      
    /* set projection */
    if (!is_void(proj)) {

      /* add a wrap point to the data, outline, mask and coordinate arrays, 
         if x appears to be nearly periodic */
      dx = 720./nx;
      if ((xl(1)+360-dx < xl(0))  && (xl(0) < xl(1)+359.99)) {
	xlst = indgen(nx+1);
	xlst(0) = 1;
	data = data(xlst,..);
	if (!is_void(outline)) outline = outline(xlst,..);
	if (!is_void(mask)) mask = mask(xlst,..);
	xl  = xl(xlst);
	xl(0) += 360;
	nx = nx + 1;
	xc = xl(*)(,-:1:ny);
	yc = yl(*)(-:1:nx,);
      }

      /* set up grid arrays, if requested */
      NONVIEW = -9999.
      if (!is_void(xgrid)) {
	lx  = int(360/xgrid)+1;
	ly  = 179;
	xx = (xgrid * (indgen(lx) -1))(,-:1:ly);
	yx = span (-89, 89, ly)(-:1:lx,);
      }
      if (!is_void(ygrid)) {
	lx  = 361;
	ly  = int(180/ygrid) -1
	xy = span(0,360,lx)(,-:1:ly);
	yy = (-90 + ygrid * indgen(ly))(-:1:lx,);
      }

      /* project coordinate arrays */
      if (proj == "mollweide") {
	trans_mw, yc, xc, yc, xc;
	if (!is_void(xgrid)) trans_mw, yx, xx ,yx, xx;
	if (!is_void(ygrid)) trans_mw, yy, xy ,yy, xy; 

      } else if (proj == "polar") {
	mask = trans_ps(yc, xc, yc, xc, cirlat=cirlat, mask=mask);
	if (!is_void(xgrid)) trans_ps, yx, xx ,yx, xx;
	if (!is_void(ygrid)) trans_ps, yy, xy ,yy, xy; 
	/* set the plot limits to the bounding circle, if requested */
	if (!is_void(cirlat)) {
	  trans_ps(cirlat, 0, limy, limx);
	  lim = abs (limy, limx);
	  limits, -lim, lim, -lim, lim;
	}

      } else if (proj == "orthographic") {
	mask = trans_or(yc, xc, yc, xc, mask=mask);
	if (!is_void(xgrid)) trans_or, yx, xx ,yx, xx, nonview=NONVIEW;
	if (!is_void(ygrid)) trans_or, yy, xy ,yy, xy, nonview=NONVIEW; 

      } else {
	print, "projection not supported yet: ", + proj;
      }
    } 
  }

  /* check that initialization has taken place and a palette is available */
  local r, g, b;
  palette, r, g, b, query=1;
  if (numberof(r) == 0) palette,"earth.gp";

  /* label positioning variables */
  get_style, landscape, systems, legends, clegends;
  if (landscape == 0) {
    pgtop  = 1.0;
  } else {
    pgtop  = 0.775;
  }
  port  = viewport();
  frcen = port(avg:1:2);
  frtop = port(4);

  /* Set region masks if necessary */ 
  if (!is_void(mask)) {
    dmask = dimsof(mask);
    ddata = dimsof(data);
    if ((dmask(1) >= 2) && (ddata(1) >=2) && allof(dmask(2:3) == ddata(2:3))) {
    } else {
      print, "mask and plotting data dimensions do not match: ",
	dmask, ddata
      return;
    }
    if (nocont == 0) {
      iregc  = array(long, dmask(2), dmask(3));
//      if (dmask(1)==2) iregc(2:,2:) = mask(1:-1,1:-1);
      if (dmask(1)==2) iregc = mask;
    }
  }

  /* check for outline array */
  if (!is_void(outline)) {
    if (nameof(structof(outline)) == "NC_variable") {
      
    } else {
      doutl = dimsof(outline);
      ddata = dimsof(data);
      if ((doutl(1) == 2) && (ddata(1) >=2) && allof(doutl(2:3) == ddata(2:3))) {
      } else {
	print, "outline and plotting data dimensions do not match: ",
	  doutl, ddata
	  return;
      }
      xwo = xc(pcen,pcen);
      ywo = yc(pcen,pcen);
      irego  = array(long, doutl(2)+1, doutl(3)+1);
      irego(2:,2:) = outline;
    }
  }

  /* Set contours and min/max values for colour fill */
  if (is_void(levs)) levs = set_contours(data, cint=cint, nlevs=nlevs,spv=spv);
  if (is_void(width)) width=1;
  /* labelled contours */
  if (is_void(nlab)) {
    nlabf=1;
    nlabi=1;
  } else if (numberof(nlab)==1) {
    nlabf=1;
    nlabi=nlab;
  } else {
    nlabf = nlab(1);
    nlabi = nlab(2);
    if (nlabf==nlabi) nlabf= 0;
  }

  /* separarate contours */
  if (nocont == 0) {
    /* labelled and unlabelled contours */
    local lcont, ucont, lclt, uclt, lcge, ucge
    ucont = [];
    if (nlabi == 1) {
      lcont = levs;
    } else {
      for (i=1; i<=numberof(levs); i++) {
	if ((i-nlabf)%nlabi) {
	  ucont = grow(ucont,levs(i));
	} else
	  lcont = grow(lcont,levs(i));
      }
    }
    /* positive and negative contours */
    if (is_void(neg)) {
      lclt=[];
      uclt=[];
      lcge=lcont;
      ucge=ucont;
    } else {
      if (numberof(lcont)>0) {
	lclt=lcont(where(lcont <  0));
	lcge=lcont(where(lcont >= 0));
      }
      if (numberof(ucont)>0) {
	uclt=ucont(where(ucont <  0));
	ucge=ucont(where(ucont >= 0));
      }
    }
  }

  /* contour widths */
  lwidth = width(1);
  uwidth = width(0);

  /* set default mesh for potential use outside this routine */
  plmesh, yc, xc;

  /* draw top label first */
  if (noname == 0) {
    plt,  nc_var.long_name + " (" + nc_var.units + ")",
      frcen, frtop+.02, justify="CA", height=18.;
  }

  /* set up and draw corner label */
  if (nocornlab == 0) {
    cornlab  = "Variable: " + nc_var.name + "\nFile: " + nc_var.filename;
    if (strlen(nc_var.title)>0) cornlab += "\nTitle: " + nc_var.title;
    if (strlen(nc_var.case) >0) cornlab += "\nCase: " + nc_var.case;
    if (!is_void(*nc_var.t)) {
      if (is_void(t)) list = indgen(numberof(*nc_var.t));
      else list = tlist;
      cornlab += "\nTime: " + pr1((*nc_var.t)(list(1)));
      if (numberof(list) > 1) {
	cornlab += " - " + pr1((*nc_var.t)(list(0)));
      }
    }
    if (!is_void(*nc_var.date)) {
      if (is_void(t)) list = indgen(numberof(*nc_var.date));
      else list = tlist;
      date = (*nc_var.date)(list(1));
      yl = int(date/10000); m = int(date/100)%100; d = date%100;
      cornlab += "\nDate: " + pr1(yl) + "/" + pr1(m) + "/" + pr1(d);
      if (numberof(list) > 1) {
	date = (*nc_var.date)(list(0));
	yl = int(date/10000); m = int(date/100)%100; d = date%100;
	cornlab +=" - " +  pr1(yl) + "/" + pr1(m) + "/" + pr1(d);
      }
    }
    plt, cornlab, .03, pgtop, justify="LT", height=12.;
  }

  // set the palette if requested and colors indices not input
  if (nofill == 0) {
    if (numberof(colors)==0) {
      if (setpal < 0) {
	reverse=1;
	setpal = abs(setpal);
      }
      if (setpal > 1) {
	colors = plnpal(levs,extend=1,reverse=reverse);
      } else if (setpal == 1) {
	colors = plnpal(levs,reverse=reverse);
      }
    }

    // draw the color bar
    if((nocbar <= 0)) {
      vert = 1;
      if (nocbar < 0) vert=0;
      if (numberof(colors)>0) colinds = char(colors(1:numberof(levs)+1));
      if (nlabi >1) {
	color_bar,levs,colinds,vert=vert,adjust=cbar,labs=[nlabf,nlabi];
      } else {
	color_bar,levs,colinds,vert=vert,adjust=cbar;
      }
    }
  }

  /* turn on animation mode */
  if (is_void(anim)) {
    nanimframe = 1;
    nanimloop  = 1;
  } else {
    fma;
    animate, 1;
    data3 = data;
    nanimframe = (dimsof(data))(0);
    nanimloop  = anim;
  }

  /* loop for animation */
  for (nan=1; nan<=nanimloop; nan++){
    for (n=1; n<=nanimframe; n++) {
  
      /* get data slice for animation */ 
      if (!is_void(anim)) {
//	print, "frame: ", n, "of ", nanimframe, "loop", nan, "of", anim;
	data = data3(..,n);
	if (!is_void(mask)) {
//	  if ((nocont==0) && (dmask(1)==3)) iregc(2:,2:) = mask(1:-1,1:-1,n);
	  if ((nocont==0) && (dmask(1)==3)) iregc = mask(,,n);
	}
      }

      /* color fill */
      if (nofill == 0) {
	if (numberof(colors)>0) colinds = char(colors);
	plfcs, data, yc, xc, iregc, levs=levs, colors=colinds,region=1,spv=spv;
      }

      /* draw contours */ 
      if (nocont == 0) {
	if (numberof(lclt) > 0) 
	  plc_svn, data, yc, xc, iregc, spv=spv, levs=lclt, 
	    width=lwidth, concol=concolor, type=neg, lab=1,
	    ngradmax=1.E9, opaque=opaque, height=con_lab_height;
	if (numberof(uclt) > 0) 
	  plc_svn, data, yc, xc, iregc, spv=spv, levs=uclt, 
	    width=lwidth, concol=concolor, type=neg, lab=0;
	if (numberof(lcge) > 0) 
	  plc_svn, data, yc, xc, iregc, spv=spv, levs=lcge, 
	    width=lwidth, concol=concolor, lab=1,
	    ngradmax=1.E9, opaque=opaque, height=con_lab_height;
	if (numberof(ucge) > 0) 
	  plc_svn, data, yc, xc, iregc, spv=spv, levs=ucge, 
	    width=lwidth, concol=concolor, lab=0;
      }

      // draw extrema labels if requested
      if (mxmn != 0) {
	maxnum = [];
	if (mxmn != 1) maxnum = mxmn;
	lab_maxmin,data,xc,yc,spv=spv,maxnum=maxnum;
      }

      /* draw outline */
      if(!is_void(outline))plm, ywo, xwo, irego, boundary=1, width=4.;

      /* draw xgrid */
      if(!is_void(xgrid)) {
	for (l=1; l<= dimsof(xx)(2); l++) {
	  yxl = yx(l,);
	  xxl = xx(l,);
	  if (anyof(yxl == NONVIEW)) {
	    list = where(yxl != NONVIEW);
	    if(numberof(list)>1) {
	      brks = where(list(2:) != list(:-1)+1);
	      i1 = list(1); 
	      i2 = list(0);
	      if (numberof(brks)>0) {
		i1 = grow(i1,list(brks+1));
		i2 = grow(list(brks),i2);
	      }
	      for (n=1; n<=numberof(i1); n++) {
		plg, yxl(i1(n):i2(n)), xxl(i1(n):i2(n)), marks=0, type="dash";
	      }
	    }
	  } else {
	    plg, yxl, xxl, marks=0, type="dash";
	  }
	}
      }

      /* draw ygrid */
      if(!is_void(ygrid)) {
	for (l=1; l<= dimsof(yy)(3); l++) {
	  yyl = yy(,l);
	  xyl = xy(,l);
	  if (anyof(yyl == NONVIEW)) {
	    list = where(yyl != NONVIEW);
	    if(numberof(list)>1) {
	      brks = where(list(2:) != list(:-1)+1);
	      i1 = list(1); 
	      i2 = list(0);
	      if (numberof(brks)>0) {
		i1 = grow(i1,list(brks+1));
		i2 = grow(list(brks),i2);
	      }
	      for (n=1; n<=numberof(i1); n++) {
		plg, yyl(i1(n):i2(n)), xyl(i1(n):i2(n)), marks=0, type="dash";
	      }
	    }
	  } else {
	    plg, yyl, xyl, marks=0, type="dash";
	  }
	}
      }

      /* do frame advance, in animation mode only */
      if (!is_void(anim)) {
	if (!is_void(*nc_var.t)) 
	  plt, "Time: " + pr1((*nc_var.t)(n)), port(1), frtop-.02, 
	    justify="LA", height=14;
	if (!is_void(*nc_var.date)) 
	  plt, "Date: " + pr1((*nc_var.date)(n)), port(2), frtop-.02, 
	    justify="RA", height=18;
	fma;
      }
    }
  }
  /* turn off animation mode */
  if (!is_void(anim)) animate, 0;
  
  return [];
}

func make_sublist (list, name, dim, rng=)
{
  if (dim == 1) return [1];

  if (is_void(rng)) {
    sublist = indgen(numberof(list));
  } else {
    if (numberof(rng) == 1 ) {
      sublist = [abs(list - rng)(mnx)];
    } else if (numberof(rng) == 2 ) {
      if (numberof(list) == 1) {
	sublist = [1];
      } else {
	rngmin = rng(min);
	rngmax = rng(max);
	sublist = where ((rngmin<=list) & (list <= rngmax));
      }
    }
  }
  print, name, "values:", list(sublist);
  if (numberof(sublist)==0) {
    print, "empty coordinate range:", rng, "coordinate values", list;
  }
  
  return sublist;
}

func set_contours(data,cint=,nlevs=,spv=)
/* DOCUMENT levs = set_contours (data, cint=CINT, nlevs=NLEV, spv=SPV)
   return contour levels between the min and max values in the data array.
     if CINT specified, use for the interval,
     if NLEVS specified and CINT not specified, generate approximately 
       NLEVS levels (default 10), requiring an interval of 1,2,or 5
     if SPV specified exclude from min/max determination
*/
{
  // find the min and max values in input data
  if (is_void(spv)) {
    fmax = max(data);
    fmin = min(data);
  } else {
    list = where(data!=spv);
    if (numberof(list) > 0) {
      fmax = max(data(list));
      fmin = min(data(list));
    } else {
      print, "set_contours: only special values in data";
    }
  }

  // check for uniform field
  if (fmax == fmin) {
    print, "set_contours: uniform input data, value=", fmin;
    if (fmax != 0) {
      return [0.9*fmax,fmax,1.1*fmax];
    } else {
      return [-0.1,0.0, 0.1];
    }
  }

  // determine a contour interval, if not specified
  if (is_void(cint)) {
    if(is_void(nlevs)) nlevs = 10;
    okints = [0.1, 0.2, 0.5,1,2,5,10,20]; // acceptable contour intervals
    zinc = (fmax-fmin)/nlevs;
    pow = int(log10(zinc));
    cints = okints*(10.^pow);
    nlevsout = (fmax-fmin)/cints;
    nlevbest = abs(nlevsout-nlevs)(mnx);
    cint = cints(nlevbest);
  }

  // min and max contour levels and number of levels
  cmin = int(fmin/cint)*cint; if (fmin > 0) cmin = cmin+cint;
  cmax = int(fmax/cint)*cint;
  nlev = int((cmax - cmin) / cint + .499999) + 1;

  // this happens if cint was specified inappropriately
  if (nlev <= 1 || nlev >= 70) {
    print, "inappropriate contour interval";
    print,fmin,fmax,cint,cmin,cmax,nlev;
    nlev = 10; cmin = fmin; cmax = fmax;
  }

  // generate the contour levels and make sure zero is zero
  flev = span(cmin,cmax,nlev);
  list = where(abs(flev/cint) < 0.000001);
  if(numberof(list)==1)flev(list) = 0.;
  return flev;
}


func trans_mw (y, x, &yt, &xt)
/* DOCUMENT trans_mw (y, x, &yt, &xt)

   Transform a pair of arrays (latitude, longitude) defining a 2d mesh to screen
   coordinates, using a mollweide projection

*/
{
  local dtr;
  dtr = pi/180;
  xt = (x-180.)/90.*cos(y*dtr);
  yt = cos((90.-y)*dtr);
}

func trans_ps (y, x, &yt, &xt, cirlat=, mask=)
/* DOCUMENT trans_ps (y, x, &yt, &xt)
   or mask =trans_ps (y, x, &yt, &xt, cirlat=, mask=)

   Transform a pair of arrays (latitude, longitude) defining a 2d mesh to screen
   coordinates, using a polar steriographic projection whose centre is defined by 
   (polat, polon). Is cirlat -s not void it defines a bounding latitude outside 
   which the returned mask is set to zero.

*/
{
  local cirrad, dtr, list, rad, phi;
  if (is_void(polat)) polat = 90.;
  if (is_void(polon)) polon =  0.;
  dtr = pi/180;
  rad = tan(abs(y - polat)*dtr/2);
  phi = (x-polon)*dtr;
  if (polat < 0) phi = -phi;
  xt  = rad * cos(phi);
  yt  = rad * sin(phi);
  /* mask out the area outside of a bounding circle, if requested */
  if (!is_void(cirlat)) {
    if (is_void(mask)) mask = array (1, dimsof(x));
    cirrad = abs(tan((polat-cirlat)*dtr/2));
    list = where(rad>cirrad);
    if (numberof(list)>0) mask(list) = 0;
  }
  return mask;
}

func trans_or (y, x, &yt, &xt, mask=, nonview=)
/* DOCUMENT trans_or (y, x, &yt, &xt, nonview=))
   or mask =trans_or (y, x, &yt, &xt, mask=)

   Transform a pair of arrays (latitude, longitude) defining a 2d mesh to screen
   coordinates, using an orthographic projection whose centre is defined by 
   (polat, polon). If nonview is nonnil, then yt and xt for nonviewable points
   are returned as its value. If nonview is nil, then the function returns a mask 
   set to zero where the point is not viewable.

*/
{
  local c, dtr, list, xr, yr, xr0, yr0;
  if (is_void(polat)) polat = 90.;
  if (is_void(polon)) polon =  0.;
  dtr = pi/180;
  xr0 = polon*dtr;
  yr0 = polat*dtr;
  xr = (x-polon)*dtr;
  if (polat < 0) xr = -xr;
  yr = y*dtr;
  xt = cos(yr)*sin(xr);
  yt = cos(yr0)*sin(yr) - sin(yr0)*cos(yr)*cos(xr);
  c = sin(yr0)*sin(yr) + cos(yr0)*cos(yr)*cos(xr);
  list = where(c<=0.01);
  if (is_void(nonview)) {
    /* mask the area which is not viewable */
    if (is_void(mask)) {
      msk = array (1, dimsof(x));
    } else {
      msk = mask;
    }
    if (numberof(list)>0) msk(list) = 0;
    return msk;
  } else {
    /* reset the nonviewable points */
    if (numberof(list)>0) yt(list)=xt(list)=nonview;
  }
}

func plnpal(levs,extend=,white=,spvcol=,reverse=)
/* DOCUMENT cols=plnpal(levs,center=,extend=,white=,spvcol=,reverse=)
   
   Generate a blue-white-red palette if the contours levels span zero
   otherwise, generate a green/yellow/orange palette.
   Return the list of color indices added to the palette.
   If extend=1, the current palette is extended, otherwise it is replaced.
   If white nonnil and nonzero, it is used as the color value 
   (usually ~220 for grey)
   If white =0, there is a direct transition from blue to red
   If spvcol nonnil, the special value color will be set to spvcol, 
      which may be either a single value or rgb triplet in the range [0,255]
   If reverse is nil, reverse the palette if all levels <=0, otherwise reverse if true.
*/
{
  if ((levs(1)<0) && (levs(0)>0)) {
    fcol = bwrpal(levs,extend=extend,center=0,white=white,spvcol=spvcol,reverse=reverse);
  } else {
    if (levs(0) <=0) reverse=1;
    nf = numberof(levs)+1;
    if (nf<3) {
      nn = (nf-1);
      fcol=glpal([255,255,255],[255,255,255],1,
		 [230,230,255],nn,extend=extend,reverse=reverse,spvcol=spvcol);
    } else if (nf<5) {
      nn = (nf-2);
      fcol=glpal([255,255,255],[255,255,255],1,
		 [230,230,255],1,[100,200,100],nn,extend=extend,reverse=reverse,spvcol=spvcol);
    } else if (nf<7) {
      nn = (nf-3)/2;
      nx = (nf-3 - nn);
      fcol=glpal([255,255,255],[255,255,255],1,
		 [230,230,255],1,[100,200,100],nn,
		 [255,150,000],1,[255,220, 80],nx,extend=extend,reverse=reverse,spvcol=spvcol);
    } else {
      nn = (nf-4)/3;
      nx = (nf-4 - 2*nn);
      fcol=glpal([255,255,255],[255,255,255],1,
		 [230,230,255],1,[100,230,100],nx,
		 [170,210,000],1,[255,255,100],nn,
		 [255,220, 80],1,[255,150,000],nn,extend=extend,reverse=reverse,spvcol=spvcol);
      //		 [255,150,000],1,[255,220, 80],nx,extend=extend,reverse=reverse);
    }
  }
  return fcol;
}

func bwrpal(levs,center=,extend=,white=,spvcol=,reverse=)
/* DOCUMENT cols=bwrpal(levs,center=,extend=,white=,spvcol=,reverse=)
   
   Generate a blue-white-red palette for the contours levels in levs and
   return the list of color indices added to the palette. The palette will 
   contain 2 "white" entries between blue and red.
   If center is nonnil, white is centered on its value
   If extend=1, the current palette is extended, otherwise it is replaced.
   If white nonnil and nonzero, it is used as the color value 
   (usually ~220 for grey)
   If white =0, there is a direct transition from blue to red
   If spvcol nonnil, the special value color will be set to spvcol, 
      which may be either a single value or rgb triplet in the range [0,255]

   If reverse!=0, reverse the output palette
*/
{
  // Set number of blue (n1) and red (n2) colours and 
  // generate the colour index list (starts with zero).
  // The number of fill intervals is 1 longer than the number of levels 
  // and an extra color is included for special values
  nlev=numberof(levs);
  i=indgen(nlev+2)-1;
  if(is_void(center)) {
    n1 = int(nlev/2);
    n2 = nlev-n1-1;
  } else {
    nz = abs(levs-center)(mnx);
    n1 = nz-1;
    n2 = nlev-nz;
  }
  /* set default "white" value or adjust size of blue and red values */
  if(is_void(white)) {
    white= 255;
  } else if(white==0) {
    n1 += 1;
    n2 += 1;
  }

  /* generate equally spaced blue entries */
  local ra,ga,ba,rb,gb,bb,rc,gc,bc;
  if (n1 >= 2) {
    ra = span( 75,230,n1); 
    ga = span( 75,230,n1);
    ba = [255](-:1:n1);
  } else if (n1 == 1) {
    ra = 75; ga = 75; ba = 255;
  }
  /* 2 white entries */
  rb = [white,white]; 
  gb = [white,white];
  bb = [white,white];
  /* generate equally spaced red entries */
  if (n2 >= 2) {
    rc = [255](-:1:n2);
    gc = span(230,75,n2);
    bc = span(230,75,n2);
  } else if (n2 ==1) {
    rc = 255; gc = 75; bc = 75;
  }

  /* assemble the full list from the 2 or 3 sublists */
  if (white != 0) {
    r=char(grow(ra,rb,rc));
    g=char(grow(ga,gb,gc));
    b=char(grow(ba,bb,bc));
  } else {
    r=char(grow(ra,rc));
    g=char(grow(ga,gc));
    b=char(grow(ba,bc));
  }    

  // reverse the palette if requested
  if (!is_void(reverse) && reverse) {
    r=r(0:1:-1); g=g(0:1:-1); b=b(0:1:-1);
  }

  // add the special value color
  if (is_void(spvcol)) {
    spvc = [150,150,150];
  } else if (numberof(spvcol==1)) {
    spvc = [spvcol,spvcol,spvcol];
  } else {
    spvc = spvcol;
  }
  r = grow(r,char(spvc(1)));
  g = grow(g,char(spvc(2)));
  b = grow(b,char(spvc(3)));

  /* if extending the palette, retrieve the existing palette, grow the 
     rgb lists, and shift the colour index list */
  if (extend==1) {
    palette,ro,go,bo,query=1;
    r=grow(ro,r); g=grow(go,g); b=grow(bo,b);
    i+=numberof(ro);
  }

  /* set the palette and return the colour indexes */
  palette,r,g,b;
  return i;
}

func glpal (.., extend=,spvcol=,reverse=)
/* DOCUMENT cols = glpal(rgb0, rgb1, n1, [rgb2, n2,...,] 
                         extend=,spvcol=,revers=)

   Generate a new palette which scales linearly, using nN values,
   between the rgbN and rgbN-1 values and return the colour indexes
   used. An arbitrary number of [rgbN, nN] pairs may be input and the
   palette will scale linearly from one to the next. The total number
   of entries will be sum(nN). Note the total number of arguments must
   be odd. The rgbN must be length 3 arrays of integers (long or char)
   representing [red, green, blue] intensities in the range [0, 255].

   If extend is void or zero, replace the existing palette, otherwise
   add the new entries to the end.

   If spvcol nonnil, the special value color will be set to spvcol, 
      which may be either a single value or rgb triplet in the range [0,255]

   If reverse!=0, reverse the output palette
*/
{
  /* check the number of input variable arguments */
  narg = more_args();
  if (narg < 3) {
    print, "glpal: at least 3 input arguments required", narg, "found";
    error;
  }
  if (narg%2 != 1) {
    print, "glpal: the number of arguments must be odd, but =", narg;
  }
  nrgb = narg / 2;

  /* loop over the input rgb values, generating colour values */
  local r, g, b;
  rgb = next_arg();
  r2 = rgb(1); g2 = rgb(2); b2 = rgb(3);
  for (irgb=1; irgb<=nrgb; irgb++) {
    r1 = r2; g1 = g2; b1 = b2;
    rgb = next_arg();
    r2 = rgb(1); g2 = rgb(2); b2 = rgb(3);
    n = next_arg();
    if (n == 1) {
      r = grow(r, r2);
      g = grow(g, g2);
      b = grow(b, b2);
    } else {
      r = grow(r, span(r1,r2,n+1)(2:));
      g = grow(g, span(g1,g2,n+1)(2:));
      b = grow(b, span(b1,b2,n+1)(2:));
    }
  }
  r = char(r); g = char(g); b = char(b);

  // reverse the palette if requested
  if (!is_void(reverse) && reverse) {
    r=r(0:1:-1); g=g(0:1:-1); b=b(0:1:-1);
  }

  // add the special value color
  if (is_void(spvcol)) {
    spvc = [150,150,150];
  } else if (numberof(spvcol==1)) {
    spvc = [spvcol,spvcol,spvcol];
  } else {
    spvc = spvcol;
  }
  r = grow(r,char(spvc(1)));
  g = grow(g,char(spvc(2)));
  b = grow(b,char(spvc(3)));

  // generate the index list
  i = indgen(numberof(r))-1;

  /* if extending the palette, retrieve the existing palette, grow the 
     rgb lists, and shift the colour index list */
  if (!is_void(extend) && extend!=0) {
    palette,ro,go,bo,query=1;
    r=grow(ro,r); g=grow(go,g); b=grow(bo,b);
    i+=numberof(ro);
  }

  /* set the palette and return the colour indexes */
  palette,r,g,b;
//  "r"; r;
//  "g"; g;
//  "b"; b;
  return i;
}

func lab_maxmin (d,x,y,cmin=,cmax=,spv=,maxnum=) {
/* DOCUMENT lab_maxmin,d,x,y,cmin=CMIN,cmax=CMAX,spv=SPV,maxnum=NUM;
     Labels relative maxima and minima in the interior of the 
     2 dimensional input array d. x and y are coordinate arrays
     of the same dimensions of d.

     Extrema greater than CMIN or less than CMAX will not be labelled.
     Points with values (or neighbours with values) of SPV will not be
     considered extrema. 

     At most NUM max and mins will be labelled (default 10). If more 
     are found, then the NUM largest maxima and smallest minima will be 
     labelled. If NUM<0 then abs(NUM) values equally spaced through 
     the array will be labelled. This should space the labels around
     the plot.
*/

  local dd, ddims,nxd,nyd;
  dd = d(2:-1,2:-1);  // only test on the interior points

  ddims = dimsof(d);
  nxd = ddims(2);
  nyd = ddims(3);

  if (is_void(cmin)) cmin = max(d); // label all maxs and mins
  if (is_void(cmax)) cmax = min(d); // by default
  if (is_void(spv )) spv  = 2*max(d); // set an out of range value
  if (is_void(maxnum))maxnum = 10    //

  // make a list of maxima > cmax
  list = where2( (dd > cmax) & dd != spv &
		 (dd > d(1:-2,2:-1)) & (d(1:-2,2:-1) != spv) & 
		 (dd > d(3: 0,2:-1)) & (d(3: 0,2:-1) != spv) & 
		 (dd > d(2:-1,1:-2)) & (d(2:-1,1:-2) != spv) & 
		 (dd > d(2:-1,3:-0)) & (d(2:-1,3:-0) != spv) & 
		 (dd > d(1:-2,1:-2)) & (d(1:-2,1:-2) != spv) & 
		 (dd > d(3: 0,1:-2)) & (d(3: 0,1:-2) != spv) & 
		 (dd > d(1:-2,3: 0)) & (d(1:-2,3: 0) != spv) & 
		 (dd > d(3: 0,3:-0)) & (d(3: 0,3:-0) != spv) );
  // if the list is not empty, plot labels on the extrema
  if(numberof(list) > 0) {
    list +=1;  // shift the indexing to the original array
    dlst = d((list(2,)-1)*nxd + list(1,)); // extrema values
    //    dlst;
    nlst = numberof(dlst);
    if(nlst>abs(maxnum)) {    // if the list is too long, sort it
      nlab = abs(maxnum);
      if (maxnum>0) {
	lst2 = sort(dlst)(0:1-nlab:-1); // nlab largest maxima
      } else {
	lst2 = long(span(1,nlst,nlab)); // space through list
	//	lst2;
      }
      dlst = dlst(lst2);
      list = list(,lst2);
      //      dlst;
    } else {
      nlab = nlst;
    }
    for (n=1; n<=nlab; n++) {
      lab = swrite(format="%.5g",dlst(n));
      i1 = list(1,n); i2 = list(2,n);
      plt,lab,x(i1,i2),y(i1,i2),tosys=1,justify="CH",
	height=10,color="red",opaque=0;
    }
  }
  
  // make a list of minima < cmaxmx
  list = where2( (dd < cmin) & dd != spv &
		 (dd < d(1:-2,2:-1)) & (d(1:-2,2:-1) != spv) & 
		 (dd < d(3: 0,2:-1)) & (d(3: 0,2:-1) != spv) & 
		 (dd < d(2:-1,1:-2)) & (d(2:-1,1:-2) != spv) & 
		 (dd < d(2:-1,3:-0)) & (d(2:-1,3:-0) != spv) & 
		 (dd < d(1:-2,1:-2)) & (d(1:-2,1:-2) != spv) & 
		 (dd < d(3: 0,1:-2)) & (d(3: 0,1:-2) != spv) & 
		 (dd < d(1:-2,3: 0)) & (d(1:-2,3: 0) != spv) & 
		 (dd < d(3: 0,3:-0)) & (d(3: 0,3:-0) != spv) );
  // if the list is not empty, plot labels on the extrema
  if(numberof(list) > 0) {
    list +=1;  // shift the indexing to the original array
    dlst = d((list(2,)-1)*nxd + list(1,)); // extrema values
    //    dlst;
    nlst = numberof(dlst);
    if(nlst>abs(maxnum)) {    // if the list is too long, sort it
      nlab = abs(maxnum);
      if (maxnum>0) {
	lst2 = sort(dlst)(1:nlab);      // nlab smallest minima
      } else {
	lst2 = long(span(1,nlst,nlab)); // space through list
	//	lst2;
      }
      dlst = dlst(lst2);
      list = list(,lst2);
      //      dlst;
    } else {
      nlab = nlst;
    }
    for (n=1; n<=nlab; n++) {
      lab = swrite(format="%.5g",dlst(n));
      i1 = list(1,n); i2 = list(2,n);
      plt,lab,x(i1,i2),y(i1,i2),tosys=1,justify="CH",
	height=10,color="blue",opaque=0;
    }
  }
}

func plg3 (v1, v2, v3, x, lab1=, lab2=, lab3=, labx=, color=, width=)
/* DOCUMENT plg3, v1, v2, v3, x, lab1=, lab2=, lab3=, labx, color=, width=

   Makes line plots on 3 panels of 3 input variables, versus x, with an optional label
   on each panel and on the bottom axis. The color keyword is passed to plg.

*/
{
  if (is_void(color)) color="black"
  plsys,1;
  plg,v1,x,marks=0,color=color,width=width;
  if (is_void(lab1)) lab1 = "all surfaces";
  plt, lab1, .1, .91, height=18.;

  plsys,2;
  plg,v2,x,marks=0,color=color,width=width;
  if (is_void(lab2)) lab2 = "land only";
  plt, lab2, .1, .63, height=18.;

  plsys,3;
  plg,v3,x,marks=0,color=color,width=width;
  if (is_void(lab3)) lab3 = "ocean/ice";
  plt, lab3, .1, .35, height=18.;

  if (is_void(labx)) labx = "year";
  plt, labx, .4, .09, justify="CA",height=18;
}

func plg1_init (viewport=, landscape=)
/* DOCUMENT plg1_init 
   read the style file and make 1 systems (viewports) in landscape mode
*/
{
  read_style, "boxed.gs", glandscape, systems, legends, clegends;
  if (is_void(landscape)) landscape = 1;
  if (is_void(viewport))  viewport  = [.15,.95,.1,.6];
  systems(1:).ticks.horiz.nMajor=20;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=40;
  systems.viewport = viewport;
  set_style, landscape, systems, legends, clegends;
}

func plg2_init
/* DOCUMENT plg2_init 
   read the style file and make 2 systems (viewports) in portrait mode
*/
{
  read_style, "boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:2);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=40;
  systems(1).viewport = [.15,.75,.53,.90];
  systems(2).viewport = [.15,.75,.11,.48];
  set_style, 0, systems, legends, clegends;
}


func plg3_init
/* DOCUMENT plg3_init 
   read the style file and make 2 systems (viewports) in portrait mode
*/
{
  read_style, "boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:3);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems(1).viewport = [.15,.75,.66,.89];
  systems(2).viewport = [.15,.75,.38,.61];
  systems(3).viewport = [.15,.75,.10,.33];
  set_style, landscape, systems, legends, clegends;
}

func pln_init (viewport=, landscape=)
{
  /* Reset current style parameters */
  read_style, "boxed.gs", glandscape, systems, legends, clegends;
  if (is_void(landscape)) landscape = glandscape;
  if (is_void(viewport))  viewport  = systems.viewport;
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems.viewport = viewport;
  set_style, landscape, systems, legends, clegends;
}
func pln_init_map (viewport=, landscape=)
{
  /* Reset current style parameters */
  read_style, "boxed.gs", glandscape, systems, legends, clegends;
  if (is_void(landscape)) landscape = 1;
  if (is_void(viewport))  viewport  = [.075,.875,.2,.6];
  systems.viewport = viewport;
  systems(1).ticks.horiz.nMajor=10;
  systems(1).ticks.horiz.nMinor=100;
  systems(1).ticks.vert.nMajor=10;
  systems(1).ticks.vert.nMinor=100;
  set_style, landscape, systems, legends, clegends;
}
func pln2_init
{
  /* read the style file and make 2 systems (viewports) */
  read_style, "boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:2);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems(1).viewport = [.08,.68,.540,.840];
  systems(2).viewport = [.08,.68,.140,.440];
  set_style, landscape, systems, legends, clegends;
}
func pln3_init
{
  /* read the style file and make 3 systems (viewports) */
  read_style, "boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:3);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems(1).viewport = [.1,.6,.675,.925];
  systems(2).viewport = [.1,.6,.375,.625];
  systems(3).viewport = [.1,.6,.075,.325];
  set_style, landscape, systems, legends, clegends;
}
func pln4_init
{
  /* read the style file and make 4 systems (viewports) */
//  read_style, "~boville/Gist/boxed-noticks.gs", landscape, systems, legends, clegends;
  read_style, "~boville/Gist/boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:4);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems(1).viewport = [.12, .42, .40, .70];
  systems(2).viewport = [.12, .42, .06, .36];
  systems(3).viewport = [.52, .82, .40, .70];
  systems(4).viewport = [.52, .82, .06, .36];
  set_style, 1, systems, legends, clegends;
}
func pln6_init
{
  /* read the style file and make 6 systems (viewports) */
  read_style, "boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:6);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems(1).viewport = [.0909, .3409, .6591, .9091];
  systems(2).viewport = [.0909, .3409, .3636, .6136];
  systems(3).viewport = [.0909, .3409, .0682, .3182];
  systems(4).viewport = [.4318, .6818, .6591, .9091];
  systems(5).viewport = [.4318, .6818, .3636, .6136];
  systems(6).viewport = [.4318, .6818, .0682, .3182];
  set_style, landscape, systems, legends, clegends;
}

func pln6_init_nocbar
{
  /* read the style file and make 6 systems (viewports) */
  read_style, "boxed.gs", landscape, systems, legends, clegends;
  systems=systems(-:1:6);
  systems(1:).ticks.horiz.nMajor=10;
  systems(1:).ticks.horiz.nMinor=100;
  systems(1:).ticks.vert.nMajor=10;
  systems(1:).ticks.vert.nMinor=100;
  systems(1).viewport = [.10, .39, .6591, .9091];
  systems(2).viewport = [.10, .39, .3636, .6136];
  systems(3).viewport = [.10, .39, .0682, .3182];
  systems(4).viewport = [.44, .73, .6591, .9091];
  systems(5).viewport = [.44, .73, .3636, .6136];
  systems(6).viewport = [.44, .73, .0682, .3182];
  set_style, landscape, systems, legends, clegends;
}