diff --git a/methods/islet/figures/amb3.hy b/methods/islet/figures/amb3.hy
index 6382bc9..4075817 100644
--- a/methods/islet/figures/amb3.hy
+++ b/methods/islet/figures/amb3.hy
@@ -1,7 +1,7 @@
 ;;; Collection of utils.
 
 ;;(require [hy.contrib.walk [let]])
-(import sys copy math os re)
+(import sys copy math os re [importlib [reload]])
 
 ;;   when passing kwargs to another function like pl.plot, the dictionary should
 ;; be like {'option value}, not {:option value}.
@@ -23,7 +23,6 @@
 (defn symbol [string] (HySymbol string))
 
 (defmacro/g! svifn [&rest args]
-  "setv if none"
   `(do ~@(map (fn [e]
                 `(if (none? ~(first e)) (setv ~(first e) ~(last e))))
               (zip (cut args 0 (len args) 2)
@@ -81,7 +80,7 @@
   `(print (.format ~@args) :end "" :file ~fid))
 
 (defmacro prc [sym]
-  `(prf ~(+ (name sym) " {}") ~sym))
+  `(print ~sym "|" (eval (read-str ~sym))))
 
 (defmacro mapply [func &rest args]
   "Apply func. Works for keyword args, unlike apply. (Probably broken in many
@@ -149,47 +148,6 @@
        (setv last? (= ~g!i ~g!last))
        ~@body)))
 
-;; Use this instead of macroexpand to get output stripped of the Hy object
-;; ctors.
-(defn ppme [quoted-form]
-  (sv sym-dict {} b (Box) b.sym-num 0
-      b.after-open True)
-  (defn sym [e]
-    (unless (in e sym-dict)
-      (assoc sym-dict e (.format "sym-{:d}" b.sym-num))
-      (inc! b.sym-num))
-    (get sym-dict e))
-  (defn prl [e ldelim rdelim]
-    (prfno "{:s}" (+ (if b.after-open "" " ") ldelim))
-    (sv b.after-open True)
-    (for [li e] (rec li))
-    (prfno "{:s}" rdelim))
-  (defn atom? [e]
-    (in "quote" (first e)))
-  (defn rec [e]
-    (setv t (type e))
-    (case/in t
-             [[hy.models.HyFloat HyInteger] (print (+ " " (str e)) :end "")]
-             [[HyExpression]
-              (if (atom? e)
-                (for [li e] (rec li))
-                (prl e "(" ")"))]
-             [[HyList] (prl e "[" "]")]
-             [[HyString] (print (.format " \"{:s}\"" e) :end "")]
-             [:else
-              (unless b.after-open (print " " :end ""))
-              (sv b.after-open False)
-              (cond [(in "keyform" e)
-                     (print (sym e) :end "")]
-                    [(in "quote" e)
-                     (print "'" :end "")
-                     (sv b.after-open True)]
-                    [:else (print e :end "")])]))
-  (setv h (macroexpand quoted-form))
-  (print h)
-  (rec h)
-  (print))
-
 (defclass Box []
   "A box to hold values to be written in closures."
   (defn --repr-- [me]
@@ -357,6 +315,11 @@
   (/ (sum (map (fn [e] (** (- e mu) 2)) coll))
      (len coll)))
 
+(defn cumsum [coll]
+  (setv c (list (copy.deepcopy coll)))
+  (for [i (range 1 (len c))] (+= (get c i) (nth c (dec i))))
+  c)
+
 (defn cross-prod [x y]
   (defn one [i0 i1]
     (- (* (get x i0) (get y i1)) (* (get x i1) (get y i0))))
@@ -397,27 +360,35 @@
    (except [e Exception]
      (print ln))))
 
-(defn ooa [y &optional [xfac 2] [x None]]
+(defn calc-ooa [y &optional [xfac 2] [x None]]
+  (defn bad [e] (or (npy.isnan e) (npy.isinf e)))
   (setv r [])
   (for [i (range (dec (len y)))]
+    (when (or (bad (nth y i)) (bad (nth y (inc i))))
+      (prf "calc-ooa: i {} y(i) {} y(i+1) {}" i (nth y i) (nth y (inc i))))
     (.append r (/ (- (math.log (get y i)) (math.log (get y (inc i))))
                   (if x
-                    (- (math.log (get x i)) (math.log (get x (inc i))))
+                    (- (math.log (get x (inc i))) (math.log (get x i)))
                     (math.log xfac)))))
   r)
 
-(defn ooa-from-file [fname fieldno &optional anchor]
+(defn ooa-from-file [fname xanchor xfldno yanchor yfldno]
   "Read text from file fname, optionally scan for lines beginning with anchor,
   and read symbol fieldno, starting from 0. Return a list of OOAs."
   (sv txt (.split (readall fname) "\n")
-      errs [])
+      x [] y [])
   (for [ln txt]
-    (unless (none? anchor)
-      (when (or (< (len ln) (len anchor))
-                (!= anchor (cut ln 0 (len anchor))))
-        (continue)))
-    (.append errs (float (get (.split ln) fieldno))))
-  (, errs (ooa errs)))
+    (cond [(and (>= (len ln) (len xanchor))
+                (= xanchor (cut ln 0 (len xanchor))))
+           (.append x (float (get (.split ln) xfldno)))]
+          [(and (>= (len ln) (len yanchor))
+                (= yanchor (cut ln 0 (len yanchor))))
+           (.append y (float (get (.split ln) yfldno)))]))
+  (, x y (ooa y :x x)))
+
+(defn sypd [procs calls time &optional [calls-per-day 48]]
+  (/ (/ calls procs calls-per-day 365)
+     (/ time (* 24 3600))))
 
 (defn single? [coll]
   (or (not (coll? coll))
@@ -521,8 +492,10 @@
  (do
   (import [numpy :as npy] ctypes)
 
+  (defmacro getc [vs i] `(get ~vs (, (s-all) ~i)))
+
   (defn array-range [&rest args]
-    (npy.array (list (apply range args)) :dtype int))
+    (npy.array (list (range #*args)) :dtype int))
 
   (defn array-if-not [A &optional [dtype float]]
     (unless (= (type A) npy.ndarray)
@@ -540,6 +513,10 @@
   (defn ones-row-vec [n] (npy.ones (, 1 (pod-or-len n))))
   (defn ones-col-vec [n] (npy.ones (, (pod-or-len n) 1)))
 
+  (defn xy->XY [x y]
+    (, (npy.dot (col-vec (npy.ones (len y))) (row-vec x))
+       (npy.dot (col-vec y) (row-vec (npy.ones (len x))))))
+
   (defn sort-with-p [ai]
     "Return sorted ai and permutation array. Each entry of a must have the same
     type."
@@ -651,6 +628,28 @@
 
   (defn antidiag [v]
     (get (npy.diag (get (npy.array v) (s-all-rev))) (s-all-rev)))
+
+  (defn get-row-2norms [vs]
+    (assert (and (= (len vs.shape) 2)))
+    (sv nrm (second vs.shape)
+        den 0)
+    (for [i (range nrm)]
+      (+= den (** (getc vs i) 2)))
+    (sv den (npy.sqrt den))
+    den)
+
+  (defn scale-rows! [scale vs]
+    (assert (and (< (len scale.shape) 2)
+                 (= (len scale) (len vs))
+                 (= (len vs.shape) 2)))
+    (sv dim (second vs.shape))
+    (for [i (range dim)]
+      (sv (getc vs i) (* scale (getc vs i))))
+    vs)
+
+  (defn normalize-rows! [vs]
+    (scale-rows! (/ (get-row-2norms vs)) vs))
+
   )
  (except []
    (do
@@ -663,9 +662,9 @@
  (do
   (import matplotlib [matplotlib.pyplot :as pl])
 
-  (defn my-grid [&optional ls zorder]
-    (svifn ls "-" zorder -1)
-    (pl.grid True :lw 0.5 :ls ls :color (, 0.8 0.8 0.8) :zorder zorder
+  (defn my-grid [&optional ls]
+    (svifn ls "-")
+    (pl.grid True :lw 0.5 :ls ls :color (, 0.8 0.8 0.8) :zorder -1
              :which "both")
     (.set_axisbelow (pl.gca) True))
 
@@ -753,118 +752,26 @@
 
   (defn make-reference-slope-triangle [x-span y-span slope pattern
                                        &optional [kwargs-plot None]
-                                       [kwargs-text None]
-                                       [opposite None]]
+                                       [kwargs-text None]]
     (assert (= 2 (len x-span)))
     (assert (= 2 (len y-span)))
     (svifn kwargs-plot {})
     (svifn kwargs-text {})
-    (svifn opposite False)
     (sv (, x0 x1) x-span (, y0 y1) y-span
-        dx (- x1 x0) dy (- y1 y0))
-    (if opposite
-        (do (sv x [x0 x1 x1 x0]
-                y [y1 y1 y0 y1]))
-        (do (sv x [x0 x1 x0 x0]
-                y [y0 y0 y1 y0])))
+        dx (- x1 x0) dy (- y1 y0)
+        x [x0 x1 x0 x0] y [y0 y0 y1 y0])
     (pl.plot #*[x y pattern] #**kwargs-plot)
-    (when opposite
-      (if (none? kwargs-text) (sv kwargs-text {}))
-      (assoc kwargs-text "horizontalalignment" "right" "verticalalignment" "top"))
-    (pl.text #*[(if opposite (- x1 (* 0.1 dx)) (+ x0 (* 0.1 dx)))
-                (if opposite (- y1 (* 0.1 dy)) (+ y0 (* 0.1 dy))) (str slope)]
+    (pl.text #*[(+ x0 (* 0.1 dx)) (+ y0 (* 0.1 dy)) (str slope)]
              #**kwargs-text))
 
-  ) (except [] ))
+  (defn frame-init []
+    (pl.show :block False))
+  (defn frame-flip [U &optional [clim None] [pause 0.05]]
+    (pl.clf)
+    (pl.imshow U :interpolation "none" :origin "lower")
+    (unless (none? clim)
+      (pl.clim clim))
+    (pl.tight-layout)
+    (pl.pause pause))
 
-;;; More extensive unit tests.
-(if-main
-  (expect
-    (case/eq 'hi
-             ('hit 'nope)
-             ('hi 'bark 'yep))
-    'yep)
-  (expect
-    (case/in 'hello
-             ['(bye hello) 'nope]
-             ('(hi) 'yep)
-             ('(4) 'another-nope))
-    'nope)
-  (expect
-    (case/in 'hi
-             ('(bye hello) 'nope)
-             ('(hi) 'yep)
-             ('(4) 'another-nope))
-    'yep)
-  (expect
-    (case/in 'hi
-             ('(bye hello) 'nope)
-             ('hi 'yep)
-             ('(4) 'another-nope))
-    'yep)
-  (expect
-    (case/test in 'hi
-               ('(bye hello) 'nope)
-               ('hi 'yep)
-               ('(4) 'another-nope))
-    'yep)
-  (expect (case/eq 4 [5 'bye] [4 'hi])
-          'hi)
-  (expect
-    (case/in 'hi
-             ('(bye hello) 'nope))
-    None)
-  (expect
-    (case/in 'hi
-             ('(bye hello) 'nope)
-             (:else 'woot))
-    'woot)
-  (expect
-    (sdo (setv key 'hi)
-         (case/in key
-                  ('(bye hello) 'nope)
-                  ([key] 'yup)))
-    'yup)
-
-  (expect
-    (do
-      (sv a None b 4 c "hi" d None)
-      (svifn a "hi" b 3 c "bye")
-      (and (= a "hi") (= b 4) (= c "hi") (none? d))))
-  (expect
-    (do
-      (sv a None)
-      (svifn a "hi")
-      (= a "hi")))
-    (expect
-    (do
-      (sv a 5)
-      (svifn a "hi")
-      (= a 5)))
-
-    (expect
-     (do (sv a 1 b "hi" c 'd)
-         (svb (ab 1) (bb "hi") (cb 'd))
-         (and (= a ab) (= b bb) (= c cb))))
-
-  (when-inp ["test-pretty-print"]
-    (ppme '(case/in 'hi
-                    (['hi 'bye] 'hi)
-                    ('[foo] (for [d dinos]
-                              (print (.format "{:s} goes {:d}"
-                                              d.name d.sound))))
-                    ('[bar]
-                     (defn axis-tight-pad [&optional [pad 0.05] [mult False]
-                                           [foo 3]]
-                       (pl.axis "tight")
-                       (setv xl (pl.xlim) yl (pl.ylim))
-                       (pl.xlim (pad-lim xl pad mult))
-                       (pl.ylim (pad-lim yl pad mult)))
-                     (axis-tight-pad)))))
-
-  (when-inp ["ooa-from-file" {:fname str :anchor str :fieldno int}]
-    (sv (, errs ooas) (ooa-from-file fname fieldno :anchor anchor))
-    (for [i (range (len errs))]
-      (prf "{:10.3e} {}" (get errs i)
-           (if (zero? i) " n/a" (.format "{:6.3f}" (get ooas (dec i)))))))
-)
+  ) (except [] ))
diff --git a/methods/islet/figures/figs-adv-diag.hy b/methods/islet/figures/figs-adv-diag.hy
index 02af0b4..9ceab27 100644
--- a/methods/islet/figures/figs-adv-diag.hy
+++ b/methods/islet/figures/figs-adv-diag.hy
@@ -1,11 +1,13 @@
 (require [amb3 [*]])
 (import amb3 [amb3 [*]]
         [figsutils [*]]
-        math glob struct)
+        math glob struct os)
 
 (assoc matplotlib.rcParams "savefig.dpi" 300)
 (do (pl-require-type1-fonts))
 
+(sv gmd2 True)
+
 ;;; parse cmd, L, M, C slmmir output
 
 (defn acc-parse [fname &optional map-nstepfac]
@@ -119,7 +121,7 @@
   (defn praccv [pne ne vp v]
     (if (or (none? vp) (zero? vp) (zero? v))
         (prfno " {:7.1e} ( ----)" v)
-        (prfno " {:7.1e} ({:5.2f})" v (first (ooa [vp v] :x [ne pne])))))
+        (prfno " {:7.1e} ({:5.2f})" v (first (calc-ooa [vp v] :x [ne pne])))))
   (defn pracc [diagnostic norm pne pe ne e]
     (sv v (geton e diagnostic norm))
     (when (none? v) (return))
@@ -254,27 +256,28 @@
                                      (:measure o))
                       " relative error")
                    :fontsize (:fs o))
+        (pl.xlabel "Dynamics grid resolution" :fontsize (:fs o))
         (cond [(= (:yform o) :log2)
                (pl.yticks (list (range -40 10)) :fontsize (:fs o))])
         (when (and (:ref-cos-033 o) (= (:ic o) "cos") (= (:ode o) "nondivergent"))
-          (pl-plot xtform (* 0.033 (npy.ones (len xtform))) "k-."
+          (pl-plot xtform (* 0.033 (npy.ones (len xtform))) "-."
                    :zorder -10 :lw (:lw o) :color gray)))
       (when (:ooa-text o)
         (sv i (- (len x) 2))
-        (pl.text (- (get xtform (inc i)) 0.1) (* (get ytform (inc i)) 2)
+        (pl.text (+ (get xtform (inc i)) -0.08) (* (get ytform (inc i)) 2)
                  (.format "{:1.1f}"
-                          (- (first (ooa (cut y i (+ i 2))
-                                         :x (cut x i (+ i 2))))))))
+                          (first (calc-ooa (cut y i (+ i 2))
+                                           :x (cut x i (+ i 2)))))))
       (when (= np (last (:nps o)))
+        (when (:ref-ooa-2 o)
+          (sv ytref (y-tfn (* 0.7 (first y-ext) (** (/ (last x) (npa x)) 2))))
+          (pl-plot xtform ytref ":" :color gray))
         (when (= (:yform o) :semilogy)
           (sv (, y ys) (pl.yticks))
           (when (> (/ (last y) (first y)) 1e3)
             (for [i (range (len ys))]
               (sv (get ys i) (.format "{}" (int (math.log10 (get y i))))))
-            (pl.yticks y ys)))
-        (when (:ref-ooa-2 o)
-          (sv ytref (y-tfn (* 0.7 (first y-ext) (** (/ (last x) (npa x)) 2))))
-          (pl-plot xtform ytref "k:" :color gray))))
+            (pl.yticks y ys)))))
     y-ext)
 
   (defn legend [me ax entries &optional [o None] [nps-legend True] bbox]
@@ -321,9 +324,15 @@
   (+ main "\n"
      (flow-short2long (:ode o)) ", "
      (ic-short2long (cut (:ic o) 0 3)) ", "
-     (nstepfac->word (:nstepfac o)) " steps,\n"
-     (if (= (:prefine o) 5) "$p$-refinement, " "")
-     (+ (if (= (:cdrglb o) "none") "no " "") "property preservation")
+     (nstepfac->word (:nstepfac o)) " steps"
+     (if gmd2
+         ""
+         (+ ",\n"
+            (if (= (:prefine o) 5) "$p$-refinement, " "")
+            (+ (if (= (:cdrglb o) "none") "no " "") "property preservation")))
+     (if (and gmd2 (= (:cdrglb o) "none"))
+         ",\nno property preservation"
+         "")
      extra))
 
 (defn fig-stab-cmp [c d]
@@ -343,7 +352,7 @@
     (my-grid)
     (p.legend ax (, (, "k-" "Islet 1 cycle") (, "k--" "Islet 100 cycles")
                     (, "k-." "Natural 1 cycle") (, "k:" "OOA 2")) :o o)
-    (p.title (make-title "Islet stability:" o) o)))
+    (p.title (make-title "Islet method stability:" o) o)))
 
 (defn nextpow10 [f] (** 10.0 (math.ceil  (math.log10 f))))
 (defn prevpow10 [f] (** 10.0 (math.floor (math.log10 f))))
@@ -354,7 +363,8 @@
     general-prefine 5 show-linf True pp True)
   (sv p (AccFig)
       o (p.get-defaults c))
-  (defn plot [o title plot-fn &optional [ref-ooa-2 False]]
+  (defn plot [o title plot-fn &optional [ref-ooa-2 False] xlabel]
+    (svifn xlabel "Dynamics grid resolution")
     (sv fname (+ prefix "acc-" (:ode o) "-" (:ic o)
                  "-" (:timeint o) "-" (if (= (:cdrglb o) "none")
                                           "nopp" "pp")
@@ -368,8 +378,9 @@
       (when ref-ooa-2 (.append legs (, "k:" "OOA 2")))
       (when legend (p.legend ax legs :o o))
       (pl.ylabel "$\log_{10}$ relative error")
+      (pl.xlabel xlabel)
       (p.title (make-title title o) o)))
-  (assoc o :ode "nondivergent" :ic "gau" :nstepfac 1
+  (assoc o :ode "nondivergent" :ic "gau" :nstepfac 5
          :yform :semilogy :cdrglb "none" :cdrlcl "none"
          :timeint "exact" :prefine 0 :filter-floor 1e-11)
   (plot o "Islet empirical order of accuracy:"
@@ -380,7 +391,8 @@
             (p.plot ax d o)
             (pl.ylim (, 1e-11 1))
             (sv e (npy.array [0 -2 -4 -6 -8 -10]))
-            (pl.yticks (** 10.0 e) e))))
+            (pl.yticks (** 10.0 e) e)))
+        :xlabel "Reference resolution")
   (for [nstepfac (, 1 5)
         ic (, "gau" "cos" "slo")
         ode (, "nondivergent" "divergent" "rotate")]
@@ -388,7 +400,7 @@
            :cdrglb (if pp "caas-node" "none") :cdrlcl (if pp "caas" "none")
            :ic ic :filter-floor None :timeint general-timeint)
     (dont when (and (= ode "divergent") (= ic "gau")) (continue))
-    (plot o "Islet accuracy:"
+    (plot o "Islet method accuracy:"
           (fn [ax d o]
             (sv ye [1e10 -1e10])
             (defn update-ye [ye1]
@@ -418,7 +430,7 @@
   (assoc o :ic "gau" :ode "nondivergent" :cdrglb "none" :cdrlcl "none"
          :measure :l2 :yform :semilogy :nstepfac 1)
   (sv nps-str (make-nps-string (:nps o)))
-  (with [(pl-plot (:figsize o) (+ c.fig-dir "midpoint-check"))]
+  (with [(pl-plot (, 4 4.5) (+ c.fig-dir "midpoint-check"))]
     (sv ax (pl.subplot 1 1 1))
     (assoc o :pat-line "-" :C-line :C :ooa-text True :prefine 0 :timeint "exact" :nps [4])
     (p.plot ax d o)
@@ -430,10 +442,13 @@
     (p.plot ax d o)
     (assoc o :prefine 5 :timeint "interp" :nps c.nps)
     (p.plot ax d o)
+    (print (pl.xlim))
     (pl.ylim (, 5e-10 1))
     (my-grid)
     (p.legend ax (, (, "k-" "1 cycle") (, "k--" "1/2 cycle")) :o o)
-    (p.title (make-title "Trajectory interpolation:" o) o)))
+    (p.title (make-title "Trajectory interpolation:" o) o)
+    (sv xl (pl.xlim))
+    (pl.xlim (, (first xl) (* 1.02 (second xl))))))
 
 (defn fig-acc-mimic-src-term-midpoint [c d &optional [np-minus-2 False]]
   (sv p (AccFig)
@@ -454,6 +469,7 @@
       (when ref-ooa-2 (.append legs (, "k:" "OOA 2")))
       (p.legend ax legs :o o)
       (pl.ylabel "$\log_{10}$ $l_2$ relative error")
+      (pl.xlabel "Dynamics grid resolution")
       (p.title title o)))
   (for [nstepfac (, 5)
         ic (, "gau")
@@ -779,7 +795,7 @@
                  "nondivergent flow, 1$^\circ$, 576 time steps/cycle\n"
                  "$p$-refinement, property preservation")
               :fontsize (:fs o))
-    (p.legend ax (, (, "k-" "$l_2$") (, "k:" "$l_\infty$"))
+    (p.legend ax (, (, "k-" "$c_2$") (, "k:" "$c_\infty$"))
               :o o :bbox (, 0 0.8))
     (pl.ylim (, (** 10.0 -16) (** 10.0 -3)))
     (pl.xlim (, -0.2 10.2))
@@ -922,6 +938,171 @@
                     "Nondivergent flow")
                  :ha "center" :fontsize (:fs o))))))
 
+;;; images of instability
+
+(defn img-instab [c fname-meta-pairs outname]
+  (sv ncol 2 nrow 2 npanel (* ncol nrow) img-idx 7
+      w (/ 1 ncol) h (/ 1 nrow)
+      vmin -0.05 vmax 1.15)
+  (with [(pl-plot (, (* 3.2 ncol) (+ (* 2.1 nrow) 1)) outname :tight False)]
+    (sv axs [])
+    (for [spi (range npanel)]
+      (sv row (// spi ncol) col (% spi ncol))
+      (.append axs (pl.axes [(* w col) (- 1 (* h row)) (* 0.95 w) h]))
+      (sv mp (nth fname-meta-pairs spi)
+          img (nth (read-slmmir-io-arrays (first mp)) img-idx)
+          n (second img.shape)
+          img (get img (, (s-all)
+                          (+ (list (range (// n 2) n))
+                             (list (range 0 (// n 2))))))
+          img (get img (, (slice 60 450) (slice 120 750))))
+      (print "img range" (npy.min img) (npy.max img))
+      (draw-slmmir-image img :switch-halves False :grid False
+                         :vmin vmin :vmax vmax)
+      (pl.text 20 (- (first img.shape) 50)
+               (second mp)
+               :fontsize 15 :bbox {"facecolor" "white"}))
+    (when (= (inc spi) npanel)
+      (sv ticks (npy.linspace 0 1.1 12)
+          c (pl.colorbar :ax axs :orientation "horizontal" :aspect 20
+                         :shrink 0.8 :pad 0.05 :anchor (, 0.5 1.3)
+                         :ticks ticks))
+      (c.ax.tick-params :labelsize 14))))
+
+;;; stability grinding fig
+
+(defn file-exists? [fname] (os.path.exists fname))
+
+(defn make-out-fname [ode nonuni ne np timeint prefine cdrglb cdrlcl basis
+                      nstepfac ncycle]
+  (+ "ode_" ode "-nonuni_" (str (if nonuni 1 0)) "-ne_" (str ne) "-np_" (str np)
+     "-timeint_" timeint "-prefine_" (str prefine)
+     "-cdrglb_" cdrglb "-cdrlcl_" cdrlcl "-basis_" basis
+     "-nstepfac_" (str nstepfac) "-ncycle_" (str ncycle) ".out"))
+
+(defn stabgrind-make-spec [batchdir ode basis np nstepfac nes
+                           &optional timeint prefine cdrs nonuni ncycle]
+  (svifn timeint "exact" prefine 0 cdrs (, "none" "none") nonuni False ncycle 100)
+  (sv s (box-slots 'batchdir 'ode 'basis 'np 'nstepfac 'nes 'timeint 'prefine
+                   'cdrs 'nonuni 'ncycle)
+      s.type 'stabgrind-spec)
+  s)
+
+(defn stabgrind-keys [s]
+  (, s.timeint s.ode s.nstepfac s.basis (first s.cdrs) (second s.cdrs)
+     s.prefine s.np))
+
+(defn stabgrind-parse [s ne &optional d]
+  (svifn d {})
+  (sv fname (->> (make-out-fname s.ode s.nonuni ne s.np s.timeint s.prefine
+                                 (first s.cdrs) (second s.cdrs) s.basis s.nstepfac
+                                 s.ncycle)
+                 (+ s.batchdir "/")))
+  (unless (file-exists? fname)
+    (prf "Can't read {}" fname)
+    (return d))
+  (sv txt (.split (readall fname) "\n")
+      keys (stabgrind-keys s))
+  (for [ln txt]
+    (unless (= (cut ln 0 2) "C ") (continue))
+    (cond [(in "cycle" ln) (sv (, - - cyc) (sscanf ln "s,s,i"))]
+          [(or (in "PASS" ln) (in "FAIL" ln))
+           (sv (, - pf) (sscanf ln "s,s"))
+           (when (and (= pf "FAIL") (!= (first s.cdrs) "none") (<= cyc 10))
+             (prf "FAIL {}" cmd))]
+          [:else
+           (sv cl (parse-C ln))
+           (assoc-nested-append d (+ keys (, (:ic cl) ne)) cl)]))
+  d)
+
+(defn stabgrind-read [s]
+  (assert (= s.type 'stabgrind-spec))
+  (sv d (Box)
+      d.type 'stabgrind-read
+      d.s s
+      d.data {})
+  (for [ne s.nes]
+    (sv d.data (stabgrind-parse s ne :d d.data)))
+  d)
+
+(defn stabgrind-fig [c ds outname &optional ic]
+  (svifn ic "gau")
+  (sv nes (, 5 10 20 40 60 80)
+      yl (, 1e-11 1)
+      fsz (, 8 4) sbs (, 1 2)
+      ;;fsz (, 4 8) sbs (, 2 1)
+      fs 12
+      np (. (first ds) s np))
+  (defn xticks-ne []
+    (sv xa (nes->degstrs nes))
+    (pl.xticks (npy.log2 (:ne xa)) (:deg xa) :fontsize fs))
+  (defn xticks-cnt []
+    (sv xa (npy.array (list (, 1 10 100))))
+    (pl.xticks xa xa :fontsize fs))
+  (defn yticks []
+    (sv ya (npy.array (list (range -13 1))))
+    (pl.yticks (** 10.0 ya) ya :fontsize fs))
+  (defn make-pat [basis nstepfac &optional [marker False]]
+    (+ (get {1 "k" 3 "g" 5 "r"} nstepfac)
+       (if marker
+         (get {1 "o" 3 "v" 5 "s"} nstepfac)
+         "")
+       (get {"Gll" ":" "GllNodal" "-"} basis)))
+  (for [d ds]
+    (assert (= d.type 'stabgrind-read))
+    (sv s d.s keys (stabgrind-keys s)
+        x [] cls [])
+    (for [ne nes]
+      (sv cl (geton d.data #*(+ keys (, ic ne))))
+      (when (none? cl) (continue))
+      (.append x ne)
+      (.append cls cl))
+    (sv d.nes x d.cls cls))
+  (with [(pl-plot fsz outname)]
+    (pl.subplot #*sbs 1)
+    (for [d ds]
+      (sv s d.s keys (stabgrind-keys s)
+          ncyc (len (get d.data #*(+ keys (, ic (first d.nes))))))
+      (for [cyc (, 1 10 100)]
+        (when (> cyc ncyc) (break))
+        (sv y [])
+        (for [cl d.cls]
+          (.append y (:l2 (nth cl (dec cyc)))))
+        (prf "{} {:1d} {:3d} {:1.4e}" s.basis s.nstepfac cyc (last y))
+        (pl.semilogy (npy.log2 d.nes) y
+                     (make-pat s.basis s.nstepfac :marker True)
+                     :fillstyle "none")))
+    (do
+     (sv xl (pl.xlim))
+     (for [basis (, "Gll" "GllNodal") nstepfac (, 1 5)]
+       (pl.plot 0 10 (make-pat basis nstepfac)
+                :label (+ (get {"Gll" "GLL" "GllNodal" "Islet"} basis)
+                          " $n_p$ " (str np) " "
+                          (get {1 "long" 5 "short"} nstepfac))))
+     (pl.legend :loc "upper right" :fontsize fs :framealpha 0)
+     (pl.xlim xl))
+    (sv letter-y (* (second yl) 0.5)
+        f (pl.gcf))
+    (f.text 0.02 0.05 "(a)" :fontsize fs)
+    (xticks-ne) (yticks) (pl.ylim yl) (my-grid)
+    (pl.xlabel "Reference resolution" :fontsize fs)
+    (pl.ylabel "log$_{10}$ $l_2$ relative error" :fontsize fs)
+    (pl.title "Error vs. resolution at 1, 10, and 100 cycles"
+              :fontsize fs)
+    (pl.subplot #*sbs 2)
+    (for [d ds]
+      (sv s d.s)
+      (for [cls d.cls]
+        (sv y [])
+        (for [cl cls] (.append y (:l2 cl)))
+        (pl.loglog (range 1 (inc (len y))) y (make-pat s.basis s.nstepfac))))
+    (f.text 0.52 0.05 "(b)" :fontsize fs)
+    (pl.xlabel "Cycle" :fontsize fs)
+    (pl.ylabel "log$_{10}$ $l_2$ relative error" :fontsize fs)
+    (pl.title "Error vs. cycle at resolutions in (a)"
+              :fontsize fs)
+    (xticks-cnt) (yticks) (pl.ylim yl) (my-grid)))
+
 ;;; miscellaneous figs likely not to go in the paper
 
 (defn img-slo-cyl-tracer-grid [c direc outname &optional [ne 10]]
@@ -937,8 +1118,9 @@
                (, 16 "interp" "caas-node" 5)
                (, 16 "interp" "caas-node" 1))]
       (.append nps (first c))
-      (.append imgss (read-slmmir-io-arrays (make-fname direc ne (first c)
-                                              nstepfac (second c) (nth c 2) (last c)))))
+      (.append imgss (read-slmmir-io-arrays
+                      (make-fname direc ne (first c) nstepfac (second c)
+                                  (nth c 2) (last c)))))
     (sv spi 0)
     (for [idx (, 1 3 5)
           (, i imgs) (enumerate imgss)]
@@ -1033,3 +1215,48 @@
   (sv c (get-context)
       d (parse-footprint (+ c.data-dir fname)))
   (fig-comm-footprint c d))
+
+(when-inp ["img-filament-tracer-grid" {:direc str}]
+  (sv c (get-context))
+  (for [ne (, 5 10)]
+    (img-slo-cyl-tracer-grid
+      c direc (+ c.fig-dir (.format "img-filament-tracer-grid-ne{}" ne)) :ne ne)))
+
+(when-inp ["figs-double" {:fname str}]
+  (sv c (get-context)
+      c.nps (, 4 6 7 8 9 10 11 12 13)
+      d (acc-parse (+ c.data-dir fname)
+                   :map-nstepfac (if 0
+                                     (fn [n] (get {0.5 1 1 1 2.5 5 5 5} n))
+                                     (fn [n] (if (<= n 1) 1 5)))))
+  (figs-acc c d :prefix "ne2x-dt1x-" :ref-ooa-2 False :legend False
+            :general-timeint "exact" :general-prefine 0 :show-linf False
+            :pp True))
+
+(when-inp ["fig-instab"]
+  (sv c (get-context)
+      pairs [["none-cyc1-Gll" "GLL"]
+             ["caas-node-cyc1-Gll" "GLL-CAAS-node"]
+             ["none-cyc1-GllNodal" "Islet"]
+             ["caas-node-cyc1-GllNodal" "Islet-CAAS-node"]])
+  (for [p pairs]
+    (sv (get p 0) (+ c.data-dir "instab-imgs/instab-nondiv-ne20np6-"
+                     (first p) ".bin")))
+  (img-instab c pairs (+ c.fig-dir "img-instab")))
+
+(when-inp ["fig-stabgrind"]
+  (sv c (get-context)
+      ode "nondivergent"
+      np 11
+      nes (, 5 10 20 40 80)
+      ds [])
+  (for [basis (, "Gll" "GllNodal")
+        nstepfac (, 1 5)
+        ]
+    (sv s (stabgrind-make-spec
+           (+ "data/mar21/batch1" (if (= basis "Gll") "6" "5"))
+           ode basis np nstepfac nes
+           :ncycle (if (= basis "Gll") 20 100))
+        d (stabgrind-read s))
+    (unless (none? d) (.append ds d)))
+  (stabgrind-fig c ds (+ c.fig-dir "fig-stabgrind")))
diff --git a/methods/islet/figures/figs-methods.hy b/methods/islet/figures/figs-methods.hy
index cbd9107..b6753d5 100644
--- a/methods/islet/figures/figs-methods.hy
+++ b/methods/islet/figures/figs-methods.hy
@@ -1,7 +1,9 @@
 (require [amb3 [*]])
-(import amb3 [amb3 [*]] [figsutils :as futils]
+(import amb3 [amb3 [*]]
+        [figsutils :as futils]
+        [matplotlib :as mpl]
         [scipy.linalg :as linalg]
-        math re poly)
+        math re poly islet)
 
 (assoc matplotlib.rcParams "savefig.dpi" 300)
 (do (pl-require-type1-fonts))
@@ -64,7 +66,7 @@
   (with [f (open fname "w")]
     (defn write [s] (.write f s))
     (write "\\begin{center}\n")
-    (write "\\begin{tabular}{r|c|l|l}\n")
+    (write "\\begin{tabular}{rcll}\n")
     (write (+ (str.join " & " col-hdrs) " \\\\\n\\hline\n"))
     (when (= xnodes 'gll)
       (write "4 & 2 & see text & see text \\\\\n\\hline\n"))
@@ -85,18 +87,21 @@
 
 ;;; illustrations
 
+(defn draw-gll-dots [np marker &optional [dx 0] [dy 0] [markersize 12]
+                     fillstyle [alpha 1]]
+  (sv xgll (col-vec (poly.get-gll-x np))
+      e (npy.ones (, np 1))
+      X (npy.dot e (.transpose xgll))
+      Y (npy.dot xgll (.transpose e)))
+  (pl.plot (+ dx X) (+ dy Y) marker :markersize markersize :fillstyle fillstyle))
+
+(sv element-blue "#3838FF")
+
 (defn illustrate-grids [img-fname]
-  (defn draw-gll-dots [np marker]
-    (sv xgll (col-vec (poly.get-gll-x np))
-        e (npy.ones (, np 1))
-        X (npy.dot e (.transpose xgll))
-        Y (npy.dot xgll (.transpose e)))
-    (pl.plot X Y marker :markersize 12))
   (sv nf 6 elc "k")
   (with [(pl-plot (, 5 5) img-fname :format "pdf")]
     (pl.plot [-1 1] [1 1] elc [-1 1] [-1 -1] elc [1 1] [-1 1] elc [-1 -1] [-1 1] elc
-               :linewidth 2
-               :color (if 0 "#1D4BFA" "#3838FF"))
+             :linewidth 2 :color element-blue)
     (sv lw 2
         line "g--"
         d 0.99 nd (- d))
@@ -110,18 +115,17 @@
             (pl.plot [x x] [nd d] line :linewidth lw)
             (pl.plot [nd d] [x x] line :linewidth lw))))
     (draw-gll-dots 4 (+ elc "o"))
-    (draw-gll-dots 8 "r.")
+    (draw-gll-dots 8 "ro" :markersize 7)
     (pl.axis "equal")
     (pl.axis "off")))
 
 (defn draw-np4-schematic [img-fname]
-  (import islet)
   (sv np 4 fs 12
       x-gll (poly.get-gll-x np)
       nx 256
       x (npy.linspace -1 1 nx)
       isl (islet.Islet)
-      (, yn yo yb) (lfor m (, 0 1 3) (.transpose (isl.eval m np x)))
+      (, yn yo yb) (lfor m (, 0 2 1) (.transpose (isl.eval m np x)))
       clrs "kgrb"
       d 1.04)
   (with [(pl-plot (, 4 4) img-fname :tight True)]
@@ -160,7 +164,144 @@
     (pl.xlim (, (- d) d))
     (pl.xticks xt :fontsize fs)
     (pl.yticks (npy.linspace -0.2 1 7) :fontsize fs)
-    (pl.xlabel "Reference coordinate" :fontsize fs)))
+    (pl.xlabel "Reference-element coordinate" :fontsize fs)))
+
+(defn draw-islet-timestep-schematic [img-fname &optional [npv 4] [npt 6]]
+  (defn make-gll-grid [np]
+    (sv x (poly.get-gll-x np)
+        ones (npy.ones np)
+        X (npy.dot (col-vec x) (row-vec ones))
+        Y (npy.dot (col-vec ones) (row-vec x)))
+    (, X Y))
+
+  (defn rotate [X Y theta]
+    (sv c (npy.cos theta) s (npy.sin theta))
+    (, (+ (* (- s) Y) (* c X))
+       (+ (* c Y) (* s X))))
+  (defn scale [X Y p]
+    (sv rad (npy.sqrt (+ (** X 2) (** Y 2)))
+        s (+ 1 (* p (- rad 1))))
+    (, (* s X) (* s Y)))
+  (defn perturb [X Y p]
+    (sv np (first X.shape)
+        r (fn [] (as-> (npy.random.rand np np) x
+                       (* 2 x)
+                       (- x 1)
+                       (* p x))))
+    (, (+ X (r)) (+ Y (r))))
+  (defn get-v-pts [np]
+    (sv (, X0 Y0) (make-gll-grid np))
+    (+= X0 0.3)
+    (+= Y0 -0.2)
+    (sv (, X0 Y0) (as-> (, X0 Y0) xy
+                        (scale #*xy 0.3)
+                        (perturb #*xy 0.05) xy
+                        (rotate #*xy 0.2)))
+    (+= X0 0.95)
+    (+= Y0 -0.1)
+    (, X0 Y0))
+
+  (defn draw-v-elem [ax X Y npt v-marker v-marker-sz t-marker t-marker-sz]
+    (sv vgll (poly.get-gll-x npv)
+        f (fn [xi idxs]
+            (sv x (get X idxs) y (get Y idxs))
+            (-> (zip (poly.eval-lagrange-poly vgll x xi)
+                     (poly.eval-lagrange-poly vgll y xi))
+                (list)))
+        g (fn [x]
+            (+ (f x (, (s-all) -1))
+               (f x (, -1 (s-all-rev)))
+               (f x (, (s-all-rev) 0))
+               (f x (, 0 (s-all)))))
+        x (npy.linspace -1 1 51)
+        vtxs (g x))
+    (sv p (mpl.patches.Polygon vtxs :facecolor "w" :fill True :alpha 0.7
+                               :zorder 10))
+    (.add-patch ax p)
+    (sv vtxs (npy.array (unzip vtxs)))
+    (pl.plot (get vtxs 0) (get vtxs 1) "--" :color element-blue :zorder 11 :lw 3)
+    (pl.plot X Y v-marker :markersize 14 :zorder 12
+             :markerfacecolor "w" :markeredgecolor "k")
+    (sv tgll (poly.get-gll-x npt)
+        t-pts [])
+    (for [i (range npt) j (range npt)]
+      (sv vi (poly.eval-lagrange-poly-basis vgll (nth tgll i))
+          vj (poly.eval-lagrange-poly-basis vgll (nth tgll j))
+          basis (npy.dot (col-vec vi) (row-vec vj)))
+      (.append t-pts (, (npy.sum (* basis X)) (npy.sum (* basis Y)))))
+    (sv t-pts (npy.array (unzip t-pts))
+        t-in-mask (npy.logical-and (>= (get t-pts 0) 1)
+                                   (<= (get t-pts 1) 0))
+        t-out-mask (npy.logical-not t-in-mask))
+    (for [(, mask marker sz) (, (, t-in-mask "g*" 12)
+                                (, t-out-mask "go" 8))]
+      (pl.plot (get t-pts 0 mask) (get t-pts 1 mask) marker
+               :markersize sz :zorder 13)))
+
+  (defn plot-basis [np]
+    (sv isl (islet.Islet)
+        gll-best 1
+        pats {0 "--" gll-best "-"}
+        clrs "krbcgm"
+        x (npy.linspace -1 1 213)
+        fs 18)
+    (for [method (, 0 gll-best)]
+      (sv y (.transpose (isl.eval method np x))
+          pat (get pats method))
+      (for [i (range np)]
+        (pl.plot x (get y i) (+ (nth clrs (% i (len clrs))) pat)
+                 :label (if (zero? i)
+                            (get {0 "Natural GLL" gll-best "Islet GLL"}
+                                 method))))
+      (my-grid)
+      (sv d 1.04)
+      (pl.yticks (, 0 1) :fontsize fs)
+      (pl.xticks (, -1 0 1) :fontsize fs)
+      (pl.xlim (, (- d) d))
+      (pl.ylim (, -0.25 1.05)))
+    (pl.xlabel "Reference-element coordinate" :fontsize fs)
+    (pl.ylabel "Basis function value" :fontsize fs)
+    (pl.text -1.03 1.1 (.format "Basis functions, $n_p$ = {}" np) :fontsize fs)
+    (pl.figlegend :loc (, 0.42 0.36) :fontsize (- fs 1) :ncol 2))
+
+  (sv nelx 3 nely 2
+      f (fn [n] (- (/ n 2) 0.5))
+      hnelx (f nelx) hnely (f nely)
+      elw 2 f (fn [n] (* elw (- (/ n 2))))
+      x0 (f nelx) y0 (f nely)
+      (, X0 Y0) (get-v-pts npv)
+      one (npy.ones 2))
+
+  (with [(pl-plot (, 8 8) img-fname)]
+    (sv ax (pl.axes (, 0 0.35 1 0.65))
+        lw 3)
+    (for [i (range (inc nelx))]
+      (pl.plot (* (+ x0 (* i elw)) one) [y0 (+ y0 (* elw nely))] "-"
+               :color element-blue :lw lw))
+    (for [i (range (inc nely))]
+      (pl.plot [x0 (+ x0 (* elw nelx))] (* (+ y0 (* i elw)) one) "-"
+               :color element-blue :lw lw))
+    (for [i (range nelx) j (range (dec nely) -1 -1)]
+      (sv pats (if (and (= i 2) (= j 0))
+                   (, (, "ko" "full" 11) (, "ro" "full" 6))
+                   (, (, "ro" "full" 8))))
+      (for [pat pats]
+        (draw-gll-dots npt (first pat)
+                       :dx (* elw (- i hnelx)) :dy (* elw (- j hnely))
+                       :markersize (last pat) :fillstyle (second pat)
+                       :alpha 0.5)))
+    (draw-gll-dots npv "ko" :markersize 12 :dx -2 :dy 1)
+    (draw-v-elem ax X0 Y0 npt "ko" 12 "ro" 8)
+    (pl.axis "tight")
+    (pl.axis "equal")
+    (pl.axis "off")
+    (pl.text -3.045 -2.3 "(a)" :fontsize 18)
+    (sv x0 0.9 x1 3.1 x1m 2.9 y0 -2.1 y1 -1.9 y2 0.1)
+    (pl.plot [x0 x1 x1 x0 x0] [y0 y0 y1 y1 y0] "r--")
+    (pl.plot [x1m x1 x1 x1m x1m] [y0 y0 y2 y2 y0] "g--")
+    (pl.axes (, 0.1 0 0.85 0.3))
+    (plot-basis 6)
+    (pl.text -1.19 -0.5 "(b)" :fontsize 18)))
 
 ;;; utils for search
 
@@ -299,44 +440,68 @@
 (defn norm->str [n]
   (get {:l1 "$l_1$" :l2 "$l_2$" :li "$l_{\infty}$"} n))
 
+(defn plot1-slmmir-vs-heuristic [c d nps prop-preserve ic norm pum-thr lebesgue]
+  (sv npa npy.array fs 11
+      plot (if lebesgue pl.semilogy pl.loglog))
+  (for [np nps]
+    (sv es (get d #*(basic-key np :prop-preserve prop-preserve)))
+    (when (none? es) (continue))
+    (sv x [] y [])
+    (for [e es]
+      (sv cls (:cls e) b (:b e))
+      (when (> (:pum b) pum-thr) (continue))
+      (.append x (nth (get b (if lebesgue :lebesgue :npm))
+                      (get {:l1 0 :l2 1 :li 2} norm)))
+      (.append y (get cls ic norm))
+      (when (nodes= (:np b) (:nodes b) (:np b) (get-slmmir-builtin (:np b)))
+        (plot (last x) (last y) "ro" :fillstyle "none" :markersize 12
+              :zorder 20)))
+    (plot (npa x) (npa y) (+ (get c.npclrs np) (get c.npmarks np))
+          :fillstyle "none" :zorder (- 20 np) :label (.format "$n_p$ {}" np)))
+  (if lebesgue
+    (do (set-semilogy-ticks :fs fs)
+        (pl.xlabel (+ (norm->str norm) " heuristic") :fontsize fs))
+    (do (set-loglog-ticks :fs fs)
+        (pl.xlabel (+ "$\log_{10}$ " (norm->str norm) " heuristic")
+                   :fontsize fs)))
+  (my-grid)
+  (pl.title (.format (+ "Test problem vs.~heuristic:\n"
+                        "nondivergent flow, 1.5$^\circ$, {}, long steps\n"
+                        "$p$-refinement, {}")
+                     (futils.ic-short2long ic)
+                     (+ (if prop-preserve "" "no ") "property preservation"))
+            :fontsize fs)
+  (pl.legend :loc "best" :fontsize fs)
+  (pl.axis "tight")
+  (pl.ylabel (+ "$\log_{10}$ " (norm->str norm) " relative error")
+             :fontsize fs))
+
 (defn plot-slmmir-vs-heuristic [c d img-fname
-                                &optional nps prop-preserve ic norm pum-thr lebesgue]
+                                &optional nps prop-preserve ic norm pum-thr
+                                lebesgue]
   (svifn nps (list (range 4 14)) prop-preserve False ic "gau" norm :l2
          pum-thr 1e-6 lebesgue False)
   (sv npa npy.array fs 11
       plot (if lebesgue pl.semilogy pl.loglog))
   (print img-fname)
   (with [(pl-plot (, 4 4.2) img-fname :format "pdf")]
-    (for [np nps]
-      (sv es (get d #*(basic-key np :prop-preserve prop-preserve)))
-      (when (none? es) (continue))
-      (sv x [] y [])
-      (for [e es]
-        (sv cls (:cls e) b (:b e))
-        (when (> (:pum b) pum-thr) (continue))
-        (.append x (nth (get b (if lebesgue :lebesgue :npm))
-                        (get {:l1 0 :l2 1 :li 2} norm)))
-        (.append y (get cls ic norm))
-        (when (nodes= (:np b) (:nodes b) (:np b) (get-slmmir-builtin (:np b)))
-          (plot (last x) (last y) "ro" :fillstyle "none" :markersize 12
-                :zorder 20)))
-      (plot (npa x) (npa y) (+ (get c.npclrs np) (get c.npmarks np))
-            :fillstyle "none" :zorder (- 20 np) :label (.format "$n_p$ {}" np)))
-    (if lebesgue
-        (do (set-semilogy-ticks :fs fs)
-            (pl.xlabel (+ (norm->str norm) " heuristic") :fontsize fs))
-        (do (set-loglog-ticks :fs fs)
-            (pl.xlabel (+ "$\log_{10}$ " (norm->str norm) " heuristic") :fontsize fs)))
-    (my-grid)
-    (pl.title (.format (+ "Test problem vs.~heuristic:\n"
-                          "nondivergent flow, 1.5$^\circ$, {}, long steps\n"
-                          "$p$-refinement, {}")
-                       (futils.ic-short2long ic)
-                       (+ (if prop-preserve "" "no ") "property preservation"))
-              :fontsize fs)
-    (pl.legend :loc "best" :fontsize fs)
-    (pl.axis "tight")
-    (pl.ylabel (+ "$\log_{10}$ " (norm->str norm) " relative error") :fontsize fs)))
+    (plot1-slmmir-vs-heuristic c d nps prop-preserve ic norm pum-thr lebesgue)))
+
+(defn plot-slmmir-vs-heuristic-ab [c d img-fname norm-pp-ic-tuples
+                                   &optional nps pum-thr]
+  (svifn nps (list (range 4 14)) pum-thr 1e-6)
+  (sv lebesgue False)
+  (sv npa npy.array fs 11
+      plot pl.loglog)
+  (with [(pl-plot (, 7.8 4.2) img-fname :format "pdf")]
+    (for [i (range 2)]
+      (sv t (nth norm-pp-ic-tuples i))
+      (pl.subplot 1 2 (inc i))
+      (plot1-slmmir-vs-heuristic c d nps (nth t 1) (nth t 2) (nth t 0) pum-thr
+                                 lebesgue)
+      (sv f (pl.gcf))
+      (f.text (nth (, 0.02 0.52) i) 0.05 (nth (, "(a)" "(b)") i)
+              :fontsize fs))))
 
 ;;; pum vs perturb
 
@@ -356,6 +521,23 @@
                    (assoc-nested d (, basis np) (, perturb meam1))))]))
   d)
 
+(defn pum-make-reference-slope-triangle
+  [x-span y-span slope pattern
+   &optional [opposite False] [kwargs-plot None] [kwargs-text None]]
+    (assert (= 2 (len x-span)))
+    (assert (= 2 (len y-span)))
+    (svifn kwargs-plot {})
+    (svifn kwargs-text {})
+    (sv (, x0 x1) x-span (, y0 y1) y-span
+        dx (- x1 x0) dy (- y1 y0)
+        x [x0 x1 x1 x0]
+        y [y1 y1 y0 y1])
+    (pl.plot #*[x y pattern] #**kwargs-plot)
+    (pl.text #*[(* 0.7 x1)
+                (* 0.05 y0)
+                (str slope)]
+             #**kwargs-text))
+
 (defn plot-pum-vs-perturb [c d fname]
   (sv f identity fs 11)
   (with [(pl-plot (, 4 4) fname)]
@@ -368,18 +550,18 @@
                    (+ (get c.npclrs np) (if uni "--" "-"))
                    :label (+ "$n_p$ " (str np) (if uni "'" "")))))
     (sv f 1.7 slope 4)
-    (make-reference-slope-triangle [(/ 2.8e-3 f) (* 2.8e-2 f)]
-                                   [(* 3e-8 (** f slope)) (/ 3e-12 (** f slope))]
-                                   slope "k-"
-                                   :opposite True
-                                   :kwargs-text {"fontsize" (+ fs 4)})
+    (pum-make-reference-slope-triangle
+     [(/ 2.8e-3 f) (* 2.8e-2 f)]
+     [(* 3e-8 (** f slope)) (/ 3e-12 (** f slope))]
+     slope "k-"
+     :kwargs-text {"fontsize" (+ fs 4)})
     (set-loglog-ticks)
     (my-grid)
     (pl.legend :loc "upper left" :fontsize (dec fs) :ncol 1)
     (pl.xlabel "$\log_{10}$ Element size relative random perturbation $\delta$"
                :fontsize fs)
     (pl.ylabel "$\log_{10}$ (max $|\lambda|$ - 1)" :fontsize fs)
-    (pl.title "Perturbed uniform mesh metric" :fontsize fs)
+    (pl.title "Perturbed uniform grid metric" :fontsize fs)
     (pl.xlim (, 1e-4 0.1))
     (pl.ylim (, 1e-14 0.1))))
 
@@ -436,7 +618,7 @@
       clrs {(nth ms 0) "r" (nth ms 1) "g" (nth ms 2) "k"}
       mrks {(nth ms 0) "." (nth ms 1) "x" (nth ms 2) "o"}
       lbls {(nth ms 0) "GLL natural" (nth ms 1) "Uniform offset nodal subset"
-            (nth ms 2) "GLL nodal subset"})
+            (nth ms 2) "GLL offset nodal subset"})
   (with [(pl-plot (, 4 4) fname)]
     (for [method ms]
       (sv e (get dpum method))
@@ -453,7 +635,7 @@
     (set-semilogy-ticks)
     (pl.ylim (, 1e-16 1))
     (pl.ylabel "$\log_{10}$ (max $|\lambda|$ - 1)" :fontsize fs)
-    (pl.xlabel "Translation, $\Delta x$, relative to element size 1" :fontsize fs)
+    (pl.xlabel "Translation distance $\Delta x$ relative to element size 1" :fontsize fs)
     (pl.title (.format "$n_p$ {} methods" np) :fontsize fs)
     (pl.legend :loc (, 0.02 0.15) :fontsize fs)
     (pl.xticks (npy.linspace 0 1 11))
@@ -500,9 +682,17 @@
       :nps [6 7 8 9 10]
       :prop-preserve pp :norm norm :ic ic :lebesgue lebesgue)))
 
-(when-inp ["illustrations"]
-  (sv c (futils.get-context))
-  (illustrate-grids (+ c.fig-dir "illustrate-grids")))
+(when-inp ["plot-slmmir-vs-heuristic-ab"]
+  (sv fnames (, "slmmir-on-basis-lines-2.txt")
+      c (futils.get-context)
+      lebesgue False
+      d {})
+  (for [fname fnames]
+    (sv d (parse-slmmir-output (+ "data/" fname) :d d :lebesgue lebesgue)))
+  (plot-slmmir-vs-heuristic-ab
+   c d (+ c.fig-dir "slmmir-vs-heuristic-ab")
+   (, (, :l2 False "gau") (, :l2 True "cos"))
+   :nps [6 7 8 9 10]))
 
 (when-inp ["tables"]
   (for [xnodes (, 'gll)]
@@ -536,28 +726,31 @@
 (defn plot-basis-schematic [np &optional [annotate False]]
   (import islet)
   (sv c (futils.get-context)
-      pats {0 "--" 3 "-"}
+      gll-best 1
+      pats {0"--" gll-best "-"}
       clrs "krbcgm"
       x (npy.linspace -1 1 512)
       isl (islet.Islet)
       fs 12)
-  (with [(pl-plot (, 6 (if annotate 4 3))
+  (with [(pl-plot (, 5 (if annotate 3.33 2.5))
                   (+ c.fig-dir "basis-schematic-np" (str np)
                      (if annotate "-annotated" "")))]
-    (for [method (, 0 3)]
+    (for [method (, 0 1)]
       (sv y (.transpose (isl.eval method np x))
           pat (get pats method))
       (for [i (range np)]
         (pl.plot x (get y i) (+ (nth clrs (% i (len clrs))) pat)
-                 :label (if (zero? i) (get {0 "Natural GLL" 3 "Islet GLL"} method))))
+                 :label (if (zero? i)
+                            (get {0 "Natural GLL" gll-best "Islet GLL"}
+                                 method))))
       (my-grid)
       (sv d 1.04)
       (pl.xlim (, (- d) d))
       (pl.ylim (, (if annotate -0.64 -0.22) 1.03)))
-    (pl.xlabel "Reference coordinate" :fontsize fs)
+    (pl.xlabel "Reference-element coordinate" :fontsize fs)
     (pl.ylabel "Basis function value" :fontsize fs)
-    (pl.text -1 1.1 (.format "Basis functions, $n_p$ = {}" np) :fontsize fs)
-    (pl.figlegend :loc (, 0.49 (if annotate 0.91 0.88)) :fontsize fs :ncol 2)
+    (pl.text -1.22 1.1 (.format "Basis functions, $n_p$ = {}" np) :fontsize fs)
+    (pl.figlegend :loc (, 0.44 0.9) :fontsize (- fs 1) :ncol 2)
     (when annotate
       (sv npa npy.array
           xgll (npa (poly.get-gll-x np))
@@ -567,7 +760,7 @@
           clr "g" y -0.3 w 0.006 lw 2 ones (npa [1 1]))
       (for [i (range 2)]
         (pl.arrow (nth xgll (+ ireg i)) y 0 0.1 :width w :color clr))
-      (pl.plot xs (* y ones) (+ clr "-") (* xc ones) [-0.36 y] :color clr :lw lw)
+      (pl.plot xs (* y ones) "-" (* xc ones) [-0.36 y] :color clr :lw lw)
       (pl.text xc -0.45 (.format "Region {}" ireg)
                :color clr :ha "center" :fontsize fs)
       (when (= np 6)
@@ -581,4 +774,450 @@
                  :color clr :ha "center" :fontsize fs)))))
 
 (when-inp ["basis-schematic" {:np int}]
-  (plot-basis-schematic np :annotate True))
+  (plot-basis-schematic np :annotate (= np 6)))
+
+(when-inp ["illustrations"]
+  (sv c (futils.get-context))
+  (illustrate-grids (+ c.fig-dir "illustrate-grids"))
+  (draw-islet-timestep-schematic (+ c.fig-dir "islet-timestep-schematic")))
+
+(when-inp ["islet-timestep-schematic"]
+  (sv c (futils.get-context))
+  (draw-islet-timestep-schematic (+ c.fig-dir "islet-timestep-schematic")))
+
+(defn ms-draw-block [np x-beg y-beg ulen clr lw fs]
+  (sv n (dec np)
+      xl (* np ulen)
+      yl (* n ulen)
+      pat (+ clr "-"))
+  (pl.plot [x-beg (+ x-beg (* np ulen))]
+           [y-beg y-beg]
+           pat :lw lw)
+  (pl.plot [x-beg (+ x-beg (* np ulen))]
+           [(+ y-beg yl) (+ y-beg yl)]
+           pat :lw lw)
+  (pl.plot [x-beg x-beg]
+           [y-beg (+ y-beg yl)]
+           pat :lw lw)
+  (pl.plot [(+ x-beg xl) (+ x-beg xl)]
+           [y-beg (+ y-beg yl)]
+           pat :lw lw)
+  (sv x (+ x-beg (- xl ulen)))
+  (pl.plot [x x]
+           [y-beg (+ y-beg yl)]
+           (+ clr "--")  :lw lw)
+  (pl.text (+ x-beg (* 0.5 n ulen)) (+ y-beg (* 0.5 yl))
+           "$\\bar{B}$"
+           :fontsize fs :ha "center" :va "center" :color clr)
+  (pl.text (+ x-beg (* (+ n 0.5) ulen)) (+ y-beg (* 0.5 yl))
+           "$b$"
+           :fontsize fs :ha "center" :va "center" :color clr))
+
+(defn ms-write-numbers [np down x-beg y-beg dx dy fs]
+  (sv n (+ down (* 2 (dec np))))
+  (for [i (range n)]
+    (pl.text (+ x-beg (* i dx))
+             (+ y-beg (* i dy))
+             (str i)
+             :fontsize fs :ha "center" :va "center")))
+
+(defn ms-draw-elements [np x-beg y ulen line-pat node-pat fill-style]
+  (sv n (dec np)
+      e-len (* n ulen)
+      x-len (* ulen (* 2 n))
+      x-end (+ x-beg x-len)
+      e-bdys [x-beg (+ x-beg (* n ulen)) x-end]
+      xgll (poly.get-gll-x np))
+  (pl.plot [x-beg x-end] [y y] line-pat
+           [x-end (* (dec (* 2 np)) ulen)] [y y] (+ (first line-pat) ":")
+           e-bdys [y y y]
+           :lw 4 node-pat :fillstyle fill-style :markersize 16 :mew 4)
+  (sv xs e-bdys)
+  (for [i (range 2)]
+    (sv xphys (cut (+ x-beg (* i e-len) (* (/ (+ 1 (npy.array xgll)) 2) e-len))
+                   1 -1))
+    (.extend xs (list xphys))
+    (pl.plot xphys (* y (npy.ones (len xphys)))
+             node-pat :fillstyle fill-style :markersize 11 :mew 4))
+  (.sort xs)
+  xs)
+
+(defn ms-annotate-elements [xs y ylbl lbl fs]
+  (for [i (range (len xs))]
+    (pl.text (nth xs i) y (str i) :fontsize fs :ha "center" :va "center"))
+  (pl.text (first xs) ylbl lbl :fontsize fs :ha "left" :va "center"))
+
+(defn make-matrix-schematic [&optional np down filename]
+  (sv c (futils.get-context))
+  (svifn np 4 down 1 filename (+ c.fig-dir "matrix-schematic"))
+  (assert (>= np 4))
+  (assert (>= down 0))
+
+  (sv n (dec np)
+      ulen 1 hulen (/ ulen 2)
+      L (* (+ down (* 2 (dec np))) ulen)
+      fs 20 fsb (int (* fs 2))
+      xgll (poly.get-gll-x np)
+      xlate (- (first xgll) (nth xgll down)))
+  
+  (with [(pl-plot (, 6 7) filename)]
+    (pl.plot [0 L] [L L] "k-"
+             [0 0] [0 L] "k-")
+    (ms-write-numbers np down hulen (+ L hulen) 1 0 fs)
+    (ms-write-numbers np down (- hulen) (- L hulen) 0 -1 fs)
+    (ms-draw-block np 0 (- L (* ulen (+ down n))) ulen "r" 2 fsb)
+    (ms-draw-block np n 0 ulen "r" 2 fsb)
+    (sv xs-src (ms-draw-elements np hulen (* -2 ulen) ulen
+                                 "b-" "bo" "full")
+        xs-tgt (ms-draw-elements np (+ hulen xlate) (* -1.35 ulen) ulen
+                                 "g--" "go" "none"))
+    (ms-annotate-elements xs-src (* -2.5 ulen) (* -3 ulen)
+                          "Source (Eulerian)" fs)
+    (ms-annotate-elements xs-tgt (* -0.85 ulen) (* -0.35 ulen)
+                          "Target (Lagrangian)" fs)
+    (pl.axis "off")))
+
+(when-inp ["matrix-schematic"]
+  (make-matrix-schematic))
+
+(defn get-unit-cube-corners []
+  (npy.array [[1 -1 -1] [1 1 -1] [-1 1 -1] [-1 -1 -1]
+              [1 -1 1] [1 1 1] [-1 1 1] [-1 -1 1]]))
+
+(defn normalize [v] (/ v (npy.linalg.norm v)))
+
+(defn cs-get-cubedsphere-ne2-lines [n]
+  (assert (>= n 2))
+  (sv npa npy.array
+      c (get-unit-cube-corners))
+  (defn avg [&rest args]
+    (sv p (.copy (nth c (first args))))
+    (for [i (range 1 (len args))]
+      (+= p (nth c (nth args i))))
+    (/ p (len args)))
+  (defn arc [a b]
+    (npy.concatenate (, (col-vec (npy.linspace (nth a 0) (nth b 0) n))
+                        (col-vec (npy.linspace (nth a 1) (nth b 1) n))
+                        (col-vec (npy.linspace (nth a 2) (nth b 2) n)))
+                     :axis 1))
+  (defn add [sq d]
+    (sv sqd (.copy sq))
+    (for [e sqd] (+= e d))
+    sqd)
+  (defn xform [v i is j js k ks]
+    (sv x [])
+    (for [e v]
+      (.append x (npy.concatenate
+                  (, (* is (col-vec (getc e i)))
+                     (* js (col-vec (getc e j)))
+                     (* ks (col-vec (getc e k))))
+                  :axis 1)))
+    x)
+  (sv sq (npy.concatenate (, (arc (nth c 0) (avg 0 1))
+                             (cut (arc (avg 0 1) (avg 0 1 4 5)) 1)
+                             (cut (arc (avg 0 1 4 5) (avg 0 4)) 1)
+                             (cut (arc (avg 0 4) (nth c 0)) 1))
+                          :axis 0))
+  (sv side [] v [])
+  (for [d (, (add sq [0 0 0])
+             (add sq [0 1 0])
+             (add sq [0 1 1])
+             (add sq [0 0 1]))]
+    (.append side d))
+  (.extend v (xform side 0 1 1 1 2 1))
+  (.extend v (xform side 1 -1 0 1 2 1))
+  (.extend v (add (xform side 0 1 1 -1 2 1) [-2 0 0]))
+  (.extend v (xform side 1 1 0 -1 2 1))
+  (.extend v (xform side 1 1 2 1 0 1))
+  (.extend v (add (xform side 1 1 2 1 0 1) [0 0 -2]))
+  (for [e v] (normalize-rows! e))
+  v)
+
+(defn cs-make-separate-lines [lines]
+  ;; make each arc a bunch of separate lines.
+  (sv sl [])
+  (for [e lines]
+    (for [i (range (dec (len e)))]
+      (.append sl (get e (, [i (inc i)] (s-all))))))
+  sl)
+
+(defn cs-project-lines [nml lines]
+  (sv vfront [] vback []
+      xhat (normalize (npy.cross [0 0 1] nml))
+      yhat (col-vec (normalize (npy.cross nml xhat)))
+      xhat (col-vec xhat)
+      nml (col-vec (normalize nml)))
+  (for [e lines]
+    (sv p (npy.concatenate (, (npy.dot e xhat)
+                              (npy.dot e yhat))
+                           :axis 1)
+        s (npy.dot e nml)
+        fmask (vectorize (>= s 0))
+        bmask (npy.logical-not fmask))
+    (.append vfront (get p (, fmask (s-all))))
+    (.append vback  (get p (, bmask (s-all)))))
+  (, vfront vback xhat yhat nml))
+
+(defn cs-ref-to-sphere [crnrs a b]
+  (sv oma (- 1 a)
+      omb (- 1 b))
+  (normalize (+ (* a b (nth crnrs 0))
+                (* oma b (nth crnrs 1))
+                (* oma omb (nth crnrs 2))
+                (* a omb (nth crnrs 3)))))
+
+(defn make-gll-pts [crnrs np]
+  (sv r (poly.get-gll-x np)
+      ps [])
+  (for [i (range np)]
+    (sv a (/ (+ 1 (nth r i)) 2))
+    (for [j (range np)]
+      (sv b (/ (+ 1 (nth r j)) 2))
+      (.append ps (cs-ref-to-sphere crnrs a b))))
+  ps)
+
+(defn cs-make-physgrid [crnrs nf n]
+  (sv lns []
+      zon (npy.linspace 0 1 n))
+ (for [i (range (inc nf))]
+    (sv a (/ i nf)
+        ln (npy.zeros (, n 3)))
+    (for [(, j b) (enumerate zon)]
+      (sv (get ln j) (cs-ref-to-sphere crnrs a b)))
+    (.append lns ln))
+ (for [i (range (inc nf))]
+    (sv b (/ i nf)
+        ln (npy.zeros (, n 3)))
+    (for [(, j a) (enumerate zon)]
+      (sv (get ln j) (cs-ref-to-sphere crnrs a b)))
+    (.append lns ln))
+  lns)
+
+(defn cubed-sphere-subelem-grid-schematic [&optional npv npt filename nf]
+  (sv c (futils.get-context))
+  (svifn npv 4 npt 6
+         filename (+ c.fig-dir "cubed-sphere-subelem-grid-schematic"))
+  (assert (and (>= npv 4) (>= npt npv)))
+
+  (sv vlines (cs-get-cubedsphere-ne2-lines 50)
+      vlines (cs-make-separate-lines vlines)
+      nml-ref [0.3 -0.7 0.2]
+      (, vfront vback xhat yhat nml) (cs-project-lines nml-ref vlines))
+
+  (defn project-pts [g]
+    (sv p (npy.zeros (, (len g) 2)))
+    (for [(, i e) (enumerate g)]
+      (sv (get p (, i (s-all))) [(npy.float64 (npy.dot e (col-vec xhat)))
+                                 (npy.float64 (npy.dot e (col-vec yhat)))]))
+    p)
+  (sv crnrs (lfor e [[1 -1 1] [0 -1 1] [0 -1 0] [1 -1 0]]
+                  (row-vec (normalize e)))
+      gv (project-pts (make-gll-pts crnrs npv))
+      gt (project-pts (make-gll-pts crnrs npt)))
+
+  (with [(pl-plot (, 6 6) filename)]
+    (sv ax (pl.gca)
+        lim (, -1.1 1.1)
+        g 0.8
+        lcb (mpl.collections.LineCollection vback  :color [g g g] :lw 2)
+        lcf (mpl.collections.LineCollection vfront :color [0 0 0] :lw 2)
+        theta (npy.linspace 0 (* 2 math.pi) 200))
+    (ax.add-collection lcb)
+    (sv g 0.4)
+    (pl.plot (npy.cos theta) (npy.sin theta) "--" :lw 2
+             :color (if (none? nf) "g" [g g g]))
+    (ax.add-collection lcf)
+    (unless (none? nf)
+      (sv pgls (cs-make-separate-lines (cs-make-physgrid crnrs nf 50))
+          (, pgls - - - -) (cs-project-lines nml-ref pgls)
+          lpg (mpl.collections.LineCollection pgls :color "g" :lw 1))
+      (ax.add-collection lpg))
+    (pl.plot (getc gv 0) (getc gv 1) "bo" :markersize 12)
+    (pl.plot (getc gt 0) (getc gt 1) "ro" :markersize 6)
+    (pl.xlim lim) (pl.ylim lim)
+    (pl.axis "square")
+    (pl.axis "off")))
+
+(when-inp ["cubed-sphere-subelem-grid-schematic"]
+  (cubed-sphere-subelem-grid-schematic :nf 2))
+
+(defn make-gll-grid [xb np &optional [end False]]
+  (sv ne (dec (len xb))
+      d (dec np)
+      n (* ne d)
+      xgll (npy.zeros (if end (inc n) n))
+      ref (* 0.5 (+ 1 (npy.array (cut (poly.get-gll-x np) 0 -1)))))
+  (for [ie (range ne)]
+    (sv dx (- (get xb (inc ie)) (get xb ie))
+        (get xgll (slice (* ie d) (* (inc ie) d))) (+ (get xb ie) (* dx ref))))
+  (when end (sv (get xgll n) (+ (- (last xb) (first xb)) (first xgll))))
+  xgll)
+
+(defn ones [x] (npy.ones x.shape))
+
+(defn wrap-01 [xi]
+  (sv mask (< xi 0)
+      (get xi mask) (+ (get xi mask) 1)
+      mask (> xi 1)
+      (get xi mask) (- (get xi mask) 1))
+  xi)
+
+(defn interp-sup [x y isup xi]
+  (sv xsup (npy.array (lfor i isup (get x i)))
+      ysup (npy.array (lfor i isup (get y i))))
+  (poly.eval-lagrange-poly xsup ysup xi))
+
+(defn csl-cubic-interp [x y xi]
+  (sv nc (dec (len x))
+      yi (npy.zeros xi.shape)
+      xi (wrap-01 xi)
+      nodes [])
+  (for [ic (range nc)]
+    (sv mask (npy.logical-and (>= xi (get x ic)) (< xi (get x (inc ic)))))
+    (.append nodes (list (first (npy.nonzero mask))))
+    (unless (npy.any mask) (continue))
+    (sv isup [(% (dec ic) nc) ic (inc ic) (% (+ ic 2) nc)]
+        (get yi mask) (interp-sup x y isup (get xi mask))))
+  (, yi nodes))
+
+(defn gll-interp [np x y xi]
+  (sv d (dec np)
+      nc (// (dec (len x)) d)
+      yi (npy.zeros xi.shape)
+      xi (wrap-01 xi)
+      nodes [])
+  (for [ic (range nc)]
+    (sv xlo (get x (* ic d))
+        xhi (get x (* (inc ic) d))
+        mask (npy.logical-and (>= xi xlo) (< xi xhi)))
+    (.append nodes (list (first (npy.nonzero mask))))
+    (unless (npy.any mask) (continue))
+    (sv isup (list (range (* ic d) (inc (* (inc ic) d))))
+        (get yi mask) (interp-sup x y isup (get xi mask))))
+  (, yi nodes))
+
+(defn isl-1d-schematic [&optional ne-gll ne-csl np filename]
+  (sv c (futils.get-context))
+  (svifn ne-gll 5 np 4
+         filename (+ c.fig-dir "isl-1d-schematic"))
+  (svifn ne-csl (* ne-gll (dec np)))
+  (assert (>= np 4))
+
+  (defn f-fn [x]
+    (sv ctr 0.4
+        fac 3
+        f (* 0.5 (+ 1 (npy.cos (* math.pi (* fac (- x ctr))))))
+        mask (npy.logical-or (< x (- ctr (/ fac))) (> x (+ ctr (/ fac))))
+        (get f mask) 0)
+    f)
+
+  (sv csl-grid (npy.linspace 0 1 (inc ne-csl))
+      el-grid (npy.linspace 0 1 (inc ne-gll))
+      gll-grid (make-gll-grid el-grid np :end True))
+
+  (defn draw-csl-grid [xos yos]
+    (pl.plot (+ xos (npy.array [0 1])) [yos yos] "k-")
+    (pl.plot (+ xos csl-grid) (* yos (ones csl-grid)) "ko"
+             :markersize 4 :fillstyle "full"))
+
+  (defn draw-gll-grid [xos yos]
+    (pl.plot (+ xos el-grid) (* yos (ones el-grid)) "ko"
+             :markersize 7 :fillstyle "full")
+    (pl.plot (+ xos (npy.array [0 1])) [yos yos] "k-")
+    (pl.plot (+ xos gll-grid) (* yos (ones gll-grid)) "ko"
+             :markersize 4 :fillstyle "full"))
+
+  (sv xlate -0.235
+      f-csl (f-fn csl-grid)
+      (, i-csl nodes-csl) (csl-cubic-interp csl-grid f-csl (+ csl-grid xlate))
+      f-gll (f-fn gll-grid)
+      (, i-gll nodes-gll) (gll-interp np gll-grid f-gll (+ gll-grid xlate)))
+  (print nodes-csl)
+
+  (with [(pl-plot (, 8 3.5) filename)]
+    (sv gap 0.1
+        yos 1.5
+        tfs 12
+        p 0.5
+        gray [p p p])
+
+    (pl.text 0.5 (+ yos 1.2) "Classical ISL"
+             :fontsize tfs :ha "center")
+    (pl.text (+ 1.5 gap) (+ yos 1.2) "Element-based ISL"
+             :fontsize tfs :ha "center")
+
+    (pl.text 0 (+ yos 0.9) "Time step n+1" :fontsize tfs)
+    (pl.text 0 0.9 "Time step n" :fontsize tfs)
+    (for [i (get nodes-csl 5)]
+      (pl.plot [(+ (get csl-grid i) xlate) (get csl-grid i)]
+               [0 yos] "g--")
+      (pl.plot (+ (get csl-grid i) xlate) 0 "gx")
+      (pl.plot (get csl-grid i) yos "gs" :fillstyle "none"))
+
+    (sv i 9
+        d -0.2)
+    (pl.plot (get csl-grid i) (+ yos (get i-csl i)) "ro"
+             :fillstyle "none")
+    (pl.plot [(+ (get csl-grid i) xlate) (get csl-grid i)]
+             [d d] "g-")
+    (pl.plot (get csl-grid i) d "gs" :fillstyle "none")
+    (pl.plot (+ (get csl-grid i) xlate) d "gx"
+             :fillstyle "none")
+    (pl.text (+ (get csl-grid i) 0.02) d
+             "Translation distance $\Delta x$"
+             :fontsize tfs :va "center")
+
+    (draw-csl-grid 0 0)
+    (sv xi (npy.linspace 0 1 111)
+        (, yi -) (csl-cubic-interp csl-grid f-csl xi))
+    (pl.plot xi yi "-" :color gray)
+    (sv isup (list (range 4 8))
+        xi (npy.linspace (get csl-grid (first isup))
+                         (get csl-grid (last isup))
+                         30)
+        yi (interp-sup csl-grid f-csl isup xi))
+    (pl.plot xi yi "r:" :lw 2)
+    (sv xi (npy.linspace (get csl-grid (get isup 1))
+                         (get csl-grid (get isup 2))
+                         10)
+        yi (interp-sup csl-grid f-csl isup xi))
+    (pl.plot xi yi "r-" :lw 2)
+    (pl.plot (get csl-grid isup) (get f-csl isup) "ro"
+             :fillstyle "none")
+    (pl.plot csl-grid f-csl "k.")
+    
+    (draw-csl-grid 0 yos)
+    (pl.plot csl-grid (+ yos i-csl) "k.")
+
+    (sv ni 1 dx (+ 1 gap))
+    (draw-gll-grid (+ dx) 0)
+    (draw-gll-grid (+ dx) yos)
+
+    (for [i (get nodes-gll ni)]
+      (pl.plot [(+ (get gll-grid i) xlate dx) (+ (get gll-grid i) dx)]
+               [0 yos] "g--")
+      (pl.plot (+ (get gll-grid i) xlate dx) 0 "gx")
+      (pl.plot (+ (get gll-grid i) dx) yos "gs" :fillstyle "none"))
+    
+    (sv xi (npy.linspace 0 1 111)
+        (, yi -) (gll-interp np gll-grid f-gll xi))
+    (pl.plot (+ dx xi) yi "-" :color gray)
+    (sv isup (list (range (* (dec np) ni) (inc (* (dec np) (inc ni)))))
+        xi (npy.linspace (get gll-grid (first isup))
+                         (get gll-grid (last isup))
+                         30)
+        yi (interp-sup gll-grid f-gll isup xi))
+    (pl.plot (+ dx xi) yi "b-" :lw 2)
+    (pl.plot (+ dx (get gll-grid isup)) (get f-gll isup) "bo"
+             :fillstyle "none")
+    (pl.plot (+ dx gll-grid) f-gll "k.")
+    
+    (pl.plot (+ dx gll-grid) (+ yos i-gll) "k.")
+    (sv idxs (get nodes-gll ni))
+    (pl.plot (+ dx (get gll-grid idxs)) (+ yos (get i-gll idxs)) "bo"
+             :fillstyle "none")
+    
+    (pl.axis "off")))
+
+(when-inp ["isl-1d-schematic"]
+  (isl-1d-schematic))
diff --git a/methods/islet/figures/figsutils.hy b/methods/islet/figures/figsutils.hy
index 69e830e..f359731 100644
--- a/methods/islet/figures/figsutils.hy
+++ b/methods/islet/figures/figsutils.hy
@@ -14,7 +14,7 @@
       c.nes (, 5 10 20 40 80)
       c.nps (, 4 6 8 9 12) ;(, 4 5 6 7 8 9 10 11 12 13 16)
       c.ics (, "gau" "cos" "slo")
-      c.npclrs  {4 "g" 5 "m" 6 "r" 7 "c" 8 "k" 9 "b" 10 "g" 11 "c" 12 "r" 13 "m" 16 "g"}
+      c.npclrs  {4 "g" 5 "m" 6 "r" 7 "c" 8 "k" 9 "b" 10 "g" 11 "c" 12 "m" 13 "m" 16 "g"}
       c.npmarks {4 "o" 5 "x" 6 "s" 7 "x" 8 "p" 9 "+" 10 "." 11 "^" 12 "." 13 "*" 16 "."})
   c)
 
@@ -63,8 +63,9 @@
       (inc! i)))
   d)
 
-(defn draw-slmmir-image [f &optional vmin vmax ncolor colorsym switch-halves]
-  (svifn vmin -0.05 vmax 1.15 ncolor 24 colorsym False switch-halves True)
+(defn draw-slmmir-image [f &optional vmin vmax ncolor colorsym switch-halves grid]
+  (svifn vmin -0.05 vmax 1.15 ncolor 24 colorsym False switch-halves True
+         grid True)
   (sv (, m n) f.shape
       lon-idx (if switch-halves
                   (+ (list (range (// n 2) n)) (list (range 0 (// n 2))))
@@ -86,7 +87,7 @@
                  :vmin vmin :vmax vmax))
   (pl.xlim (, (first x) (last x))) (pl.xticks x xticks :fontsize fs)
   (pl.ylim (, (first y) (last y))) (pl.yticks y yticks :fontsize fs)
-  (my-grid :ls ":"))
+  (when grid (my-grid :ls ":")))
 
 ;;; parse slmmir text output
 
diff --git a/methods/islet/figures/islet.hy b/methods/islet/figures/islet.hy
index 0fca082..d5fdcb2 100644
--- a/methods/islet/figures/islet.hy
+++ b/methods/islet/figures/islet.hy
@@ -1,20 +1,24 @@
 (require [amb3 [*]])
 (import amb3 [amb3 [*]]
+        poly util
         [scipy.linalg :as linalg]
         scipy.integrate
         re math sys ctypes)
 
+(if-main (sv amb3.*when-inp-verbosity* 1))
+
 (defn nelem [xb yb]
   (* (dec (len xb)) (dec (len yb))))
 
 (defn ndof [method ne np]
-  (cond [(< method 2) (* ne (** (dec np) 2))]
+  (cond [(or (< method 2) (= method 5)) (* ne (** (dec np) 2))]
         [(= method 2) ne]
         [:else (raisefmt "nope")]))
 
 (defclass Islet []
   (defn --init-- [me]
     (try (sv lib (npy.ctypeslib.load-library "libislet" ".")
+             lib.power-bound.restype ctypes.c-double
              me.lib lib)
          (except [e [Exception]]
            (print e)
@@ -27,6 +31,139 @@
     (me.lib.eval-interpolant (c-int method) (c-int np) (c-int (len x))
                              (as-ctypes x) (as-ctypes y))
     y)
+  (defn make-space-time-op-general [me xb method np tx]
+    ;; xb: eulerian mesh
+    ;; tx: general target points
+    ;; method: 0: natural, 1: offset nodal subset
+    (sv ne (dec (len xb))
+        n (* ne (dec np))
+        A (npy.zeros (, (len tx) n)))
+    (me.lib.make-space-time-op-for-nonuni-mesh
+     (as-ctypes xb) (ctypes.c-int ne)
+     (ctypes.c-int method) (ctypes.c-int np)
+     (as-ctypes tx) (ctypes.c-int (len tx))
+     (as-ctypes A))
+    A)
+  (defn make-space-time-op [me xb method np txgll]
+    (sv ne (dec (len xb))
+        n (* ne (dec np)))
+    (assert (= (len txgll) n))
+    (me.make-space-time-op-general xb method np txgll))
+  (defn make-space-time-op-rons-general [me xb np subnp offst tx]
+    ;; xb: eulerian mesh
+    ;; tx: general target points
+    ;; subnp, offst: description of offset nodal basis.
+    (assert (= (len subnp) (len offst)))
+    (assert (<= (len subnp) (// (inc np) 2)))
+    (sv ne (dec (len xb))
+        n (* ne (dec np))
+        A (npy.zeros (, (len tx) n)))
+    (me.lib.make-space-time-op-for-nonuni-mesh-rons
+     (as-ctypes xb) (ctypes.c-int ne)
+     (ctypes.c-int np) (ctypes.c-int (len subnp))
+     (as-ctypes subnp) (as-ctypes offst)
+     (as-ctypes tx) (ctypes.c-int (len tx))
+     (as-ctypes A))
+    A)
+  (defn make-space-time-op-rons [me xb np subnp offst txgll]
+    (sv ne (dec (len xb))
+        n (* ne (dec np)))
+    (assert (= (len txgll) n))
+    (me.make-space-time-op-rons-general xb np subnp offst txgll))
+  (defn make-space-time-op-rls-general [me xb np subnp tx]
+    ;; xb: eulerian mesh
+    ;; tx: general target points
+    ;; subnp: description of LS basis.
+    (assert (<= (len subnp) (// (inc np) 2)))
+    (sv ne (dec (len xb))
+        n (* ne (dec np))
+        A (npy.zeros (, (len tx) n)))
+    (me.lib.make-space-time-op-for-nonuni-mesh-rls
+     (as-ctypes xb) (ctypes.c-int ne)
+     (ctypes.c-int np) (ctypes.c-int (len subnp))
+     (as-ctypes subnp)
+     (as-ctypes tx) (ctypes.c-int (len tx))
+     (as-ctypes A))
+    A)
+  (defn make-space-time-op-rls [me xb np subnp txgll]
+    (sv ne (dec (len xb))
+        n (* ne (dec np)))
+    (assert (= (len txgll) n))
+    (me.make-space-time-op-rls-general xb np subnp txgll))
+  (defn make-space-time-op-nondiv2d [me xb yb method np phi-param dt]
+    (assert (in method (, 0 1 2 5)))
+    (sv ne (nelem xb yb)
+        n (ndof method ne np)
+        txgll (npy.zeros n)
+        tygll (npy.zeros n)
+        A (npy.zeros (, n n))
+        phi-param (array-if-not phi-param))
+    (me.lib.make-space-time-op-nondiv2d
+     (as-ctypes xb) (ctypes.c-int (dec (len xb)))
+     (as-ctypes yb) (ctypes.c-int (dec (len yb)))
+     (ctypes.c-int method) (ctypes.c-int np)
+     (as-ctypes phi-param) (ctypes.c-int (len phi-param))
+     (ctypes.c-double dt)
+     (as-ctypes txgll) (as-ctypes tygll) (as-ctypes A))
+    (, txgll tygll A))
+  (defn apply-space-time-op-nondiv2d [me xb yb method np phi-param dt src
+                                      &optional [nstep 1] [use-caas False]]
+    (assert (in method (, 0 1 2 5)))
+    (sv ne (nelem xb yb)
+        n (ndof method ne np)
+        tgt (npy.zeros n)
+        phi-param (array-if-not phi-param))
+    (assert (= n (len src)))
+    (me.lib.apply-space-time-op-nondiv2d
+     (as-ctypes xb) (ctypes.c-int (dec (len xb)))
+     (as-ctypes yb) (ctypes.c-int (dec (len yb)))
+     (ctypes.c-int method) (ctypes.c-int np)
+     (as-ctypes phi-param) (ctypes.c-int (len phi-param))
+     (ctypes.c-double dt) (ctypes.c-int nstep) (ctypes.c-bool use-caas)
+     (as-ctypes src) (as-ctypes tgt))
+    tgt)
+  (defn apply-space-time-density-op-div1d [me method manifold trajint np
+                                           xb t0 dt nstep src]
+    (assert (in method (, 0 1 2)))
+    (sv ne (dec (len xb))
+        n (len src)
+        tgt (npy.zeros n)
+        c-int ctypes.c-int c-double ctypes.c-double)
+    (me.lib.apply-space-time-density-op-div1d
+      (c-int method) (c-int manifold) (c-int trajint) (c-int np) (as-ctypes xb)
+      (c-int ne) (c-double t0) (c-double dt) (c-int nstep) (as-ctypes src)
+      (as-ctypes tgt))
+    tgt)
+  (defn div1d-se-interp [me xb np y xi]
+    (sv yi (npy.zeros (len xi)))
+    (me.lib.div1d-se-interp (as-ctypes xb) (dec (len xb)) (ctypes.c-int np)
+                            (as-ctypes y) (as-ctypes xi) (ctypes.c-int (len xi))
+                            (as-ctypes yi))
+    yi)
+  (defn div1d-cc-interp [me xb y xi]
+    (sv yi (npy.zeros (len xi)))
+    (me.lib.div1d-cc-interp (as-ctypes xb) (dec (len xb)) (as-ctypes y)
+                            (as-ctypes xi) (ctypes.c-int (len xi)) (as-ctypes yi))
+    yi)
+  (defn div1d-unittest [me]
+    (me.lib.div1d-unittest))
+  (defn calc-power-bound [me B mu-arg]
+    (when (none? me.lib) (return -1))
+    (sv np (second B.shape))
+    (me.lib.power-bound
+     (ctypes.c-int 0)
+     (as-ctypes B)
+     (ctypes.c-int np)
+     (ctypes.c-double mu-arg)))
+  (defn apply-mixer [me xb np subnp X]
+    (sv ne (dec (len xb))
+        n (* ne (dec np)))
+    (assert (= (first X.shape) n))
+    (assert (<= subnp np))
+    (me.lib.apply-mixer
+     (as-ctypes xb) (ctypes.c-int ne) (ctypes.c-int np) (ctypes.c-int subnp)
+     (as-ctypes X) (ctypes.c-int (second X.shape)))
+    X)
   (defn get-xnodes [me method np]
     (sv xnodes (npy.zeros np))
     (me.lib.get-xnodes (ctypes.c-int method) (ctypes.c-int np) (as-ctypes xnodes))
@@ -35,7 +172,136 @@
     (sv metrics (npy.zeros 3 :dtype float))
     (me.lib.calc-xnodes-metrics-from-basis-string
       (str-ctypes basis) (as-ctypes metrics))
+    metrics)
+  (defn calc-lebesgue-consts-from-basis-string [me basis]
+    (sv metrics (npy.zeros 3 :dtype float))
+    (me.lib.calc-lebesgue-consts-from-basis-string
+      (str-ctypes basis) (as-ctypes metrics))
     metrics))
 
 (defn diff [x] (- (cut x 1) (cut x 0 -1)))
 
+(if-main
+(when-inp ["islet-A" {:ne int :np int :xbopt int :txopt int :dx float}]
+  (sv islet (Islet)
+      xb (case/eq xbopt
+                  [0
+                   (npy.linspace 0 1 (inc ne))]
+                  [1
+                   (sort (util.np-pad 0 (npy.random.rand (dec ne)) 1))]
+                  [2
+                   (sv orig (npy.linspace 0 1 (inc ne))
+                       perturb (npy.random.rand (inc ne))
+                       perturb (* 0.75 (/ (- perturb 0.5) ne)))
+                   (+ orig perturb)])
+      L (- (last xb) (first xb))
+      n (* ne (dec np))
+      txgll (case/eq txopt
+                     [0 ;; constant offset
+                      (sv y (util.make-gll-mesh xb np)
+                          y (+ y dx))
+                      y]
+                     [1 ;; random offset at element level
+                      (sv y (sort (npy.random.rand (inc ne)))
+                          y (* y (/ L (- (last y) (first y))))
+                          y (+ dx (util.make-gll-mesh y np)))
+                      y]
+                     [2 ;; random u, halfway to crossing characteristics
+                      (sv
+                       y (util.make-gll-mesh xb np)
+                       u (npy.random.rand n)
+                       dt (- (* 0.5 (max (filter (fn [e] (< e 0))
+                                                 (/ (diff y) (diff u))))))
+                       y (+ y dx (* dt u)))
+                      y]
+                     [3 ;; sin(x)
+                      (sv y (util.make-gll-mesh xb np :end True)
+                          u (npy.sin (* 2 math.pi L y))
+                          dt 0.09
+                          y (+ y dx (* dt u)))
+                      (cut y 0 -1)]))
+  (sv delta (diff txgll)
+      t-min-diff (min delta)
+      t-max-diff (max delta))
+  (assert (> t-min-diff 0))
+  (sv delta (diff xb)
+      min-diff (min delta)
+      max-diff (max delta))
+  (when (<= ne 3)
+    (print xb)
+    (print txgll)
+    (print (diff txgll)))
+  (for [method (, 1 0)]
+    (sv A (islet.make-space-time-op xb method np txgll)
+        (, lams V) (linalg.eig A :right True)
+        meam1 (- (max (npy.abs lams)) 1)
+        condv (util.cond-2norm V))
+    (prf "{:2d} {:d} {:d} {:d} diff {:5.1f} {:4.1f} meam1: {:8.1e} condv {:8.1e} {:s}"
+         np method xbopt txopt (/ t-max-diff t-min-diff) (/ max-diff min-diff)
+         meam1 condv (if (> meam1 1e-14) "BAD" "g"))))
+
+(when-inp ["test-libislet"]
+  (sv islet (Islet))
+  (islet.unittest)
+  (sv ne 7
+      np 12
+      xb (util.make-cells ne :pattern 'random)
+      tx (util.make-gll-mesh (util.make-cells (inc ne) :pattern 'random) np)
+      A1 (islet.make-space-time-op-general xb 1 np tx)
+      int (npy.dtype "int32")
+      A2 (islet.make-space-time-op-rons-general
+          xb np
+          (npy.array [9 9 10 10 9 10] :dtype int)
+          (npy.array [0 0 0 0 1 1] :dtype int) tx))
+  (prf "reldif {:1.3e}" (reldif A1 A2))
+  (dont with [(pl-plot (, 6 9) "fig/test-libislet")]
+    (pl.subplot 2 1 1) (pl.imshow A1 :interpolation "none") (pl.colorbar)
+    (pl.subplot 2 1 2) ;;(pl.imshow A2 :interpolation "none") (pl.colorbar)
+    (pl.imshow (npy.log10 (npy.abs (- A2 A1))) :interpolation "none") (pl.colorbar)))
+
+(when-inp ["test-pbc"]
+  (sv islet (Islet))
+  (print (islet.calc-power-bound (npy.random.randn 3 4) 0.2)))
+
+(when-inp ["foo"]
+  (sv x (npy.zeros (, 2 2))
+      c (as-ctypes x))
+  (print x)
+  (print c)
+  (sv (get x 0) 3)
+  (print x)
+  (print (npy.ctypeslib.as-array c))
+  (print (as-ctypes (npy.array [12 9 10 10 9 10] :dtype int)))
+  (print (as-ctypes (npy.array [12 9 10 10 9 10] :dtype (npy.dtype "int32")))))
+
+(when-inp ["plot-bases" {:method int :np int :nx int :fname str}]
+  ;; 0 gll_natural, 1 gll_offset_nodal_subset, 2 xnodal, 3 gll_best,
+  ;; 4 uniform_offset_nodal_subset
+  (sv x (npy.linspace -1 1 nx)
+      islet (Islet)
+      y (.transpose (islet.eval method np x)))
+  (with [(pl-plot (, 4 4) fname)]
+    (for [i (range np)]
+      (pl.plot x (get y i) "-"))
+    (my-grid)
+    (axis-tight-pad :pad 0)
+    (sv d 1.04)
+    (pl.xlim (, (- d) d))))
+
+(when-inp ["plot-np4" {:nx int :fname str}]
+  (sv np 4
+      x (npy.linspace -1 1 nx)
+      islet (Islet)
+      (, yn yo yb) (lfor m (, 0 1 3) (.transpose (islet.eval m np x)))
+      clrs "bgrk")
+  (with [(pl-plot (, 4 4) fname)]
+    (for [i (range np)]
+      (sv c (nth clrs i))
+      (pl.plot x (get yn i) (+ c ":"))
+      (pl.plot x (get yo i) (+ c "--"))
+      (pl.plot x (get yb i) (+ c "-")))
+    (my-grid)
+    (axis-tight-pad :pad 0)
+    (sv d 1.04)
+    (pl.xlim (, (- d) d))))
+)
diff --git a/methods/islet/figures/poly.hy b/methods/islet/figures/poly.hy
index e1e2a2d..e629af6 100644
--- a/methods/islet/figures/poly.hy
+++ b/methods/islet/figures/poly.hy
@@ -187,4 +187,15 @@
  (when-inp ["test-gll-w"]
    (for [np (range 1 8)]
      (sv w (get-gll-w np))
-     (assert (<= (reldif 2 (sum w)) (* 1 (epsilon)))))))
+     (assert (<= (reldif 2 (sum w)) (* 1 (epsilon))))))
+
+ (when-inp ["spacing"]
+   (sv nps (list (range 4 13)) dx-min [])
+   (for [np nps]
+     (.append dx-min (npy.min (npy.diff (get-gll-x np)))))
+   (for [i (range (len nps))]
+     (sv actual (/ (first dx-min) (nth dx-min i))
+         predicted (** (/ (nth nps i) (first nps)) 2))
+     (prf "{:2d} {:1.3f} {:5.2f} {:4.1f} {:6.2f}" (nth nps i) (nth dx-min i)
+          actual predicted (* 100 (/ predicted actual)))))
+)
diff --git a/methods/islet/figures/util.hy b/methods/islet/figures/util.hy
new file mode 100644
index 0000000..9812244
--- /dev/null
+++ b/methods/islet/figures/util.hy
@@ -0,0 +1,218 @@
+(require [amb3 [*]])
+(import amb3 [amb3 [*]]
+        poly
+        [scipy.linalg :as linalg])
+
+(defn make-gll-mesh [xb np &optional [end False]]
+  (sv ne (dec (len xb))
+      d (dec np)
+      n (* ne d)
+      xgll (npy.zeros (if end (inc n) n))
+      ref (* 0.5 (+ 1 (npy.array (cut (poly.get-gll-x np) 0 -1)))))
+  (for [ie (range ne)]
+    (sv dx (- (get xb (inc ie)) (get xb ie))
+        (get xgll (slice (* ie d) (* (inc ie) d))) (+ (get xb ie) (* dx ref))))
+  (when end (sv (get xgll n) (+ (- (last xb) (first xb)) (first xgll))))
+  xgll)
+
+(defn cond-2norm [A]
+  (sv s (linalg.svd A :compute-uv False))
+  (/ (max s) (min s)))
+
+(defn np-pad [a arr b]
+  (npy.append a (npy.append arr b)))
+
+(defn make-cells [ncell &optional [pattern 'uniform] [perturb-factor 0.5]]
+  (defn dx->xb [dx]
+    (sv xb (npy.cumsum dx)
+        xb (npy.append 0 xb)
+        xb (/ xb (get xb -1)))
+    xb)
+  (defn make-sin [perturb-factor]
+    (sv dx (npy.sin (* 1 math.pi (npy.linspace 0 1 ncell)))
+        dx (+ 2 dx)
+        dx (* (+ 1 (* perturb-factor (npy.random.rand ncell))) dx))
+    (dx->xb dx))
+  (assert (> ncell 1))
+  (case/eq pattern
+           ['perturb
+            (sv r (npy.random.rand (dec ncell))
+                xb (+ (npy.linspace 0 1 (inc ncell))
+                      (np-pad 0 (* perturb-factor (- r 0.5) (/ ncell)) 0)))]
+           ['sin (sv xb (make-sin perturb-factor))]
+           ['random
+            (sv xb (np-pad 0 (npy.random.rand (dec ncell)) 1))
+            (.sort xb)]
+           ['random-refine
+            (sv xb (make-cells (// ncell 2) 'random)
+                xb (npy.concatenate (, xb (cb->cc xb))))
+            (.sort xb)]
+           ['uniform
+            (sv xb (npy.linspace 0 1 (inc ncell)))]
+           ['uniformbut1
+            (sv xb (npy.linspace 0 1 (inc ncell))
+                (get xb 1) (+ (get xb 1) (/ 0.5 ncell)))]
+           ['alternating
+            (sv nper 3
+                dx (npy.zeros (, nper (// (+ ncell nper -1) nper)))
+                (get dx (, 0 (s-all))) 0.75
+                (get dx (, 1 (s-all))) 0.65
+                (get dx (, 2 (s-all))) 0.6
+                dx (cut (.reshape (.transpose dx) (npy.prod dx.shape)) 0 ncell)
+                xb (dx->xb dx))]
+           ['alternating-perturb
+            (sv nper 3
+                dx (npy.zeros (, nper (// (+ ncell nper -1) nper))))
+            (for [i (range nper)]
+              (sv (get dx (, i (s-all)))
+                  (+ 1 (* 2 (npy.random.rand) perturb-factor))))
+            (sv dx (cut (.reshape (.transpose dx) (npy.prod dx.shape)) 0 ncell)
+                xb (dx->xb dx))]
+           [:else (raisefmt "Not a pattern: {:s}" (name pattern))])
+  (when False
+    (sv dx-max (npy.max (npy.diff xb))
+        dx-min (npy.min (npy.diff xb)))
+    (print dx-max dx-min (/ dx-max dx-min)))
+  (assert (npy.all (> (cut xb 1) (cut xb 0 -1))))
+  xb)
+
+(defn make-ic [x shape-list &optional [width 0.2]]
+  (in-require (np-array? x))
+
+  (when (in 'locality shape-list)
+    (when (> (len shape-list) 1)
+      (raisefmt "make-ic: locality can be used only on its own."))
+    (sv width 0.2
+        hw (/ width 2)))
+
+  (sv bot 0.1
+      top 1.0
+      y (* bot (npy.ones (len x)))
+      nshapes (len shape-list)
+      hs (/ (* 2 nshapes))
+      hw (/ width 2)
+      x0 (* 1 hs))
+  (defn draw-trig [ctr]
+    (sv (get y (s-all)) (+ 1 (npy.sin (* 2 math.pi x)))))
+  (defn make-mask [ctr]
+    (* (>= x (- ctr hw)) (<= x (+ ctr hw))))
+  (defn draw-delta [ctr]
+    (sv mask (>= x ctr)
+        i (first (get (npy.array (range (len x))) mask))
+        (get y i) top))
+  (defn draw-rect [ctr]
+    (sv (get y (make-mask ctr)) (* 0.8 top)))
+  (defn draw-tri [ctr]
+    (sv left  (* (>= x (- ctr hw)) (<= x ctr))
+        right (* (<= x (+ ctr hw)) (>= x ctr))
+        lx (get x left) rx (get x right)
+        l0 (- ctr hw) r0 ctr
+        slope (/ (- top bot) hw)
+        (get y left)  (+ bot (* slope (- lx l0)))
+        (get y right) (+ top (* slope (- r0 rx)))))
+  (defn draw-sin [ctr]
+    (sv mask (make-mask ctr)
+        xm (get x mask)
+        cosxm (npy.cos (* math.pi (/ (- xm ctr) hw)))
+        (get y mask) (+ bot
+                        (* 0.5 (- (* 0.95 top) bot)
+                           (+ 1 cosxm)))))
+  (defn draw-gauss [ctr]
+    (sv (cut y 0)
+        (+ y
+           (* (- top bot)
+              (npy.exp (- (/ (** (- x ctr) 2)
+                            (** (* 0.27 width) 2))))))))
+  (defn draw-locality [ctr]
+    (sv mask (make-mask ctr)
+        not-mask (= mask False)
+        (get y mask) (* 1 top)
+        (get y not-mask) (* 0.25
+                            (+ 1 (npy.sin
+                                  (* 2 math.pi (get x not-mask)))))))
+  ;; treat the line as a circle embedded in a plane from which a 2D gaussian is
+  ;; sampled. this gives a truly C^inf periodic function with substantial
+  ;; bandwidth (the second unlike trig).
+  (defn draw-gauss2d [ctr]
+    (sv theta (* 2 math.pi x)
+        theta-ctr (* 2 math.pi ctr)
+        u (npy.cos theta) v (npy.sin theta)
+        u0 (npy.cos theta-ctr) v0 (npy.sin theta-ctr))
+    (sv (cut y 0)
+        (+ y
+           (* (- top bot)
+              (npy.exp (- (/ (+ (** (- u u0) 2) (** (- v v0) 2))
+                             (** (* 2 width) 2))))))))
+
+  (for [(, i shape) (enumerate shape-list)]
+    ((case/eq shape
+              ['delta draw-delta]
+              ['rect draw-rect]
+              ['tri draw-tri]
+              ['sin draw-sin]
+              ['gauss draw-gauss]
+              ['gauss2d draw-gauss2d]
+              ['trig draw-trig]
+              ['locality
+               (sv width 0.2 hw (* 0.5 width))
+               draw-locality]
+              [:else (raise (Exception "Not a shape."))])
+     (+ x0 (* i 2 hs))))
+  y)
+
+(defn cc [x]
+  (* 0.5 (+ (cut x 0 -1) (cut x 1))))
+
+(defn duplicate-dg-data-as-cg [data ne &optional [end False]]
+  (sv np (len data)
+      n (* (dec np) ne)
+      a (npy.zeros (if end (inc n) n)))
+  (for [ie (range ne)]
+    (sv (get a (slice (* ie (dec np)) (* (inc ie) (dec np)))) (cut data 0 -1)))
+  (when end
+    (sv (get a -1) (last data)))
+  a)
+
+(defn get-domain-of-dependence [A]
+  (sv n (get A.shape 1)
+      dod [])
+  (for [row A]
+    (sv nz [])
+    (for [j (range n)]
+      (if (!= (get row j) 0)
+        (.append nz j)))
+    (.sort nz)
+    (.append dod nz))
+  dod)
+
+(defn frame-init []
+  (pl.show :block False))
+(defn frame-flip [U &optional [clim None]]
+  (pl.clf)
+  (pl.imshow U :interpolation "none")
+  (unless (none? clim)
+    (pl.clim clim))
+  (pl.pause 0.05))
+
+(defn caas [dod w yp y]
+  (sv n (len yp)
+      l (npy.zeros n)
+      u (npy.zeros n))
+  (for [i (range n)]
+    (sv yp-i (get yp (get dod i))
+        lo (min yp-i)
+        (get l i) lo
+        hi (max yp-i)
+        (get u i) hi
+        (get y i) (min hi (max lo (get y i)))))
+  (sv mass-prev (sum (* w yp))
+      mass-curr (sum (* w y))
+      dm (- mass-prev mass-curr))
+  (if (>= dm 0)
+    (sv capacity (- u y)
+        alpha (/ dm (sum (* w capacity)))
+        y (+ y (* alpha capacity)))
+    (sv capacity (- l y)
+        alpha (/ dm (sum (* w capacity)))
+        y (+ y (* alpha capacity))))
+  y)
diff --git a/methods/islet/islet_studymetrics.cpp b/methods/islet/islet_studymetrics.cpp
index 8662384..7fef43a 100644
--- a/methods/islet/islet_studymetrics.cpp
+++ b/methods/islet/islet_studymetrics.cpp
@@ -49,13 +49,13 @@ class Basis : public UserInterpMethod {
     free_nodal = read_xnodes(nodes.get_np(), basis, xnodes);
   }
   bool is_ok () const { return ok; }
-  void eval(const Real& x, Real* const v) override {
+  void eval (const Real& x, Real* const v) override {
     ::eval(nodes, get_xnodes(), x, v);
   }
-  const Real* get_xnodes() const override {
+  const Real* get_xnodes () const override {
     return free_nodal ? xnodes : islet::get_x_gll_special(nodes.get_np());
   }
-  Int get_np() const override { return nodes.get_np(); }
+  Int get_np () const override { return nodes.get_np(); }
 };
 
 } // namespace anon
@@ -100,7 +100,7 @@ void run_thorough_diagnostics_from_basis_string (const char* basis) {
   {
     Real m[3];
     calc_xnodes_metrics_from_basis_string(basis, m);
-    printf("npm %1.4e %1.4e %1.4e\n", m[0], m[1], m[2]);
+    printf("    npm %1.4e %1.4e %1.4e\n", m[0], m[1], m[2]);
   }
   static const int ne_max = 11111;
   pum::Options po;
@@ -108,7 +108,7 @@ void run_thorough_diagnostics_from_basis_string (const char* basis) {
   po.ntrial = 33;
   po.mec_ne = 333;
   po.perturb = 0.01;
-  printf("ne,ndx_max %d po.ntrial %d po.mec_ne %d po.perturb %1.4f\n",
+  printf("    ne,ndx_max %d po.ntrial %d po.mec_ne %d po.perturb %1.4f\n",
          ne_max, po.ntrial, po.mec_ne, po.perturb);
   auto b = std::make_shared<Basis>(basis);
   if ( ! b->is_ok()) return;
diff --git a/methods/islet/search.cpp b/methods/islet/search.cpp
index 8c769e8..23c487e 100644
--- a/methods/islet/search.cpp
+++ b/methods/islet/search.cpp
@@ -560,12 +560,12 @@ struct NsbSearchAtom : public UserInterpMethod {
         }
 #endif
         for (Int j = 0; j < subnp[i]; ++j)
-          xsub[j] = x_gll[nodes[i][j]];
+          xsub[j] = x_gll[(int) nodes[i][j]];
         // Lagrange polynomial basis.
         std::fill(v, v + np, 0);
         eval_lagrange_poly(xsub, subnp[i], x, vsub);
         for (Int j = 0; j < subnp[i]; ++j)
-          v[nodes[i][j]] = vsub[j];
+          v[(int) nodes[i][j]] = vsub[j];
       }
       break;
     }
@@ -594,8 +594,8 @@ struct Count {
 
 void randperm (Int* const p, const Int n) {
   using islet::urand;
-  for (size_t i = 0; i < n; ++i) p[i] = i;
-  for (size_t i = 0; i < 5*n; ++i) {
+  for (Int i = 0; i < n; ++i) p[i] = i;
+  for (Int i = 0; i < 5*n; ++i) {
     const int j = urand()*n, k = urand()*n;
     std::swap(p[j], p[k]);
   }
@@ -659,7 +659,6 @@ static Restarter::Ptr read_restart (const int np) {
   if ( ! is.is_open()) return nullptr;
   const auto mt = std::make_shared<MetricsTracker>(np);
   const auto r = std::make_shared<Restarter>(mt, np, 0, 0);
-  int lnp;
   const bool ok = (read(is, r->np) && r->np == np &&
                    read(is, r->min_ooa) &&
                    read(is, r->dont_eval_if_below) &&
@@ -698,8 +697,8 @@ void eval (std::vector<NsbSearchAtom>& esa,
     }
     const auto& in = input_list[perm[ili]];
     auto& my_esa = esa[omp_get_thread_num()];
-    bool all_pve_wts;
-    Real wtr, xnodes[islet::np_max], metrics[3], pum_metric;
+    bool all_pve_wts = false;
+    Real wtr, metrics[3], pum_metric;
     const auto maxeigampm1 = my_esa.run(in, all_pve_wts, wtr, mt, metrics,
                                         pum_metric);
     if (maxeigampm1 > in.maxeigampm1) return;
@@ -757,7 +756,7 @@ void recur (const int np, const std::vector<ValidNodesList::Ptr>& vnls,
       for (int j = 0; j < b.subnp[i]; ++j)
         in.nodes[i][j] = b.nodes[i][j];
     ++input_list_pos;
-    if (input_list_pos == input_list.size()) {
+    if ((size_t) input_list_pos == input_list.size()) {
       if (restarter.eval_count >= restarter.dont_eval_if_below) {
         // Run a bunch of analyses in parallel.
         eval(esa, input_list, input_list_pos, mt, false, count);
@@ -832,6 +831,7 @@ static Int run (const int np, int min_ooa = -1,
   }
   if (input_list_pos > 0)
     eval(esa, input_list, input_list_pos, *mt, true, count);
+  return 0;
 }
 
 static void run () {
@@ -897,10 +897,137 @@ static void run_general_unittests () {
   std::cout << (nerr ? "FAIL" : "PASS") << " unit test\n";
 }
 
+namespace find_natural_bases {
+// Sweep over node locations and check the natural basis for positive weights
+// and, if positive, stability.
+
+struct NaturalInterp : public UserInterpMethod {
+  Int np = 0;
+  Real xnodes[islet::np_max] = {0};
+  
+  virtual void eval (const Real& x, Real* const v) override {
+    eval_lagrange_poly(xnodes, np, x, v);
+  }
+
+  virtual const Real* get_xnodes() const override { return xnodes; }
+
+  virtual Int get_np() const override { return np; }
+};
+
+struct NaturalXnodes {
+  NaturalXnodes (const Int inp, const Int insample)
+    : np(inp), tol(1e-13), ne_max(500), ndx_max(500), nsample(insample),
+      nslot(np/2-1), mt(np), verbose(false)
+  {
+    nodes.init(np, true);
+    std::vector<Int> subset(np);
+    for (int i = 0; i < np; ++i) subset[i] = i;
+    for (int i = 0; i < nodes.get_nh(); ++i) nodes.set(i, subset);
+    assert(nodes.ok_to_eval());
+
+    im.np = np;
+    im.xnodes[0] = -1;
+    im.xnodes[np-1] =  1;
+    if (np % 2 == 1) im.xnodes[np/2] = 0;
+  }
+
+  void set_verbose (const bool b) { verbose = b; }
+
+  void run () {
+    recur(0, 0);
+    print_stats();
+  }
+
+  void print_stats () {
+    printf("neval %d npassfilt %d npumeval %d n_mt_noaccept %d\n",
+           neval, npassfilt, npumeval, n_mt_noaccept);
+  }
+
+private:
+  Real tol;
+  Int np, ne_max, ndx_max, nsample, nslot;
+  Nodes nodes;
+  MetricsTracker mt;
+  NaturalInterp im;
+  MaxEigComputer max_eig_amp;
+  bool verbose;
+  Int neval = 0, npassfilt = 0, npumeval = 0, n_mt_noaccept = 0;
+
+  void recur (const Int slot, const Real x_base) {
+    if (slot == nslot) {
+      ++neval;
+
+      assert_stuff();
+
+      const auto* xnodes = im.get_xnodes();
+      Real xnodes_metric[3];
+      calc_xnodes_metrics(nodes, xnodes, xnodes_metric);
+      if ( ! mt.acceptable_metrics(nodes, xnodes, xnodes_metric)) {
+        ++n_mt_noaccept;
+        return;
+      }
+      bool all_pve_wts = false;
+      Real wtr = 0; {
+        Real wt[islet::np_max];
+        calc_weights(nodes, xnodes, wt);
+        calc_wts_metrics(np, wt, all_pve_wts, wtr);
+      }
+      if ( ! all_pve_wts) return;
+      ++npassfilt;
+
+      const Real maxeigampm1 = max_eig_amp.run(np, ne_max, ndx_max, tol, true,
+                                               &im);
+      const auto good = maxeigampm1 <= tol;
+      Real pum_metric = 1;
+      if (good) {
+        ++npumeval;
+        const Real pum_to_accept = mt.pum_to_accept(nodes, xnodes, xnodes_metric);
+        pum_metric = calc_pum_metric(im, pum_to_accept);
+        mt.update(xnodes_metric, pum_metric);
+      }
+
+      if (good || verbose) {
+        printf(" (%1.3e, %1.3e, %1.3e) %s",
+               maxeigampm1, pum_metric, wtr, good ? "good" : "");
+        printf(" x");
+        for (int i = 0; i < np; ++i)
+          printf(" %22.15e", xnodes[i]);
+        printf("\n");
+      }
+    } else {
+      const int n = std::ceil(nsample*(1 - x_base)) + 1;
+      assert(n >= 2);
+      for (int i = 1; i < n; ++i) {
+        const Real a = Real(i)/n;
+        const Real x = (1-a)*x_base + a;
+        im.xnodes[(np+1)/2+slot] = x;
+        im.xnodes[np/2-1-slot] = -x;
+        recur(slot+1, x);
+      }
+    }    
+  }
+
+  void assert_stuff () {
+#ifndef NDEBUG
+    const auto* x = im.get_xnodes();
+    for (int i = 1; i < np; ++i)
+      assert(x[i] > x[i-1]);
+    for (int i = 0; i < np; ++i)
+      assert(x[i] == -x[np-1-i]);
+#endif
+  }
+};
+
+void run (const Int np, const Int nsample) {
+  NaturalXnodes nx(np, nsample);
+  nx.run();
+}
+} // namespace find_natural_bases
+
 struct Command {
   enum Enum { unittest, findoffsetnodal, findnodal,
               finduniform, findlegendre, findcheb,
-              findnodal_given_bestosn};
+              findnodal_given_bestosn, findnatural};
   static Enum convert (const std::string& s) {
     if (s == "unittest") return unittest;
     if (s == "findoffsetnodal") return findoffsetnodal;
@@ -909,6 +1036,7 @@ struct Command {
     if (s == "finduniform") return finduniform;
     if (s == "findlegendre") return findlegendre;
     if (s == "findcheb") return findcheb;
+    if (s == "findnatural") return findnatural;
     throw std::logic_error("Not a command.");
   }
 };
@@ -955,6 +1083,12 @@ int main (int argc, char** argv) {
     find_offset_nodal_subset_bases::runall(np, NosbBasis::legendre); break;
   case Command::findcheb:
     find_offset_nodal_subset_bases::runall(np, NosbBasis::cheb); break;
+  case Command::findnatural: {
+    int nsample = 111;
+    if (argc > 3) nsample = std::atoi(argv[3]);
+    if (argc >= 1)
+      find_natural_bases::run(np, nsample);
+  } break;
   default:
     throw std::logic_error("Not a command.");
   }
diff --git a/methods/slmm/Makefile b/methods/slmm/Makefile
index cc73036..c620d66 100644
--- a/methods/slmm/Makefile
+++ b/methods/slmm/Makefile
@@ -3,7 +3,7 @@ include make.inc
 SIQK=../../siqk
 BUILD_DIAG=1
 
-CXXFLAGS += -DRELAX_TIME -DSLMM_NP_GT_4 -DSLMM_NP_MAX=16 -fopenmp -std=c++11 -I$(SIQK) -I$(PWD) -I$(KOKKOS)/include -DSIQK_TIME -Wno-unused-function -fPIC
+CXXFLAGS += -DRELAX_TIME -DSLMM_NP_GT_4 -DSLMM_NP_MAX=16 -fopenmp -std=c++14 -I$(SIQK) -I$(PWD) -I$(KOKKOS)/include -DSIQK_TIME -Wno-unused-function -fPIC
 LDFLAGS += -fopenmp
 
 OS_NAME := $(shell uname -s)
diff --git a/methods/slmm/make.depends b/methods/slmm/make.depends
index b938f11..9333968 100644
--- a/methods/slmm/make.depends
+++ b/methods/slmm/make.depends
@@ -106,21 +106,21 @@ slmmir_remap_data.o: slmmir_remap_data.cpp slmmir_remap_data.hpp \
  ../../siqk/siqk_search.hpp ../../siqk/siqk_intersect.hpp \
  ../../siqk/siqk_sqr.hpp slmmir.hpp slmm_spf.hpp slmmir_mesh.hpp \
  slmm_basis.hpp slmmir_util.hpp slmm_gll.hpp
-slmmir_remapper.o: slmmir_remapper.cpp slmm_mesh.hpp slmm_defs.hpp \
- ../../siqk/siqk.hpp ../../siqk/siqk_geometry.hpp \
- ../../siqk/siqk_defs.hpp ../../siqk/siqk_quadrature.hpp \
- ../../siqk/siqk_search.hpp ../../siqk/siqk_intersect.hpp \
- ../../siqk/siqk_sqr.hpp slmm_array.hpp slmm_array_tree.hpp \
- slmmir_remapper.hpp slmm_fit_extremum.hpp slmmir_p_refine.hpp \
- slmmir_mono_data.hpp slmm_gll.hpp slmm_basis.hpp slmm_gallery.hpp \
- slmm_time_int.hpp slmm_util.hpp slmmir_mesh.hpp slmmir_remap_data.hpp \
- slmm_nla.hpp slmmir.hpp slmm_spf.hpp slmmir_physgrid.hpp slmm_vis.hpp \
- slmmir_d2c.hpp slmmir_util.hpp
-slmmir_remapper_isl.o: slmmir_remapper_isl.cpp slmm_mesh.hpp \
+slmmir_remapper.o: slmmir_remapper.cpp slmm_mesh.hpp slmm_array_tree.hpp \
  slmm_defs.hpp ../../siqk/siqk.hpp ../../siqk/siqk_geometry.hpp \
  ../../siqk/siqk_defs.hpp ../../siqk/siqk_quadrature.hpp \
  ../../siqk/siqk_search.hpp ../../siqk/siqk_intersect.hpp \
- ../../siqk/siqk_sqr.hpp slmm_array.hpp slmm_array_tree.hpp \
+ ../../siqk/siqk_sqr.hpp slmm_array.hpp slmmir_remapper.hpp \
+ slmm_fit_extremum.hpp slmmir_p_refine.hpp slmmir_mono_data.hpp \
+ slmm_gll.hpp slmm_basis.hpp slmm_gallery.hpp slmm_time_int.hpp \
+ slmm_util.hpp slmmir_mesh.hpp slmmir_remap_data.hpp slmm_nla.hpp \
+ slmmir.hpp slmm_spf.hpp slmmir_physgrid.hpp slmm_vis.hpp slmmir_d2c.hpp \
+ slmmir_util.hpp
+slmmir_remapper_isl.o: slmmir_remapper_isl.cpp slmm_mesh.hpp \
+ slmm_array_tree.hpp slmm_defs.hpp ../../siqk/siqk.hpp \
+ ../../siqk/siqk_geometry.hpp ../../siqk/siqk_defs.hpp \
+ ../../siqk/siqk_quadrature.hpp ../../siqk/siqk_search.hpp \
+ ../../siqk/siqk_intersect.hpp ../../siqk/siqk_sqr.hpp slmm_array.hpp \
  slmm_accum.hpp slmmir_remapper.hpp slmm_fit_extremum.hpp \
  slmmir_p_refine.hpp slmmir_mono_data.hpp slmm_gll.hpp slmm_basis.hpp \
  slmm_gallery.hpp slmm_time_int.hpp slmm_util.hpp slmmir_mesh.hpp \
@@ -134,11 +134,11 @@ slmmir_snapshot.o: slmmir_snapshot.cpp slmmir_snapshot.hpp slmm_defs.hpp \
  slmm_time_int.hpp slmmir_mesh.hpp slmm_basis.hpp slmmir_p_refine.hpp \
  slmmir_mono_data.hpp slmm_gll.hpp slmmir_remap_data.hpp slmm_nla.hpp \
  slmmir.hpp slmm_spf.hpp slmm_accum.hpp
-slmmir_test.o: slmmir_test.cpp slmm_defs.hpp ../../siqk/siqk.hpp \
- ../../siqk/siqk_geometry.hpp ../../siqk/siqk_defs.hpp \
- ../../siqk/siqk_quadrature.hpp ../../siqk/siqk_search.hpp \
- ../../siqk/siqk_intersect.hpp ../../siqk/siqk_sqr.hpp slmm_array.hpp \
- slmm_mesh.hpp slmm_array_tree.hpp slmm_debug.hpp
+slmmir_test.o: slmmir_test.cpp slmmir.hpp slmm_defs.hpp \
+ ../../siqk/siqk.hpp ../../siqk/siqk_geometry.hpp \
+ ../../siqk/siqk_defs.hpp ../../siqk/siqk_quadrature.hpp \
+ ../../siqk/siqk_search.hpp ../../siqk/siqk_intersect.hpp \
+ ../../siqk/siqk_sqr.hpp slmm_array.hpp slmm_util.hpp slmm_spf.hpp
 slmmir_time_int.o: slmmir_time_int.cpp slmmir_time_int.hpp \
  slmm_time_int.hpp slmm_defs.hpp ../../siqk/siqk.hpp \
  ../../siqk/siqk_geometry.hpp ../../siqk/siqk_defs.hpp \
@@ -172,12 +172,12 @@ slmm_islet_string.o: slmm_islet_string.cpp slmm_islet.hpp slmm_gll.hpp \
  ../../siqk/siqk_quadrature.hpp ../../siqk/siqk_search.hpp \
  ../../siqk/siqk_intersect.hpp ../../siqk/siqk_sqr.hpp slmm_array.hpp \
  slmm_util.hpp
-slmm_mesh.o: slmm_mesh.cpp slmm_mesh.hpp slmm_defs.hpp \
- ../../siqk/siqk.hpp ../../siqk/siqk_geometry.hpp \
+slmm_mesh.o: slmm_mesh.cpp slmm_mesh.hpp slmm_array_tree.hpp \
+ slmm_defs.hpp ../../siqk/siqk.hpp ../../siqk/siqk_geometry.hpp \
  ../../siqk/siqk_defs.hpp ../../siqk/siqk_quadrature.hpp \
  ../../siqk/siqk_search.hpp ../../siqk/siqk_intersect.hpp \
- ../../siqk/siqk_sqr.hpp slmm_array.hpp slmm_array_tree.hpp slmm_gll.hpp \
- slmm_basis.hpp slmm_util.hpp slmm_islet.hpp
+ ../../siqk/siqk_sqr.hpp slmm_array.hpp slmm_gll.hpp slmm_basis.hpp \
+ slmm_util.hpp slmm_islet.hpp
 slmm_nla.o: slmm_nla.cpp slmm_util.hpp slmm_defs.hpp ../../siqk/siqk.hpp \
  ../../siqk/siqk_geometry.hpp ../../siqk/siqk_defs.hpp \
  ../../siqk/siqk_quadrature.hpp ../../siqk/siqk_search.hpp \
diff --git a/methods/slmm/slmm_defs.hpp b/methods/slmm/slmm_defs.hpp
index f0908da..6e2eaa8 100644
--- a/methods/slmm/slmm_defs.hpp
+++ b/methods/slmm/slmm_defs.hpp
@@ -1,13 +1,20 @@
 #ifndef INCLUDE_SLMM_DEFS_HPP
 #define INCLUDE_SLMM_DEFS_HPP
 
+#include "siqk.hpp"
+#include "slmm_array.hpp"
 #ifdef _OPENMP
 # include <omp.h>
 #endif
-#include "siqk.hpp"
-#include "slmm_array.hpp"
+
+#define ompparfor  _Pragma("omp parallel for")
+#define omppar     _Pragma("omp parallel")
+#define ompfor     _Pragma("omp for")
+#define ompbarrier _Pragma("omp barrier")
+#define ompsingle  _Pragma("omp single")
 
 namespace slmm {
+
 using siqk::Int;
 using siqk::Real;
 typedef Int Size;
@@ -48,6 +55,7 @@ inline ConstVec3s::HostMirror a2ConstVec3s (const Array2Dim<Real,3>& a)
 { return Vec3s::HostMirror(const_cast<Real*>(a.data()), nslices(a)); }
 inline ConstIdxs::HostMirror a2ConstIdxs (const Array2D<Int>& a)
 { return Idxs::HostMirror(const_cast<Int*>(a.data()), nslices(a), szslice(a)); }
+
 } // namespace slmm
 
 #endif
diff --git a/methods/slmm/slmm_gallery.cpp b/methods/slmm/slmm_gallery.cpp
index db62d07..8b65d00 100644
--- a/methods/slmm/slmm_gallery.cpp
+++ b/methods/slmm/slmm_gallery.cpp
@@ -6,11 +6,12 @@ namespace gallery {
 const char* InitialCondition::inputs[] =
   {"xyztrig", "gaussianhills", "cosinebells", "slottedcylinders",
    "correlatedcosinebells", "constant", "vortextracer", "toychem1", "toychem2",
-   "slotcyltrig", "smoothbelts", "cbandsc", "zero"};
+   "slotcyltrig", "smoothbelts", "cbandsc", "zero", "equatorstep",
+   "equatorsmoothstep"};
 
 const char* WindFieldType::inputs[] =
   {"dcmip1d3ll", "nondivergent", "divergent", "rotate", "nondivergenthack",
-   "movingvortices"};
+   "movingvortices", "testfn"};
 
 std::string InitialCondition::get_inputs () {
   return slmm::format_strings_as_list(inputs, 10);
@@ -19,19 +20,21 @@ std::string InitialCondition::get_inputs () {
 InitialCondition::Shape InitialCondition::from_string (const std::string& si) {
   std::string s(si);
   slmm::tolower(s);
-  if (s == inputs[0]) return XYZTrig;
-  if (s == inputs[1]) return GaussianHills;
-  if (s == inputs[2]) return CosineBells;
-  if (s == inputs[3]) return SlottedCylinders;
-  if (s == inputs[4]) return CorrelatedCosineBells;
-  if (s == inputs[5]) return Constant;
-  if (s == inputs[6]) return VortexTracer;
-  if (s == inputs[7]) return ToyChem1;
-  if (s == inputs[8]) return ToyChem2;
-  if (s == inputs[9]) return SlotCylTrig;
+  if (s == inputs[ 0]) return XYZTrig;
+  if (s == inputs[ 1]) return GaussianHills;
+  if (s == inputs[ 2]) return CosineBells;
+  if (s == inputs[ 3]) return SlottedCylinders;
+  if (s == inputs[ 4]) return CorrelatedCosineBells;
+  if (s == inputs[ 5]) return Constant;
+  if (s == inputs[ 6]) return VortexTracer;
+  if (s == inputs[ 7]) return ToyChem1;
+  if (s == inputs[ 8]) return ToyChem2;
+  if (s == inputs[ 9]) return SlotCylTrig;
   if (s == inputs[10]) return SmoothBelts;
   if (s == inputs[11]) return CBandSC;
   if (s == inputs[12]) return Zero;
+  if (s == inputs[13]) return EquatorStep;
+  if (s == inputs[14]) return EquatorSmoothStep;
   throw std::runtime_error(si + " is not an initial condition.");
 }
 
@@ -50,6 +53,8 @@ const char* InitialCondition::to_string (const Shape& shape) {
   case SmoothBelts: return inputs[10];
   case CBandSC: return inputs[11];
   case Zero: return inputs[12];
+  case EquatorStep: return inputs[13];
+  case EquatorSmoothStep: return inputs[14];
   }
   throw std::runtime_error("Should not be here.");
 }
@@ -154,6 +159,19 @@ void InitialCondition::init (const Shape shape, const Size n, const Real* const
     for (Size i = 0; i < n; ++i)
       u[i] = 0;
   } break;
+  case EquatorStep: {
+    for (Size i = 0; i < n; ++i)
+      u[i] = lat[i] >= 0 ? 1 : 0.1;
+  } break;
+  case EquatorSmoothStep: {
+    static const Real lat_thr = M_PI/4, a = 0.1, b = 1;
+    for (Size i = 0; i < n; ++i) {
+      if (std::abs(lat[i]) >= lat_thr)
+        u[i] = lat[i] >= 0 ? b : a;
+      else
+        u[i] = a + ((b-a)/2)*(1 + std::sin(M_PI/2*(lat[i]/lat_thr)));
+    }
+  } break;
   case ToyChem1: {
     for (Size i = 0; i < n; ++i) {
       Real k1, k2;
@@ -254,5 +272,250 @@ void ToyChem::tendency (const Real& lat, const Real& lon,
   cl2_f = -cl_f; // no /2 b/c there's an implicit conversion back to mixing ratio
 }
 
+// Convert from (u,v), where u is velocity along latitude and v is velocity
+// along longitude, to (x,y,z), which is velocity in the global cartesian
+// coordinate system. Add a w (local vertical) component to push the position
+// (X,Y,Z) back to the unit sphere.
+static void uv2xyz (
+  const Real X, const Real Y, const Real Z, // position
+  const Real u, const Real v, // velocity in tangent coord system
+  Real& x, Real& y, Real& z) // velocity in global coord system
+{
+  // r should be 1 but will numerically drift, so measure it ...
+  const Real r = std::sqrt(X*X + Y*Y + Z*Z);
+  // ... and then add a local vertical velocity to project back to the sphere.
+  const Real w = (1 - r)/slmm::consts::earth_radius_m;
+  Real R[9]; // Row major.
+  // The local vertical is just the position vector.
+  R[2] = X/r; R[5] = Y/r; R[8] = Z/r;
+  // The local along-latitude vector.
+  R[0] = -Y; R[3] = X; R[6] = 0;
+  const Real den = std::sqrt(R[0]*R[0] + R[3]*R[3]);
+  R[0] /= den; R[3] /= den;
+  // Local vertical x along-latitude.
+  R[1] = R[5]*R[6] - R[8]*R[3];
+  R[4] = R[8]*R[0] - R[2]*R[6];
+  R[7] = R[2]*R[3] - R[5]*R[0];
+  // Transform.
+  x = R[0]*u + R[1]*v + R[2]*w;
+  y = R[3]*u + R[4]*v + R[5]*w;
+  z = R[6]*u + R[7]*v + R[8]*w;
+}
+
+bool Dcmip1d3llOdeFn
+::eval (const Real t, const Real* const d, Real* const f) const {
+  assert ( ! use_xyz_form());
+  const Real
+    a = M_PI/6,
+    a_ref = slmm::consts::earth_radius_m,
+    tau = 1036800,
+    u0 = 2*M_PI*a_ref/tau,
+    sina = std::sin(a),
+    cosa = std::sqrt(1 - slmm::square(sina)),
+    lat = d[0],
+    lon = d[1],
+    sinp = std::sin(lat),
+    cosp = std::cos(lat),
+    sinl = std::sin(lon),
+    cosl = std::cos(lon);
+  // In what follows,
+  //     u = u0*(cosp*cosa + sinp*cosl*sina)
+  //     v = -u0*sinl*sina
+  //     w = 0
+  //     lat_t = slmm::m2radlat(v)
+  //     lon_t = slmm::m2radlon(lat, u).
+  // For numerical reasons, write this a little differently.
+  const Real v = -u0*sinl*sina;
+  f[0] = slmm::m2radlat(v);
+  // tan(phi) is singular at the pole. We could introduce a cutoff so the wind
+  // speed is not infinite, but for now it does not matter.
+  f[1] = slmm::m2radlat(u0*(slmm::sign(cosp)*cosa +
+                            sinp*cosl*sina/std::abs(cosp)));
+  return true;
+}
+
+bool NonDivergentWindField
+::eval (const Real t, const Real* const d, Real* const f) const {
+  Real theta, lambda;
+  if (use_xyz_form())
+    xyz2ll(d[0], d[1], d[2], theta, lambda);
+  else {
+    theta = d[0];  // latitude
+    lambda = d[1]; // longitude
+  }
+  const Real
+    T = slmm::day2sec(12),
+    R = slmm::consts::earth_radius_m,
+    lambda_p = lambda - 2*M_PI*t/T,
+    costh = std::cos(theta),
+    cost = std::cos(M_PI*t/T);
+  // v
+  f[0] = 10*R/T*std::sin(2*lambda_p)*costh*cost;
+  // u
+  f[1] = R/T*(10*slmm::square(std::sin(lambda_p))*std::sin(2*theta)*cost +
+              2*M_PI*costh);
+  if (use_xyz_form())
+    uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
+  else {  
+    f[0] = slmm::m2radlat(f[0]);
+    f[1] = slmm::m2radlon(theta, f[1]);
+  }
+  return true;
+}
+
+bool DivergentWindField
+::eval (const Real t, const Real* const d, Real* const f) const {
+  Real theta, lambda;
+  if (use_xyz_form())
+    xyz2ll(d[0], d[1], d[2], theta, lambda);
+  else {
+    theta = d[0];  // latitude
+    lambda = d[1]; // longitude
+  }
+  const Real
+    T = slmm::day2sec(12),
+    R = slmm::consts::earth_radius_m,
+    lambda_p = lambda - 2*M_PI*t/T,
+    costh = std::cos(theta),
+    cost = std::cos(M_PI*t/T);
+  // v
+  f[0] = 2.5*R/T*std::sin(lambda_p)*slmm::cube(costh)*cost;
+  // u
+  f[1] = R/T*(-5*slmm::square(std::sin(0.5*lambda_p))*std::sin(2*theta)*
+              slmm::square(costh)*cost + 2*M_PI*costh);
+  if (use_xyz_form())
+    uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
+  else {  
+    f[0] = slmm::m2radlat(f[0]);
+    f[1] = slmm::m2radlon(theta, f[1]);
+  }
+  return true;
+}
+
+// This test is based on
+//     Nair, R. D., & Jablonowski, C. (2008). Moving vortices on the sphere: A
+//     test case for horizontal advection problems. Monthly Weather Review,
+//     136(2), 699-711.
+// However, I have modified it as follows:
+//     1. Insert a cost(pi t/T) factor so that the test reverses at the
+//        midpoint.
+//     2. Multiply the prescribed omg by 4.
+bool MovingVortices
+::eval (const Real t, const Real* const d, Real* const f) const {
+  Real theta, lambda;
+  if (use_xyz_form())
+    xyz2ll(d[0], d[1], d[2], theta, lambda);
+  else {
+    theta = d[0];
+    lambda = d[1];
+  }
+  const Real
+    T = slmm::day2sec(12),            // rotational period, in seconds
+    R = slmm::consts::earth_radius_m, // sphere radius, in meters
+    Omega = 2*M_PI/T,                 // angular velocity of background rotation, rad/s
+    cost = std::cos(M_PI*t/T),
+    // Solid body rotation factor. 1 is one rotation per T. 0 gives stationary
+    // vortices. 1e-3 is enough to stabilize the instability of the stationary
+    // vortex center being on a cell corner.
+    fac = 1,
+    u0 = Omega*R,                     // velocity due to background rotation, m/s
+    lambda_p = lambda - fac*Omega*t,  // longitudinal center of vortex, radians
+    costh = std::cos(theta),          // cosine of latitude
+    rr = 3*std::sqrt(1 - slmm::square(costh)*slmm::square(std::sin(lambda_p))),
+    rr_denom = rr / (slmm::square(rr) + slmm::square(1.0e-10)),
+    omg = 4*(1.5*std::sqrt(3.0)*u0*std::tanh(rr)*rr_denom /
+             (slmm::square(std::cosh(rr))));
+  // v
+  f[0] = omg*std::cos(lambda_p)*cost;
+  // u
+  f[1] = omg*std::sin(lambda_p)*std::sin(theta)*cost + u0*costh*fac;
+  if (use_xyz_form())
+    uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
+  else {
+    f[0] = slmm::m2radlat(f[0]);
+    f[1] = slmm::m2radlon(theta, f[1]);
+  }
+  return true;
+}
+
+bool NonDivergentWindFieldHack
+::eval (const Real t, const Real* const d, Real* const f) const {
+  Real theta, lambda;
+  if (use_xyz_form())
+    xyz2ll(d[0], d[1], d[2], theta, lambda);
+  else {
+    theta = d[0];  // latitude
+    lambda = d[1]; // longitude
+  }
+  const Real
+    T = slmm::day2sec(12),
+    R = slmm::consts::earth_radius_m,
+    lambda_p = lambda,
+    costh = std::cos(theta),
+    cost = std::cos(M_PI*t/T);
+  // v
+  f[0] = 10*R/T*std::sin(2*lambda_p)*costh*cost;
+  // u
+  f[1] = 10*R/T*slmm::square(std::sin(lambda_p))*std::sin(2*theta)*cost;
+  if (use_xyz_form())
+    uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
+  else {  
+    f[0] = slmm::m2radlat(f[0]);
+    f[1] = slmm::m2radlon(theta, f[1]);
+  }
+  return true;
+}
+
+bool TestWindField
+::eval (const Real t, const Real* const d, Real* const f) const {
+  Real theta, lambda;
+  if (use_xyz_form())
+    xyz2ll(d[0], d[1], d[2], theta, lambda);
+  else {
+    theta = d[0];  // latitude
+    lambda = d[1]; // longitude
+  }
+  const Real
+    ztop = 12000,
+    z = 0.05*ztop,
+    T = slmm::day2sec(12),
+    R = slmm::consts::earth_radius_m,
+    T0 = 300,
+    Rd = 287.04,
+    g  = 9.80616,
+    p0 = 100000,
+    H  = Rd*T0/g,
+    u0 = 2*M_PI*R/T,
+    k0 = 10*R/T,
+    omega0 = (2*23000*M_PI)/T,
+    lambda_p = lambda - 2*M_PI*t/T,
+    costh   = std::cos(theta),
+    cost    = std::cos(M_PI*t/T),
+    sinlamp = std::sin(lambda_p),
+    p    = p0*std::exp(-z/H),
+    ptop = p0*std::exp(-ztop/H),
+    bs = 0.2;
+
+  const Real
+    s = (1.0 + std::exp((ptop-p0)/(bs*ptop)) - std::exp((p-p0)/(bs*ptop)) -
+         std::exp((ptop-p)/(bs*ptop))),
+    s_p = (-std::exp((p-p0)/(bs*ptop)) + std::exp((ptop-p)/(bs*ptop)))/(bs*ptop),
+    ud = (omega0*R)*std::cos(lambda_p)*square(costh)*cost*s_p;
+
+  // v
+  f[0] = 10*R/T*std::sin(2*lambda_p)*costh*cost;
+  // u
+  f[1] = R/T*(10*slmm::square(std::sin(lambda_p))*std::sin(2*theta)*cost +
+              2*M_PI*costh) + ud;
+
+  if (use_xyz_form())
+    uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
+  else {  
+    f[0] = slmm::m2radlat(f[0]);
+    f[1] = slmm::m2radlon(theta, f[1]);
+  }
+  return true;
+}
+
 } // namespace gallery
 } // namespace slmm
diff --git a/methods/slmm/slmm_gallery.hpp b/methods/slmm/slmm_gallery.hpp
index fe6f25c..2991d16 100644
--- a/methods/slmm/slmm_gallery.hpp
+++ b/methods/slmm/slmm_gallery.hpp
@@ -29,7 +29,7 @@ class InitialCondition {
     XYZTrig, GaussianHills, CosineBells, SlottedCylinders,
     CorrelatedCosineBells, Constant, VortexTracer, ToyChem1, ToyChem2,
     SlotCylTrig, // Slotted cylinders with a trig background
-    SmoothBelts, CBandSC, Zero
+    SmoothBelts, CBandSC, Zero, EquatorStep, EquatorSmoothStep
   };
 
   static Shape from_string(const std::string& si);
@@ -42,201 +42,40 @@ class InitialCondition {
   static std::string get_inputs();
 };
 
-// Convert from (u,v), where u is velocity along latitude and v is velocity
-// along longitude, to (x,y,z), which is velocity in the global cartesian
-// coordinate system. Add a w (local vertical) component to push the position
-// (X,Y,Z) back to the unit sphere.
-inline void uv2xyz (
-  const Real X, const Real Y, const Real Z, // position
-  const Real u, const Real v, // velocity in tangent coord system
-  Real& x, Real& y, Real& z) // velocity in global coord system
-{
-  // r should be 1 but will numerically drift, so measure it ...
-  const Real r = std::sqrt(X*X + Y*Y + Z*Z);
-  // ... and then add a local vertical velocity to project back to the sphere.
-  const Real w = (1 - r)/slmm::consts::earth_radius_m;
-  Real R[9]; // Row major.
-  // The local vertical is just the position vector.
-  R[2] = X/r; R[5] = Y/r; R[8] = Z/r;
-  // The local along-latitude vector.
-  R[0] = -Y; R[3] = X; R[6] = 0;
-  const Real den = std::sqrt(R[0]*R[0] + R[3]*R[3]);
-  R[0] /= den; R[3] /= den;
-  // Local vertical x along-latitude.
-  R[1] = R[5]*R[6] - R[8]*R[3];
-  R[4] = R[8]*R[0] - R[2]*R[6];
-  R[7] = R[2]*R[3] - R[5]*R[0];
-  // Transform.
-  x = R[0]*u + R[1]*v + R[2]*w;
-  y = R[3]*u + R[4]*v + R[5]*w;
-  z = R[6]*u + R[7]*v + R[8]*w;
-}
-
 // Integrate the ODE in lat-lon space. Not good numerically in the lon direction
 // because of the poles.
 struct Dcmip1d3llOdeFn : public OdeFn {
-  bool eval (const Real t, const Real* const d, Real* const f) const override {
-    assert ( ! use_xyz_form());
-    const Real
-      a = M_PI/6,
-      a_ref = slmm::consts::earth_radius_m,
-      tau = 1036800,
-      u0 = 2*M_PI*a_ref/tau,
-      sina = std::sin(a),
-      cosa = std::sqrt(1 - slmm::square(sina)),
-      lat = d[0],
-      lon = d[1],
-      sinp = std::sin(lat),
-      cosp = std::cos(lat),
-      sinl = std::sin(lon),
-      cosl = std::cos(lon);
-    // In what follows,
-    //     u = u0*(cosp*cosa + sinp*cosl*sina)
-    //     v = -u0*sinl*sina
-    //     w = 0
-    //     lat_t = slmm::m2radlat(v)
-    //     lon_t = slmm::m2radlon(lat, u).
-    // For numerical reasons, write this a little differently.
-    const Real v = -u0*sinl*sina;
-    f[0] = slmm::m2radlat(v);
-    // tan(phi) is singular at the pole. We could introduce a cutoff so the wind
-    // speed is not infinite, but for now it does not matter.
-    f[1] = slmm::m2radlat(u0*(slmm::sign(cosp)*cosa +
-                              sinp*cosl*sina/std::abs(cosp)));
-    return true;
-  }
+  bool eval(const Real t, const Real* const d, Real* const f) const override;
 };
 
 // Also from Lauritzen et al.
 struct NonDivergentWindField : public OdeFn {
-  bool eval (const Real t, const Real* const d, Real* const f) const override {
-    Real theta, lambda;
-    if (use_xyz_form())
-      xyz2ll(d[0], d[1], d[2], theta, lambda);
-    else {
-      theta = d[0];  // latitude
-      lambda = d[1]; // longitude
-    }
-    const Real
-      T = slmm::day2sec(12),
-      R = slmm::consts::earth_radius_m,
-      lambda_p = lambda - 2*M_PI*t/T,
-      costh = std::cos(theta),
-      cost = std::cos(M_PI*t/T);
-    // v
-    f[0] = 10*R/T*std::sin(2*lambda_p)*costh*cost;
-    // u
-    f[1] = R/T*(10*slmm::square(std::sin(lambda_p))*std::sin(2*theta)*cost +
-                2*M_PI*costh);
-    if (use_xyz_form())
-      uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
-    else {  
-      f[0] = slmm::m2radlat(f[0]);
-      f[1] = slmm::m2radlon(theta, f[1]);
-    }
-    return true;
-  }
+  bool eval(const Real t, const Real* const d, Real* const f) const override;
 };
 
 // Also from Lauritzen et al.
 struct DivergentWindField : public OdeFn {
-  bool eval (const Real t, const Real* const d, Real* const f) const override {
-    Real theta, lambda;
-    if (use_xyz_form())
-      xyz2ll(d[0], d[1], d[2], theta, lambda);
-    else {
-      theta = d[0];  // latitude
-      lambda = d[1]; // longitude
-    }
-    const Real
-      T = slmm::day2sec(12),
-      R = slmm::consts::earth_radius_m,
-      lambda_p = lambda - 2*M_PI*t/T,
-      costh = std::cos(theta),
-      cost = std::cos(M_PI*t/T);
-    // v
-    f[0] = 2.5*R/T*std::sin(lambda_p)*slmm::cube(costh)*cost;
-    // u
-    f[1] = R/T*(-5*slmm::square(std::sin(0.5*lambda_p))*std::sin(2*theta)*
-                slmm::square(costh)*cost + 2*M_PI*costh);
-    if (use_xyz_form())
-      uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
-    else {  
-      f[0] = slmm::m2radlat(f[0]);
-      f[1] = slmm::m2radlon(theta, f[1]);
-    }
-    return true;
-  }
+  bool eval(const Real t, const Real* const d, Real* const f) const override;
 };
 
 struct MovingVortices : public OdeFn {
-  bool eval( const Real t, const Real* const d, Real* const f) const override {
-    Real theta, lambda;
-    if (use_xyz_form()) {
-      xyz2ll(d[0], d[1], d[2], theta, lambda);
-    }
-    else {
-      theta = d[0];
-      lambda = d[1];
-    }
-    const Real T = slmm::day2sec(12);           // rotational period, in seconds
-    const Real R = slmm::consts::earth_radius_m;// sphere radius, in meters
-    const Real Omega = 2*M_PI/T;                // angular velocity of background rotation, radians per second
-    const Real u0 = Omega*R;                    // velocity due to background rotation, meters per second
-    const Real lambda_p = lambda - Omega*t;     // longitudinal center of vortex, radians
-    const Real costh = std::cos(theta);         // cosine of latitude
-        
-    const Real rr = 3*std::sqrt(1 - slmm::square(costh)*slmm::square(std::sin(lambda_p)));
-    const Real rr_denom = rr / (slmm::square(rr) + slmm::square(1.0e-10));
-    const Real omg = 1.5*std::sqrt(3.0)*u0*std::tanh(rr) * rr_denom / (slmm::square(std::cosh(rr)));
-    // v
-    f[0] = omg * std::cos(lambda_p);
-    // u
-    f[1] = omg * std::sin(lambda_p) * std::sin(theta) + u0 * costh;
-    if (use_xyz_form())
-      uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
-    else {
-      f[0] = slmm::m2radlat(f[0]);
-      f[1] = slmm::m2radlon(theta, f[1]);
-    }
-    return true;
-  }
+  bool eval(const Real t, const Real* const d, Real* const f) const override;
 };
 
 struct NonDivergentWindFieldHack : public OdeFn {
-  bool eval (const Real t, const Real* const d, Real* const f) const override {
-    Real theta, lambda;
-    if (use_xyz_form())
-      xyz2ll(d[0], d[1], d[2], theta, lambda);
-    else {
-      theta = d[0];  // latitude
-      lambda = d[1]; // longitude
-    }
-    const Real
-      T = slmm::day2sec(12),
-      R = slmm::consts::earth_radius_m,
-      lambda_p = lambda,
-      costh = std::cos(theta),
-      cost = std::cos(M_PI*t/T);
-    // v
-    f[0] = 10*R/T*std::sin(2*lambda_p)*costh*cost;
-    // u
-    f[1] = 10*R/T*slmm::square(std::sin(lambda_p))*std::sin(2*theta)*cost;
-    if (use_xyz_form())
-      uv2xyz(d[0], d[1], d[2], f[1]/R, f[0]/R, f[0], f[1], f[2]);
-    else {  
-      f[0] = slmm::m2radlat(f[0]);
-      f[1] = slmm::m2radlon(theta, f[1]);
-    }
-    return true;
-  }
+  bool eval(const Real t, const Real* const d, Real* const f) const override;
+};
+
+// DCMIP test 1-1 with w=0 and a particular value of z.
+struct TestWindField : public OdeFn {
+  bool eval(const Real t, const Real* const d, Real* const f) const override;
 };
 
 struct WindFieldType {
   static const char* inputs[];
 public:
   enum Enum { Dcmip1d3ll, NonDivergentWindField, DivergentWindField, Rotate,
-              NonDivergentWindFieldHack, MovingVortices };
+              NonDivergentWindFieldHack, MovingVortices, TestWindField };
   static Enum from_string (const std::string& si) {
     std::string s(si);
     slmm::tolower(s);
@@ -246,6 +85,7 @@ struct WindFieldType {
     if (s == inputs[3]) return Rotate;
     if (s == inputs[4]) return NonDivergentWindFieldHack;
     if (s == inputs[5]) return MovingVortices;
+    if (s == inputs[6]) return TestWindField;
     throw std::runtime_error(si + " is not an ODE function.");
   }
   static std::string get_inputs ()
diff --git a/methods/slmm/slmm_io.cpp b/methods/slmm/slmm_io.cpp
index 68eba62..04bb436 100644
--- a/methods/slmm/slmm_io.cpp
+++ b/methods/slmm/slmm_io.cpp
@@ -1,13 +1,14 @@
 #include "slmm_io.hpp"
 #include "slmm_util.hpp"
 
+#include <fstream>
+
 #ifdef SLMM_HAVE_NETCDF
 # include <netcdf>
 #endif
 using namespace netCDF;
 
 #include <stdexcept>
-#include <fstream>
 
 namespace slmm {
 namespace io {
diff --git a/methods/slmm/slmm_islet.cpp b/methods/slmm/slmm_islet.cpp
index 6584731..768e06c 100644
--- a/methods/slmm/slmm_islet.cpp
+++ b/methods/slmm/slmm_islet.cpp
@@ -28,6 +28,55 @@ Int GllOffsetNodal::max_degree (const Int& np) const {
   return degrees[np];
 }
 
+static Real normalize_x (const Real* gll_x, const Real& x) {
+  const Real x0 = gll_x[1];
+  return (x - x0) / (1 - x0);
+}
+
+static void outer_eval (const Real* gll_x, const Real& x, Real v[4]) {
+  const Real
+    xbar = normalize_x(gll_x, gll_x[2]),
+    ooxbar = 1 / xbar,
+    ybar = 1 / (xbar - 1),
+    xn = normalize_x(gll_x, x);
+  v[0] = 0;
+  v[1] = 1 + ybar*xn*((1 - ooxbar)*xn + ooxbar - xbar);
+  v[2] = ybar*ooxbar*xn*(xn - 1);
+  v[3] = ybar*xn*(xbar - xn);
+}
+
+#if 0
+static bool np4_subgrid_eval (const Real* const x_gll, const Real& x,
+                              Real y[4]) {
+  static constexpr Real
+    alpha = 0.5527864045000416708,
+    v = 0.427*(1 + alpha),
+    x2 = 0.4472135954999579277, // 1/sqrt(5)
+    x3 = 1 - x2,
+    det = x2*x3*(x2 - x3),
+    y2 = alpha,
+    y3 = v,
+    c1 = (x3*y2 - x2*y3)/det,
+    c2 = (-x3*x3*y2 + x2*x2*y3)/det;
+  if (x < x_gll[1] || x > x_gll[2]) {
+    Real y4[4];
+    GLL::eval_lagrange_poly(4, x_gll, x, y4);
+    if (x < x_gll[1]) {
+      outer_eval(x_gll, -x, y);
+      std::swap(y[0], y[3]);
+      std::swap(y[1], y[2]);
+    } else
+      outer_eval(x_gll, x, y);
+    const Real x0 = 1 - std::abs(x);
+    const Real a = (c1*x0 + c2)*x0;
+    for (int i = 0; i < 4; ++i)
+      y[i] = a*y[i] + (1 - a)*y4[i];
+  }
+  else
+    GLL::eval_lagrange_poly(4, x_gll, x, y);
+  return true;
+}
+#else
 static bool np4_subgrid_eval (const Real* const x_gll, const Real& x,
                               Real y[4]) {
   static const Real c1 = 0.306;
@@ -45,6 +94,7 @@ static bool np4_subgrid_eval (const Real* const x_gll, const Real& x,
     GLL::eval_lagrange_poly(4, x_gll, x, y);
   return true;
 }
+#endif
 
 bool GllOffsetNodal::eval (const Int& np, const Real& x, Real* const v) const {
   const Real* xnode;
diff --git a/methods/slmm/slmm_islet_string.cpp b/methods/slmm/slmm_islet_string.cpp
index f5efea1..562b3c0 100644
--- a/methods/slmm/slmm_islet_string.cpp
+++ b/methods/slmm/slmm_islet_string.cpp
@@ -116,7 +116,7 @@ bool Nodes::init (const std::string& s) {
     ss.get();
   };
   np = read_int();
-  if (np < 2 || np >= 13) return false;
+  if (np < 2 || np > Basis::np_max || (np > 13 && np != 16)) return false;
   Int include_bdy = read_int();
   assert(include_bdy);
   init(np);
diff --git a/methods/slmm/slmm_mesh.cpp b/methods/slmm/slmm_mesh.cpp
index 78ba206..1e84077 100644
--- a/methods/slmm/slmm_mesh.cpp
+++ b/methods/slmm/slmm_mesh.cpp
@@ -105,8 +105,7 @@ static void remove_unused_vertices (AVec3s& p, AIdxs& e,
   }
 }
 
-void make_cubedsphere_mesh (AVec3s& p, AIdxs& e,
-                            const Int n) {
+void make_cubedsphere_mesh (AVec3s& p, AIdxs& e, const Int n) {
   // Transformation of the reference mesh make_planar_mesh to make each of the
   // six faces.
   const Real d = 1.0 / std::sqrt(3.0);
@@ -780,7 +779,15 @@ inline void map_sphere_coord_to_face_coord (
 }
 } // namespace
 
-Int get_cell_idx (const Int ne, Real x, Real y, Real z) {
+Int get_cell_idx (const Int ne, const Real angle, const Real* const R,
+                  Real x, Real y, Real z) {
+  if (angle != 0) {
+    // R'(x,y,z) to bring the point to the unrotated grid.
+    Real v[3] = {x,y,z};
+    x = R[0]*v[0] + R[3]*v[1] + R[6]*v[2];
+    y = R[1]*v[0] + R[4]*v[1] + R[7]*v[2];
+    z = R[2]*v[0] + R[5]*v[1] + R[8]*v[2];
+  }
   const Int face_idx = get_cube_face_idx(x, y, z);
   Real fx, fy;
   map_sphere_coord_to_face_coord(face_idx, x, y, z, fx, fy);
@@ -825,5 +832,20 @@ void make_nonuniform (AVec3s& geo_p) {
   }
 }
 
+void rotate_grid (const Real axis[3], const Real angle, Real* R, AVec3s& geo_p) {
+  form_rotation(axis, angle, R);
+  const Int n = nslices(geo_p);
+  for (Int i = 0; i < n; ++i) {
+    auto p = slice(geo_p,i);
+    Real p0[3];
+    for (int d = 0; d < 3; ++d) p0[d] = p[d];
+    for (Int d = 0; d < 3; ++d) {
+      const Real* const row = R + 3*d;
+      p[d] = siqk::SphereGeometry::dot(row, p0);
+    }
+    siqk::SphereGeometry::normalize(p);
+  }
+}
+
 } // namespace mesh
 } // namespace slmm
diff --git a/methods/slmm/slmm_mesh.hpp b/methods/slmm/slmm_mesh.hpp
index 42845f6..bd428f6 100644
--- a/methods/slmm/slmm_mesh.hpp
+++ b/methods/slmm/slmm_mesh.hpp
@@ -1,8 +1,8 @@
 #ifndef INCLUDE_SLMM_MESH_HPP
 #define INCLUDE_SLMM_MESH_HPP
 
-#include "slmm_defs.hpp"
 #include "slmm_array_tree.hpp"
+#include "slmm_mesh.hpp"
 
 namespace slmm {
 
@@ -43,6 +43,10 @@ void make_cubedsphere_mesh(
 
 void make_nonuniform(AVec3s& geo_p);
 
+// Rotate the grid given (axis, angle). Store the rotation in the row-major 3x3
+// matrix R.
+void rotate_grid(const Real axis[3], const Real angle, Real* R, AVec3s& p);
+
 // Make geometric data for a mesh that is a cubed-sphere refined into a subcell
 // mesh, by default the CGLL mesh. This is used to study p=1 transport over a
 // p>1 GLL mesh, for example. In that case, call this instead of
@@ -79,7 +83,9 @@ void get_adjacent_cells(const AIdxs& geo_c2n, AIdxArray& geo_c2cnbrs_ptr,
 
 // Given a coordinate (x,y,z) on the sphere, return the index of the cell it is
 // in for the quasiuniform cubed-sphere mesh with ne x ne faces.
-Int get_cell_idx(const Int ne, Real x, Real y, Real z);
+//   If angle != 0, then the mesh was rotated by row-major rotation matrix R.
+Int get_cell_idx(const Int ne, const Real angle, const Real* const R,
+                 Real x, Real y, Real z);
 
 // tree is a tree over geometric cells. It is stored as an Int array.
 bool make_cubedsphere_tree_over_cells(const Int ne, AIdxArray& tree,
diff --git a/methods/slmm/slmm_runtests.py b/methods/slmm/slmm_runtests.py
index e210180..f415a54 100755
--- a/methods/slmm/slmm_runtests.py
+++ b/methods/slmm/slmm_runtests.py
@@ -47,16 +47,16 @@ def print_test (cmd):
     if len(cmd) > ll:
         cmd = cmd[len(cmd)-ll+1:]
     fmt = '{{0:.<{0:d}s}}'.format(ll+1)
-    print '{:3d} '.format(print_test.ctr) + fmt.format(cmd + ' '),
+    print('{:3d} '.format(print_test.ctr) + fmt.format(cmd + ' '), end='')
     sys.stdout.flush()
 print_test.ctr = 0;
 
 def print_result (passed):
     if not passed:
-        print '***FAILED'
+        print('***FAILED')
         return 1
     else:
-        print '   PASSED'
+        print('   PASSED')
         return 0
 
 def check_passed (cmd):
@@ -77,11 +77,11 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
               and o.mo_min >= min and o.mo_max <= max)
     result = print_result(passed)
     if not passed:
-        print '   ' + cmd
-        print (('    l2 {:1.2e} cv {:1.2e} cv_gll {:1.2e} mo_min {:1.2e} mo_max {:1.2e}' +
-                ' but l2_err {:1.2e} cv {:1.2e} cv_gll {:1.2e} min {:1.2e} max {:1.2e}').
-               format(o.l2, o.cv, o.cv_gll, o.mo_min, o.mo_max,
-                      l2_err, cv, cv_gll, min, max))
+        print('   ' + cmd)
+        print(('    l2 {:1.2e} cv {:1.2e} cv_gll {:1.2e} mo_min {:1.2e} mo_max {:1.2e}' +
+               ' but l2_err {:1.2e} cv {:1.2e} cv_gll {:1.2e} min {:1.2e} max {:1.2e}').
+            format(o.l2, o.cv, o.cv_gll, o.mo_min, o.mo_max,
+                   l2_err, cv, cv_gll, min, max))
     return result
 
 p = optparse.OptionParser()
@@ -112,41 +112,41 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
 base = ('./slmmir -method {method:s} -ode divergent -ic slottedcylinders ' +
         '-ic cosinebells -ic gaussianhills -we 0 -np {np:d} -dmc f -mono {mono:s} ' +
         '-nsteps 12 -ne {ne:d}')
-nerr += check_errs(base.format(method='pcsl', np=4, ne=10, mono='qlt'),
-                   3.34e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is also done with CSL
-nerr += check_errs(base.format(method='pcsl', np=6, ne=6, mono='qlt'),
-                   3.34e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is also done with CSL
-nerr += check_errs(base.format(method='csl', np=4, ne=10, mono='qlt'),
+nerr += check_errs(base.format(method='pisl', np=4, ne=10, mono='qlt'),
+                   3.34e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is also done with ISL
+nerr += check_errs(base.format(method='pisl', np=6, ne=6, mono='qlt'),
+                   3.34e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is also done with ISL
+nerr += check_errs(base.format(method='isl', np=4, ne=10, mono='qlt'),
                    3.47e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is remapped
-nerr += check_errs(base.format(method='pcsl', np=4, ne=10, mono='qlt-pve'),
+nerr += check_errs(base.format(method='pisl', np=4, ne=10, mono='qlt-pve'),
                    3.36e-1, cv_gll=5e-14, min=0, max=2)    # >= 0 constraint only
-nerr += check_errs(base.format(method='pcsl', np=4, ne=10, mono='caas'),
-                   3.47e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is also done with CSL
-nerr += check_errs(base.format(method='csl', np=4, ne=10, mono='caas'),
+nerr += check_errs(base.format(method='pisl', np=4, ne=10, mono='caas'),
+                   3.47e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is also done with ISL
+nerr += check_errs(base.format(method='isl', np=4, ne=10, mono='caas'),
                    3.47e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is remapped
-nerr += check_errs(base.format(method='csl', np=4, ne=10, mono='mn2'),
+nerr += check_errs(base.format(method='isl', np=4, ne=10, mono='mn2'),
                    3.47e-1, cv_gll=5e-14, min=0.1, max=1)  # rho is remapped
-#   Tracer consistency test. Apply CSL to constant q but remap rho.
-nerr += check_errs('./slmmir -method csl -ode divergent -ic constant -we 0 -np 4 ' +
+#   Tracer consistency test. Apply ISL to constant q but remap rho.
+nerr += check_errs('./slmmir -method isl -ode divergent -ic constant -we 0 -np 4 ' +
                    '-dmc f -mono qlt -rit -nsteps 12 -ne 10',
                    3e-15, cv_gll=1e-13, min=0.42, max=0.42, l2_err_is_0=True)
 
 # Test ISL with p-refinement.
-base = ('./slmmir -method pcsl -ode divergent -ic gaussianhills ' +
+base = ('./slmmir -method pisl -ode divergent -ic gaussianhills ' +
         '-we 0 -np {np:d} -dmc f -mono {mono:s} ' +
         '-nsteps 12 -ne {ne:d} -timeint {timeint:s}')
 nerr += check_errs(base.format(np=12, ne=3, mono='none', timeint='interp'),
                    9.939e-3)
 nerr += check_errs(base.format(np=12, ne=3, mono='none', timeint='exact'),
                    8.793e-3)
-base = ('./slmmir -method pcsl -ode divergent -ic slottedcylinders '
+base = ('./slmmir -method pisl -ode divergent -ic slottedcylinders '
         '-we 0 -np {np:d} -dmc f -mono {mono:s} ' +
         '-nsteps 12 -ne {ne:d} -timeint interp')
 nerr += check_errs(base.format(np=12, ne=3, mono='caas', timeint='interp'),
                    2.896e-1, cv_gll=5e-14, min=0.1, max=1)
 
 # ISL with p-refinement and separate t and v meshes.
-base = ('./slmmir -method pcsl -ode divergent -ic gaussianhills ' +
+base = ('./slmmir -method pisl -ode divergent -ic gaussianhills ' +
         '-we 0 -rit -dmc {dmc:s} -mono {mono:s} -lim {lim:s} -nsteps 13 -T 12 ' +
         '-ne 6 -np 8 -timeint interp -prefine {prefine:d} -d2c')
 nerr += check_errs(base.format(prefine=0, dmc='es', mono='caas', lim='caas'),  5.968e-03, cv=2e-14)
@@ -160,6 +160,9 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
 nerr += check_errs(base.format(prefine=5, dmc='eh', mono='caas-node', lim='caas'),  5.886e-03, cv_gll=2e-14)
 # don't break the no prop preserve case
 nerr += check_errs(base.format(prefine=5, dmc='es', mono='none', lim='none'),  4.2e-03)
+# rotate the grid
+rbase = base + ' -rotate-grid'
+nerr += check_errs(rbase.format(prefine=5, dmc='eh', mono='caas-node', lim='caas'),  5.886e-03, cv_gll=2e-14)
 # GllOffsetNodal
 base += ' -basis GllOffsetNodal'
 nerr += check_errs(base.format(prefine=5, dmc='es', mono='caas', lim='caas'),  5.885e-03, cv=4e-14)
@@ -169,8 +172,8 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
 
 base = './slmmir -nsteps 12 -ne 10 -we 0 -ode divergent -ic gaussianhills '
 
-# DSS for QOF rho, CSL tracer, with QLT.
-nerr += check_errs(base + '-np 3 -d2c -method csl -dmc f -mono qlt',
+# DSS for QOF rho, ISL tracer, with QLT.
+nerr += check_errs(base + '-np 3 -d2c -method isl -dmc f -mono qlt',
                    9.05e-2, cv_gll=2e-14)
 
 # Cell-integrated method basics.
@@ -183,7 +186,7 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
                    3.18e-2, 4e-15, min=1.495e-08, max=9.518e-01)
 nerr += check_errs(base + '-np 3 -xyz -d2c', 3.64e-2, 3e-15)
 nerr += check_errs(base + '-np 4 -xyz -d2c', 1.02e-2, 8e-15)
-nerr += check_errs(base + '-np 4 -xyz -d2c -method cdg', 1.02e-2, 3e-15)
+nerr += check_errs(base + '-np 4 -xyz -d2c -method cdg', 1.02e-2, 3.5e-15)
 
 # Limiter.
 nerr += check_errs('./slmmir -nsteps 12 -ne 10 -we 0 -ode divergent ' +
@@ -207,7 +210,7 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
 # Local DMC with Homme mass definition.
 nerr += check_errs('./slmmir -nsteps 12 -ne 10 -we 0 -ode divergent ' +
                    '-ic gaussianhills -np 4 -dmc eh',
-                   9.1e-3, cv_gll=4e-15)
+                   9.1e-3, cv_gll=5e-15)
 # Global (weaker than local) DMC with Homme mass definition.
 nerr += check_errs('./slmmir -nsteps 12 -ne 10 -we 0 -ode divergent ' +
                    '-ic gaussianhills -np 4 -dmc geh',
@@ -284,9 +287,9 @@ def check_errs (cmd, l2_err, cv=10, cv_gll=10, min=-float('Inf'), max=float('Inf
 #  Test that if rho is perturbed, a constant q stays a constant.
 nerr += check_errs('./slmmir -nsteps 12 -ne 10 -np 4 -ode nondivergent ' +
                    '-ic constant -dmc ef -mono qlt -we 0 --perturb-rho 0.05',
-                   1e-14, cv_gll=5e-14, min=0.42, max=0.42, l2_err_is_0=True)
+                   1e-6, cv_gll=5e-14, min=0.42, max=0.42, l2_err_is_0=True)
 nerr += check_errs('./slmmir -nsteps 12 -ne 10 -np 4 -ode divergent ' +
                    '-ic constant -dmc ef -mono qlt -we 0 --perturb-rho 0.05',
-                   1e-14, cv_gll=5e-14, min=0.42, max=0.42, l2_err_is_0=True)
+                   1e-6, cv_gll=5e-14, min=0.42, max=0.42, l2_err_is_0=True)
 
-print '{0:d} tests failed'.format(nerr)
+print('{0:d} tests failed'.format(nerr))
diff --git a/methods/slmm/slmm_test.cpp b/methods/slmm/slmm_test.cpp
index 49268b5..6c79afa 100644
--- a/methods/slmm/slmm_test.cpp
+++ b/methods/slmm/slmm_test.cpp
@@ -17,7 +17,8 @@ struct Command {
   enum Enum {
     test_make_cubedsphere, test_make_gll_mesh, test_make_gll_subcell_mesh,
     test_gll, test_gll_2d, test_time_int, test_qp_limiter, test_face_tree,
-    test_spf, test_fit_extremum, test_nla, islet_compute, test_mass_matrix
+    test_spf, test_fit_extremum, test_nla, islet_compute,
+    test_mass_matrix
   };
   static Enum from_string (const std::string& s) {
     if (s == "test_make_cubedsphere") return test_make_cubedsphere;
@@ -86,7 +87,8 @@ static Int test_make_cubedsphere (const Input& in) {
       den = std::sqrt(den);
       for (Int j = 0; j < 3; ++j) xyz[j] /= den;
       // Check whether we recover the cell index.
-      const Int ci_calc = mesh::get_cell_idx(in.ne, xyz[0], xyz[1], xyz[2]);
+      const Int ci_calc = mesh::get_cell_idx(in.ne, 0, nullptr,
+                                             xyz[0], xyz[1], xyz[2]);
       if (ci_calc != ci) {
         ++nerr;
         success = false;
@@ -399,37 +401,22 @@ static Int test_nla (const Input& in) {
   return test_solve_kkt() + test_form_ls_op() + QrFac::unittest(11);
 }
 
-// Use sphere_tri_area method in Homme.
-template <typename CV3s>
-static Real calc_area_homme (const CV3s& v, const Int n) {
-  using sg = siqk::SphereGeometry;
-  using sg = siqk::SphereGeometry;
+static Real calc_gll_area (const Basis& b, const Int np, const Vec3s& p, const Int nv) {
+  static const Int e[] = {0,1,2,3};
+  const Real* xnode, *wt;
+  b.get_x(np, xnode);
+  b.get_w(np, wt);
   Real area = 0;
-  for (Int i = 1, ilim = n - 1; i < ilim; ++i) {
-    const Real al = sg::calc_arc_length(slice(v,0), slice(v,i));
-    const Real bl = sg::calc_arc_length(slice(v,i), slice(v,i+1));
-    const Real cl = sg::calc_arc_length(slice(v,i+1), slice(v,0));
-    const Real
-      sina = std::sin((bl+cl-al)/2),
-      sinb = std::sin((al+cl-bl)/2),
-      sinc = std::sin((al+bl-cl)/2),
-      sins = std::sin((al+bl+cl)/2),
-      a=std::sqrt((sinb*sinc)/(sins*sina)),
-      b=std::sqrt((sina*sinc)/(sins*sinb)),
-      c=std::sqrt((sina*sinb)/(sins*sinc));
-    Real
-      a1 = 2*std::atan(a),
-      b1 = 2*std::atan(b),
-      c1 = 2*std::atan(c);
-    if (a > b && a > c)
-      a1 = -2*atan(1/a);
-    else if (b > c)
-      b1 = -2*atan(1/b);
-    else 
-      c1 = -2*atan(1/c);
-    area += a1+b1+c1;
+  for (Int i = 0; i < np; ++i) {
+    for (Int j = 0; j < np; ++j) {
+      Real J[6], tmp[3];
+      siqk::sqr::impl::calc_Jacobian(p, e, xnode[i], xnode[j], J);
+      siqk::SphereGeometry::cross(J, J+3, tmp);
+      const Real jac = std::sqrt(siqk::SphereGeometry::norm2(tmp));
+      area += wt[i]*wt[j]*jac;
+    }
   }
-  return area;
+  return area/4;
 }
 
 static Int test_fit_extremum (const Input& in) {
diff --git a/methods/slmm/slmm_time_int.hpp b/methods/slmm/slmm_time_int.hpp
index 1da226a..fadc313 100644
--- a/methods/slmm/slmm_time_int.hpp
+++ b/methods/slmm/slmm_time_int.hpp
@@ -102,7 +102,7 @@ inline void aixiy (const Size n,
  *       bool eval(Real t, const Real* y, Real* f) const
  *   to evaluate f(t), the ODE at time t. Return false on failure.
  * - method
- *       record(Real t, const Real* y)
+ *       record(Real t, const Real* y) const
  *   to optionally record y(t).
  *
  * \param opts [in] Options struct.
diff --git a/methods/slmm/slmm_util.cpp b/methods/slmm/slmm_util.cpp
index 4f469a4..24961c3 100644
--- a/methods/slmm/slmm_util.cpp
+++ b/methods/slmm/slmm_util.cpp
@@ -54,4 +54,20 @@ void fill_normals (const AVec3s& p, const AIdxs& e, AVec3s& nml, AIdxs& en) {
       }
 }
 
+void form_rotation(const Real axis[3], const Real angle, Real r[9]) {
+  const Real nrm = std::sqrt(square(axis[0]) + square(axis[1]) +
+                             square(axis[2]));
+  const Real x = axis[0]/nrm, y = axis[1]/nrm, z = axis[2]/nrm, th = angle;
+  const Real cth = std::cos(th), sth = std::sin(th), omcth = 1 - cth;
+  r[0] = cth + x*x*omcth;
+  r[3] = y*x*omcth + z*sth;
+  r[6] = z*x*omcth - y*sth;
+  r[1] = x*y*omcth - z*sth;
+  r[4] = cth + y*y*omcth;
+  r[7] = z*y*omcth + x*sth;
+  r[2] = x*z*omcth + y*sth;
+  r[5] = y*z*omcth - x*sth;
+  r[8] = cth + z*z*omcth;
+}
+
 } // namespace slmm
diff --git a/methods/slmm/slmm_util.hpp b/methods/slmm/slmm_util.hpp
index 215b55a..4738225 100644
--- a/methods/slmm/slmm_util.hpp
+++ b/methods/slmm/slmm_util.hpp
@@ -4,6 +4,7 @@
 #include "slmm_defs.hpp"
 
 #include <stdexcept>
+#include <fstream>
 
 #define require(condition) do {                                         \
     if ( ! (condition)) {                                               \
@@ -77,22 +78,7 @@ inline Real reldif (const Real a, const Real b, const Real abstol = 0)
 { return std::abs(b - a)/(abstol + std::abs(a)); }
 
 // Row-major R.
-inline void form_rotation (const Real axis[3], const Real angle, Real r[9]) {
-  const Real nrm = std::sqrt(square(axis[0]) + square(axis[1]) +
-                             square(axis[2]));
-  const Real& x = axis[0] / nrm, & y = axis[1] / nrm, & z = axis[2] / nrm,
-    & th = angle;
-  const Real cth = std::cos(th), sth = std::sin(th), omcth = 1 - cth;
-  r[0] = cth + x*x*omcth;
-  r[3] = y*x*omcth + z*sth;
-  r[6] = z*x*omcth - y*sth;
-  r[1] = x*y*omcth - z*sth;
-  r[4] = cth + y*y*omcth;
-  r[7] = z*y*omcth + x*sth;
-  r[2] = x*z*omcth + y*sth;
-  r[5] = y*z*omcth - x*sth;
-  r[8] = cth + z*z*omcth;
-}
+void form_rotation(const Real axis[3], const Real angle, Real r[9]);
 
 void fill_normals(const AVec3s& p, const AIdxs& e, AVec3s& nml, AIdxs& en);
 
@@ -123,8 +109,9 @@ inline T* tin (T* const p, const char* const msg="") {
 
 inline bool
 eq (const std::string& a, const char* const b1, const char* const b2 = 0) {
-  return (a == std::string(b1) || (b2 && a == std::string(b2)) ||
-          a == std::string("-") + std::string(b1));
+  return ((a == std::string(b1)) ||
+          (b2 && a == std::string(b2)) ||
+          (a == std::string("-") + std::string(b1)));
 }
 
 std::string& tolower(std::string& s);
diff --git a/methods/slmm/slmmir.cpp b/methods/slmm/slmmir.cpp
index 393f033..9030a09 100644
--- a/methods/slmm/slmmir.cpp
+++ b/methods/slmm/slmmir.cpp
@@ -4,13 +4,13 @@
    This program runs our various SLMMIR algorithms on test cases. There are a
    number of options, as follows:
 
-   -method: Integration method, one of {ir, cdg, csl, pcsl}. Default: ir. ir is
+   -method: Integration method, one of {ir, cdg, isl, pisl}. Default: ir. ir is
      incremental remap; density Q_ki is propagated and then remapped. cdg is
      characteristic discontinuous Galerkin. csl is classical semi-Lagrangian,
      but with rho remapped; this is a technical option and in practice should
-     not be used. pcsl is pure CSL and should be used. pcsl is super fast
+     not be used. pcsl is pure ISL and should be used. pcsl is super fast
      compared with the others, so use it if you're just wanting to get a
-     solution of some sort.
+     solution of some sort. Alias pairs: (isl,csl), (pisl,pcsl).
    -ode: The ODE test case to run, one of {dcmip1d3ll, nondivergent, divergent,
      rotate, movingvortices}. Default: divergent.
    -ic: Initial condition, in {xyztrig, gaussianhills, cosinebells,
@@ -49,10 +49,12 @@
    -o: A string that is a prefix for all output files. Default: tmp/out
    -we: Output write interval, in time steps. Default: 1. Set to 0 to turn off
      output.
-   -io: I/O output type, in {netcdf, internal}. Default: netcdf. If netcdf,
+   -io-type: I/O output type, in {netcdf, internal}. Default: netcdf. If netcdf,
      write output for use in paraview; if internal, latlon output for python.
    -io-nodss: In the I/O DSS for q, inject just one value rather than summing
      according to the DSS.
+   -io-start-cycle: Wait to start output, except ICs, until this cycle. Cycle
+    count starts at 1.
    -res: Output resolution; depends on writer. If internal, half number of
      latitude points. Default: 64.
    -rit: Record scalar measurements in time. Output is a matlab file with a
@@ -62,6 +64,10 @@
    -footprint: Track ISL and CISL footprints.
    -midpoint-check: Output error analysis for midpoint of run w.r.t. to a
     high-res 1-step solution.
+   -rotate-grid: If provided, rotate the grid so that, in particular, it's not
+    aligned with the N-S poles and other features of the flows.
+   -rhot0: If provided, freeze rho at its initial value. Impl'ed for only CSL
+    methods.
 
    Debug and analysis options:
 
@@ -295,9 +301,11 @@ struct IntegrateOptions {
   // Netcdf or custom.
   IOType::Enum io_type;
   bool io_no_dss;
+  Int io_start_cycle;
   // If I/O is internal type, lat,lon grid resolution parameter.
   Int output_resolution;
   Int pg;
+  bool rhot0;
 };
 
 struct Input {
@@ -317,6 +325,7 @@ struct Input {
   bool xyz_form; // Integrate in (x,y,z) space instead of (lat,lon).
   IntegrateOptions integrate_options;
   std::string basis;
+  bool rotate_grid;
 
   // Super-duper debug/analysis options.
   struct Debug {
@@ -428,7 +437,7 @@ static void print_one_liner (const Input& in, const Output& out) {
 
 static void init_mesh (const Int np, const Int tq_order, const Int ne,
                        const Int nonunimesh, const MeshType::Enum mesh_type,
-                       Mesh& m) {
+                       Mesh::GridRotation* grid_rotation, Mesh& m) {
   m.np = MeshType::is_subcell(mesh_type) ? 2 : np;
   m.tq_order = tq_order;
   m.nonuni = nonunimesh;
@@ -440,6 +449,11 @@ static void init_mesh (const Int np, const Int tq_order, const Int ne,
     mesh::make_cubedsphere_mesh(m.geo_p, m.geo_c2n, ne);
   }
   if (m.nonuni) mesh::make_nonuniform(m.geo_p);
+  if (grid_rotation && grid_rotation->angle != 0) {
+    mesh::rotate_grid(grid_rotation->axis, grid_rotation->angle,
+                      grid_rotation->R, m.geo_p);
+    m.grid_rotation = *grid_rotation;
+  }
   mesh::make_cgll_from_geo(m.geo_p, m.geo_c2n, m.np, *m.basis, m.cgll_p, m.cgll_c2n);
   mesh::make_dgll_from_cgll(m.cgll_p, m.cgll_c2n, m.dglln2cglln, m.dgll_c2n);
   mesh::make_io_cgll_from_internal_cgll(m.cgll_p, m.cgll_c2n, m.cgll_io_c2n);
@@ -1152,7 +1166,7 @@ static void integrate (
 
   AVec3s departure_p, rho_departure_p;
   resize(departure_p, mi->nnodes());
-  if (opts.method == Method::csl)
+  if (opts.method == Method::isl)
     resize(rho_departure_p, nslices(m.geo_p));
 
   if (pg) {
@@ -1177,9 +1191,12 @@ static void integrate (
     Timer::start(Timer::ts_integrate);
     Real ts = dt*step;
     Real tf = step == last_step ? T : ts + dt;
+    const Int cycle = 1 + step / nsteps_per_12days;
 
     const bool do_io = ((ncw || iw) &&
-                        ((step+1) % write_every == 0 || step == last_step));
+                        (((step+1) % write_every == 0 &&
+                          cycle >= opts.io_start_cycle) ||
+                         step == last_step));
     const bool do_Qq = (Method::is_ci(opts.method) &&
                         (do_io || so || step == last_step));
 
@@ -1233,15 +1250,15 @@ static void integrate (
       mi->integrate(ts, tf, departure_p);
       break;
     case IntegrateOptions::bwd:
-      assert(Method::is_csl(opts.method));
+      assert(Method::is_isl(opts.method));
       mi->integrate(tf, ts, departure_p);
-      if (opts.method == Method::csl)
+      if (opts.method == Method::isl)
         mi->integrate(ts, tf, rho_departure_p);
       break;
     }
     Timer::stop(Timer::ts_integrate);
 
-    // CDG, IR, or CSL.
+    // CDG, IR, or ISL.
     Timer::start(Timer::ts_remap);
     switch (opts.method) {
     case Method::ir:
@@ -1250,31 +1267,32 @@ static void integrate (
                       tracer_p[0]->data(), tracer_p[1]->data(), nics,
                       density_p[0]->data(), density_p[1]->data());
       break;
-    case Method::csl:
-      remapper->remap(rho_departure_p, nullptr, nullptr, 0,
-                      density_p[0]->data(), density_p[1]->data());
-      remapper->csl(departure_p,
+    case Method::isl:
+      if ( ! opts.rhot0)
+        remapper->remap(rho_departure_p, nullptr, nullptr, 0,
+                        density_p[0]->data(), density_p[1]->data());
+      remapper->isl(departure_p,
                     tracer_p[0]->data(), tracer_p[1]->data(), nics,
                     density_p[0]->data(), density_p[1]->data(),
                     Filter::is_positive_only(opts.filter), false,
                     srcterm_tendency_ptr, srcterm_beg, srcterm_end, pg);
       break;
-    case Method::pcsl:
-    case Method::pcslu:
-      remapper->csl(departure_p,
+    case Method::pisl:
+    case Method::pislu:
+      remapper->isl(departure_p,
                     tracer_p[0]->data(), tracer_p[1]->data(), nics,
                     density_p[0]->data(), density_p[1]->data(),
-                    Filter::is_positive_only(opts.filter), true,
+                    Filter::is_positive_only(opts.filter), ! opts.rhot0,
                     srcterm_tendency_ptr, srcterm_beg, srcterm_end, pg);
       break;
     default: throw std::logic_error("Not a Method.");
     }
 
     if (opts.d2c)
-      dss(*d2cer, Method::is_csl(opts.method), nics, len, wrk.data(),
+      dss(*d2cer, Method::is_isl(opts.method), nics, len, wrk.data(),
           density_p[1]->data(), tracer_p[1]->data(),
-          // For pcsl* methods, csl routine takes care of rho DSS.
-          ! Method::is_pcsl(opts.method));
+          // For pisl* methods, isl routine takes care of rho DSS.
+          ! Method::is_pisl(opts.method));
 
     Timer::stop(Timer::ts_remap);
 
@@ -1309,7 +1327,7 @@ static void integrate (
     // Observer.
     if (so) {
       if (step % nsteps_per_12days == 0)
-        std::cout << "\nC cycle " << 1 + step / nsteps_per_12days << "\n";
+        std::cout << "\nC cycle " << cycle << "\n";
       const bool report = (((step + 1) % nsteps_per_12days == 0) ||
                            step == last_step);
       so->set_time(tf);
@@ -1415,7 +1433,7 @@ static Snapshot::Ptr make_midpoint_snapshot (
 {
   IntegrateOptions io;
   io.stepping = IntegrateOptions::bwd;
-  io.method = Method::pcslu; // 1 step with natural interpolant
+  io.method = Method::pislu; // 1 step with natural interpolant
   io.d2c = true;
   io.filter = Filter::none;
   io.limiter = Limiter::none;
@@ -1439,7 +1457,9 @@ static Snapshot::Ptr make_midpoint_snapshot (
   const auto& mesh = pref->m_f;
   mesh->basis = Basis::create(Basis::Type::gll);
   //const Int plus = 4, np = m.np >= 14-plus ? 16 : m.np + plus;
-  init_mesh(m.np, 4 /* unused */, ne, m.nonuni, MeshType::geometric, *mesh);
+  mesh->grid_rotation = m.grid_rotation;
+  init_mesh(m.np, 4 /* unused */, ne, m.nonuni, MeshType::geometric,
+            &mesh->grid_rotation, *mesh);
   pref->rd_f = pref->rd_t = std::make_shared<RemapData>(*mesh, io.method, io.dmc);
 
   const auto mi = MeshIntegratorFactory::create(ode, false, mesh->cgll_p);
@@ -1482,9 +1502,19 @@ static void run (const Input& in) {
   pref->experiment = in.p_refine_experiment;
 
   const auto m = std::make_shared<Mesh>();
-  m->basis = parse_basis(in.integrate_options.method, in.basis);
-  init_mesh(in.np, in.tq_order, in.ne, in.nonunimesh, in.mesh_type, *m);
-  pref->m_f = m;
+  {
+    m->basis = parse_basis(in.integrate_options.method, in.basis);
+    Mesh::GridRotation gr;
+    if (in.rotate_grid) {
+      // Some fixed but random angle.
+      gr.axis[0] = 0.11111; gr.axis[1] = -0.051515; gr.axis[2] = 1;
+      siqk::SphereGeometry::normalize(gr.axis);
+      gr.angle = 0.142314*M_PI;
+    }
+    init_mesh(in.np, in.tq_order, in.ne, in.nonunimesh, in.mesh_type,
+              in.rotate_grid ? &gr : nullptr, *m);
+    pref->m_f = m;
+  }
 
   MeshIntegrator::Ptr mi; {
     Mesh::Ptr m_coarse;
@@ -1511,6 +1541,7 @@ static void run (const Input& in) {
     case TimeInt::interpline: {
       m_coarse->np = in.timeint_v_np;
       m_coarse->nonuni = m->nonuni;
+      m_coarse->grid_rotation = m->grid_rotation;
       m_coarse->tq_order = m->tq_order;
       m_coarse->geo_p = m->geo_p;
       m_coarse->geo_nml = m->geo_nml;
@@ -1596,25 +1627,31 @@ Input::Input (int argc, char** argv)
     ne(5), nonunimesh(0), nsteps(120), write_every(1), np(4), tq_order(18),
     mesh_type(MeshType::geometric), xyz_form(false)
 {
+  assert(eq("-foo", "-foo") && eq("--foo", "-foo") && eq("--foo", "-f", "--foo") &&
+         eq("-f", "-f", "--foo"));
   program_name = argv[0];
-  integrate_options.stepping = IntegrateOptions::fwd;
-  integrate_options.d2c = false;
-  integrate_options.record_in_time = false;
-  integrate_options.check_midpoint = false;
-  integrate_options.lauritzen_diag = false;
-  integrate_options.lauritzen_diag_io = false;
-  integrate_options.perturb_rho = 0;
-  integrate_options.dmc = Dmc::none;
-  integrate_options.filter = Filter::none;
-  integrate_options.limiter = Limiter::mn2;
-  integrate_options.subcell_bounds = false;
-  integrate_options.fitext = false;
-  integrate_options.track_footprint = false;
-  integrate_options.io_type = IOType::netcdf;
-  integrate_options.io_no_dss = false;
-  integrate_options.output_resolution = 64;
-  integrate_options.pg = 0;
+  auto& iopts = integrate_options;
+  iopts.stepping = IntegrateOptions::fwd;
+  iopts.d2c = false;
+  iopts.record_in_time = false;
+  iopts.check_midpoint = false;
+  iopts.lauritzen_diag = false;
+  iopts.lauritzen_diag_io = false;
+  iopts.perturb_rho = 0;
+  iopts.dmc = Dmc::none;
+  iopts.filter = Filter::none;
+  iopts.limiter = Limiter::mn2;
+  iopts.subcell_bounds = false;
+  iopts.fitext = false;
+  iopts.track_footprint = false;
+  iopts.io_type = IOType::netcdf;
+  iopts.io_start_cycle = 0;
+  iopts.io_no_dss = false;
+  iopts.output_resolution = 64;
+  iopts.pg = 0;
+  iopts.rhot0 = false;
   p_refine_experiment = 0;
+  rotate_grid = false;
   bool tq_specified = false, method_specified = false;
   for (int i = 1; i < argc; ++i) {
     const std::string& token = argv[i];
@@ -1639,21 +1676,21 @@ Input::Input (int argc, char** argv)
       auto& ic = initial_conditions.back();
       ic.n = atoi(argv[++i]);
     } else if (eq(token, "-mono", "--monotone"))
-      integrate_options.filter = Filter::convert(argv[++i]);
+      iopts.filter = Filter::convert(argv[++i]);
     else if (eq(token, "-lim", "--limiter"))
-      integrate_options.limiter = Limiter::convert(argv[++i]);
-    else if (eq(token, "subcellbounds"))
-      integrate_options.subcell_bounds = true;
-    else if (eq(token, "fitext"))
-      integrate_options.fitext = true;
+      iopts.limiter = Limiter::convert(argv[++i]);
+    else if (eq(token, "-subcellbounds"))
+      iopts.subcell_bounds = true;
+    else if (eq(token, "-fitext"))
+      iopts.fitext = true;
     else if (eq(token, "-method", "--method")) {
       method_specified = true;
-      integrate_options.method = Method::convert(argv[++i]);
+      iopts.method = Method::convert(argv[++i]);
     }
     else if (eq(token, "-np"))
       np = atoi(argv[++i]);
     else if (eq(token, "-pg"))
-      integrate_options.pg = atoi(argv[++i]);
+      iopts.pg = atoi(argv[++i]);
     else if (eq(token, "-mesh"))
       mesh_type = MeshType::convert(argv[++i]);
     else if (eq(token, "-tq")) {
@@ -1665,35 +1702,45 @@ Input::Input (int argc, char** argv)
       nonunimesh = atoi(argv[++i]);
     else if (eq(token, "-we", "--write-every"))
       write_every = atoi(argv[++i]);
-    else if (eq(token, "-io"))
-      integrate_options.io_type = IOType::convert(argv[++i]);
+    else if (eq(token, "-io-start-cycle"))
+      iopts.io_start_cycle = atoi(argv[++i]);
+    else if (eq(token, "-io-type"))
+      iopts.io_type = IOType::convert(argv[++i]);
     else if (eq(token, "-io-nodss"))
-      integrate_options.io_no_dss = true;
+      iopts.io_no_dss = true;
     else if (eq(token, "-res"))
-      integrate_options.output_resolution = atoi(argv[++i]);
-    else if (eq(token, "-xyz", "--xyz"))
+      iopts.output_resolution = atoi(argv[++i]);
+    else if (eq(token, "-xyz"))
       xyz_form = true;
-    else if (eq(token, "-d2c", "--d2c"))
-      integrate_options.d2c = true;
-    else if (eq(token, "-dmc", "--dmc"))
-      integrate_options.dmc = Dmc::convert(argv[++i]);
+    else if (eq(token, "-d2c"))
+      iopts.d2c = true;
+    else if (eq(token, "-dmc"))
+      iopts.dmc = Dmc::convert(argv[++i]);
     else if (eq(token, "-rit", "--record_in_time"))
-      integrate_options.record_in_time = true;
-    else if (eq(token, "-midpoint-check", "--midpoint-check"))
-      integrate_options.check_midpoint = true;
-    else if (eq(token, "-lauritzen", "--lauritzen"))
-      integrate_options.lauritzen_diag = true;
-    else if (eq(token, "-lauritzen-io", "--lauritzen-io"))
-      integrate_options.lauritzen_diag_io = true;
+      iopts.record_in_time = true;
+    else if (eq(token, "-midpoint-check"))
+      iopts.check_midpoint = true;
+    else if (eq(token, "-lauritzen"))
+      iopts.lauritzen_diag = true;
+    else if (eq(token, "-lauritzen-io"))
+      iopts.lauritzen_diag_io = true;
     else if (eq(token, "-footprint", "--footprint"))
-      integrate_options.track_footprint = true;
-    else if (eq(token, "--write-binary-at"))
+      iopts.track_footprint = true;
+    else if (eq(token, "-write-binary-at"))
       debug.write_binary_at = atoi(argv[++i]);
+    else if (eq(token, "-perturb-rho"))
+      iopts.perturb_rho = atof(argv[++i]);
+    else if (eq(token, "-rotate-grid"))
+      rotate_grid = true;
+    else if (eq(token, "-rhot0"))
+      iopts.rhot0 = true;
+    else
+      SIQK_THROW_IF(true, "Unrecognized option: " << token);
   }
   if ( ! tq_specified) {
-    if (Dmc::is_facet(integrate_options.dmc))
+    if (Dmc::is_facet(iopts.dmc))
       tq_order = (np-1)*4;
-    else if ( ! Method::is_ir(integrate_options.method))
+    else if ( ! Method::is_ir(iopts.method))
       tq_order = np >= 4 ? 20 : 14;
     else
       tq_order = np == 5 ? 20 : np == 4 ? 18 : np == 3 ? 14 : 8;
@@ -1704,19 +1751,19 @@ Input::Input (int argc, char** argv)
   }
   if ( ! method_specified) {
     // Default method to the one natural for the quadrature geometry.
-    integrate_options.method = Dmc::is_facet(integrate_options.dmc) ?
+    iopts.method = Dmc::is_facet(iopts.dmc) ?
       Method::cdg : Method::ir;
   }
   if (initial_conditions.empty())
     initial_conditions.push_back(InitialCondition("xyztrig", 1));
-  if (Method::is_csl(integrate_options.method))
-    integrate_options.stepping = IntegrateOptions::bwd;
+  if (Method::is_isl(iopts.method))
+    iopts.stepping = IntegrateOptions::bwd;
   if (MeshType::is_subcell(mesh_type) &&
-      Dmc::is_facet(integrate_options.dmc) &&
-      Method::is_ir(integrate_options.method)) {
+      Dmc::is_facet(iopts.dmc) &&
+      Method::is_ir(iopts.method)) {
     std::cout << "WARNING: Switching to CDG; (QOF, IR) is not supported "
               << "for subcell mesh.\n";
-    integrate_options.method = Method::cdg;
+    iopts.method = Method::cdg;
   }
   SIQK_THROW_IF(p_refine_experiment == 1 && ! TimeInt::is_interp(timeint),
                 "For p_refine_experiment = 1, -timeint must be an interp type.");
@@ -1730,26 +1777,27 @@ Input::Input (int argc, char** argv)
       " p_refine_experiment to 0.\n";
     p_refine_experiment = 0;
   }
-  if ( ! integrate_options.record_in_time)
-    integrate_options.lauritzen_diag = false;
+  if ( ! iopts.record_in_time)
+    iopts.lauritzen_diag = false;
 #ifndef SLMMIR_LAURITZEN_DIAG
   std::cout << "Warning: Turning off Lauritzen diagnostics b/c not built.\n";
-  integrate_options.lauritzen_diag = false;
+  iopts.lauritzen_diag = false;
 #endif
   print(std::cout);
-  SIQK_THROW_IF(integrate_options.subcell_bounds &&
-                ! Method::is_csl(integrate_options.method),
-                "-subcellbounds and not CSL is not supported.");
-  SIQK_THROW_IF(integrate_options.limiter == Limiter::none &&
-                ! (integrate_options.filter == Filter::none ||
-                   integrate_options.filter == Filter::caas),
+  SIQK_THROW_IF(iopts.subcell_bounds && ! Method::is_isl(iopts.method),
+                "-subcellbounds and not ISL is not supported.");
+  SIQK_THROW_IF(iopts.limiter == Limiter::none &&
+                ! (iopts.filter == Filter::none ||
+                   iopts.filter == Filter::caas),
                 "-lim none is valid only with -mono none or caas");
-  SIQK_THROW_IF(integrate_options.pg > 0 &&
-                ( ! Dmc::use_homme_mass(integrate_options.dmc) ||
-                  ! Method::is_csl(integrate_options.method)),
-                "-pg requires Homme mass and is impl'ed for CSL only");
-  SIQK_THROW_IF(integrate_options.check_midpoint && nsteps % 2 != 0,
+  SIQK_THROW_IF(iopts.pg > 0 &&
+                ( ! Dmc::use_homme_mass(iopts.dmc) ||
+                  ! Method::is_isl(iopts.method)),
+                "-pg requires Homme mass and is impl'ed for ISL only");
+  SIQK_THROW_IF(iopts.check_midpoint && nsteps % 2 != 0,
                 "nsteps must be even to check the midpoint.");
+  SIQK_THROW_IF(iopts.rhot0 && ! Method::is_isl(iopts.method),
+                "-rhot0 is impl'ed only for ISL methods.");
 }
 
 void Input::print (std::ostream& os) const {
diff --git a/methods/slmm/slmmir.hpp b/methods/slmm/slmmir.hpp
index 6e63a8e..a9abc82 100644
--- a/methods/slmm/slmmir.hpp
+++ b/methods/slmm/slmmir.hpp
@@ -107,33 +107,37 @@ struct MeshType {
 };
 
 struct Method {
-  enum Enum { ir, cdg, csl,
-              pcsl, /* stabilized interpolation SL (ISL) */
-              pcslu /* unstable ISL for testing and comparison */ };
+  enum Enum { ir, cdg, isl,
+              pisl, /* stabilized interpolation SL (ISL) */
+              pislu /* unstable ISL for testing and comparison */ };
   static std::string convert (const Enum& e) {
     switch (e) {
     case Enum::ir: return "ir";
     case Enum::cdg: return "cdg";
-    case Enum::csl: return "csl";
-    case Enum::pcsl: return "pcsl";
-    case Enum::pcslu: return "pcslu";
+    case Enum::isl: return "isl";
+    case Enum::pisl: return "pisl";
+    case Enum::pislu: return "pislu";
     default: throw std::runtime_error("Method: Not a valid enum.");
     }
   }
   static Enum convert (const std::string& s) {
     if (eq(s, "ir")) return Enum::ir;
     if (eq(s, "cdg")) return Enum::cdg;
-    if (eq(s, "csl")) return Enum::csl;
-    if (eq(s, "pcsl")) return Enum::pcsl;
-    if (eq(s, "pcslu")) return Enum::pcslu;
+    if (eq(s, "isl")) return Enum::isl;
+    if (eq(s, "pisl")) return Enum::pisl;
+    if (eq(s, "pislu")) return Enum::pislu;
+    // Aliases.
+    if (eq(s, "csl")) return Enum::isl;
+    if (eq(s, "pcsl")) return Enum::pisl;
+    if (eq(s, "pcslu")) return Enum::pislu;
     throw std::runtime_error(std::string("Method: Not a valid string: ") + s);
   }
-  static bool is_csl (const Enum& e) { return e == csl || e == pcsl || e == pcslu; }
-  static bool is_pcsl (const Enum& e) { return e == pcsl || e == pcslu; }
-  static bool is_stabilized (const Enum& e) { return e != pcslu; }
+  static bool is_isl (const Enum& e) { return e == isl || e == pisl || e == pislu; }
+  static bool is_pisl (const Enum& e) { return e == pisl || e == pislu; }
+  static bool is_stabilized (const Enum& e) { return e != pislu; }
   static bool is_ir (const Enum& e) { return e == ir; }
   // Is cell-integrated.
-  static bool is_ci (const Enum& e) { return ! is_csl(e); }
+  static bool is_ci (const Enum& e) { return ! is_isl(e); }
 };
 
 struct Filter {
diff --git a/methods/slmm/slmmir_mesh.hpp b/methods/slmm/slmmir_mesh.hpp
index 4511f89..05e71f1 100644
--- a/methods/slmm/slmmir_mesh.hpp
+++ b/methods/slmm/slmmir_mesh.hpp
@@ -10,8 +10,14 @@ using namespace slmm;
 struct Mesh {
   typedef std::shared_ptr<Mesh> Ptr;
 
+  struct GridRotation {
+    Real axis[3], angle, R[9];
+    GridRotation () : angle(0) {}
+  };
+
   Int np, nonuni, tq_order;
   Basis::Ptr basis;
+  GridRotation grid_rotation;
   AVec3s geo_p, geo_nml, cgll_p;
   AIdxs geo_c2n, geo_c2nml, cgll_c2n, dgll_c2n, cgll_io_c2n;
   AIdxArray dglln2cglln;
diff --git a/methods/slmm/slmmir_remap_data.cpp b/methods/slmm/slmmir_remap_data.cpp
index c5c84e7..a934381 100644
--- a/methods/slmm/slmmir_remap_data.cpp
+++ b/methods/slmm/slmmir_remap_data.cpp
@@ -237,7 +237,7 @@ RemapData::RemapData (const Mesh& m, const Method::Enum method, const Dmc::Enum
   calc_gll_basis_function_integrals(m, GLL(), dgbfi_gll_);
   method_ = method;
   dmc_ = dmc;
-  if ( ! Method::is_pcsl(method)) calc_M_fwd(m, fmm_, dmc);
+  if ( ! Method::is_pisl(method)) calc_M_fwd(m, fmm_, dmc);
   init_dgbfi_mass(Dmc::use_homme_mass(dmc_) ? dgbfi_gll_ : dgbfi_);
 }
 
diff --git a/methods/slmm/slmmir_remapper.cpp b/methods/slmm/slmmir_remapper.cpp
index f52a445..6fc25ec 100644
--- a/methods/slmm/slmmir_remapper.cpp
+++ b/methods/slmm/slmmir_remapper.cpp
@@ -675,22 +675,22 @@ ::Remapper (const PRefineData::Ptr& pr, const MonoData::Ptr& md, const D2Cer::Pt
     resize(rho_src_cmbc_, dnn_);
     resize(rho_tgt_cmbc_, dnn_);
   }
-  if (Method::is_csl(rd_->method()))
-    init_csl();
+  if (Method::is_isl(rd_->method()))
+    init_isl();
 }
 
 void Remapper::
 use_subcell_bounds () {
-  SIQK_THROW_IF( ! Method::is_csl(rd_->method()),
-                 "Subcell bounds are for CSL only.");
+  SIQK_THROW_IF( ! Method::is_isl(rd_->method()),
+                 "Subcell bounds are for ISL only.");
   subcell_bounds_ = true;
   resize(q_data_, square(np_-1)*ncell_, 2);
 }
 
 void Remapper::
 use_fit_extremum () {
-  SIQK_THROW_IF( ! Method::is_csl(rd_->method()),
-                 "Quadratic extremum is impl'ed CSL only.");
+  SIQK_THROW_IF( ! Method::is_isl(rd_->method()),
+                 "Quadratic extremum is impl'ed ISL only.");
   fit_extremum_ = std::make_shared<FitExtremum>(np_);
 }
 
diff --git a/methods/slmm/slmmir_remapper.hpp b/methods/slmm/slmmir_remapper.hpp
index 324cdf9..9455567 100644
--- a/methods/slmm/slmmir_remapper.hpp
+++ b/methods/slmm/slmmir_remapper.hpp
@@ -9,6 +9,8 @@
 #include "slmmir_d2c.hpp"
 #include "slmmir.hpp"
 
+struct C2DRelations { AIdxArray ptr, r; };
+
 // Orchestrate the remap operation in one time step.
 //   "np basis" is the full nonmonotonic basis.
 //   "npmono basis" is the monotonic basis used in the cell mean boundedness
@@ -36,12 +38,15 @@ class Remapper {
     rho_src_cmbc_,      // Corrected, in the np basis.
     rho_tgt_cmbc_;      // Density in target cell given by cmbc sources.
 
-  // For CSL.
-  struct C2DRelations { AIdxArray ptr, r; };
+  // For ISL.
+public:
+  class IslImpl;
+private:
+  std::shared_ptr<IslImpl> isl_impl_;
   C2DRelations c2d_v_, c2d_f_;
   AIdxArray cn_src_cell_;
   ARealArray2 q_data_;
-  ARealArray Je_, density_factor_;
+  ARealArray Je_, Jdep_;
 
   Real perturb_rho_;
   bool record_total_mass_redistribution_;
@@ -79,11 +84,8 @@ class Remapper {
     Real* const src_tracer, Real* const tgt_tracer, const Int ntracers,
     Real* const src_density, Real* const tgt_density);
 
-  class IslImpl;
-  std::shared_ptr<IslImpl> isl_impl_;
-
-  void init_csl();
-  void init_csl_jacobian(const Mesh& m);
+  void init_isl();
+  void init_isl_jacobian(const Mesh& m);
   static void make_c2d_relations(
     const Int cnn, const AIdxArray& dglln2cglln,
     C2DRelations& c2d);
@@ -92,16 +94,16 @@ class Remapper {
     const Mesh& m, const C2DRelations& c2d, const AVec3s& advected_p,
     const Real* const src_tracer, Real* const tgt_tracer, const Int ntracers,
     const Real* const src_density, Real* const tgt_density,
-    const bool rho_csl = false);
+    const bool rho_isl = false);
 
   // Constrained density reconstruction for classical semi-Lagrangian.
-  void csl_cdr(
+  void isl_cdr(
     const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
     const MonoData& md, const spf::MassRedistributor::Ptr& mrd,
     Real* const src_density, Real* const src_tracer, const Int np_src,
     Real* const tgt_density, Real* const tgt_tracer, const Int ntracers,
-    const bool positive_only, const bool rho_csl, const Filter::Enum cdr_method);
-  void csl_cdr_rho(
+    const bool positive_only, const bool rho_isl, const Filter::Enum cdr_method);
+  void isl_cdr_rho(
     const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
     const MonoData& md, const spf::MassRedistributor::Ptr& mrd,
     Real* const src_density, const Int np2_src,
@@ -128,10 +130,10 @@ class Remapper {
              Real* const src_tracer, Real* const tgt_tracer, const Int ntracers,
              Real* const src_density, Real* const tgt_density);
     
-  void csl(const AVec3s& advected_p, Real* const src_tracer,
+  void isl(const AVec3s& advected_p, Real* const src_tracer,
            Real* const tgt_tracer, const Int ntracers,
            Real* const src_density, Real* const tgt_density,
-           const bool positive_only, const bool rho_csl,
+           const bool positive_only, const bool rho_isl,
            const Real* toychem_tendency, Int toychem_idx_beg, Int toychem_idx_end,
            const pg::PhysgridData::Ptr& pg);
 
diff --git a/methods/slmm/slmmir_remapper_isl.cpp b/methods/slmm/slmmir_remapper_isl.cpp
index cf93492..b7239fe 100644
--- a/methods/slmm/slmmir_remapper_isl.cpp
+++ b/methods/slmm/slmmir_remapper_isl.cpp
@@ -150,9 +150,9 @@ void calc_sphere_to_ref (
 } // namespace test
 
 /* The jacobian when all nodes are moved. Use the interp formula for
-   position. Let (a,b) in [-1,1]^2. Let va be the GLL basis value at a, and
+   position. Let (a,b) in [-1,1]^2. Let va be the GLL basis values at a, and
    similarly for b. Let p(:,i) be the physical location of the i'th node. Then
-       f(a,b) = sum_{i=1}^np sum_{j=1}^np p(:, np*(j-1) + i) va[i] vb[j].
+       f(a,b) = sum_{i=1}^np sum_{j=1}^np p(:, np*(j-1) + i) va[i] vb[j]
        g(a,b) = norm(f(a,b))
        q = f(a,b) / g(a,b).
    Derivatives:
@@ -160,7 +160,7 @@ void calc_sphere_to_ref (
        g_a = g_f f_a
        g_f = 1/2 (f'f)^(-1/2) 2 f = f/norm(f)
    and similarly for q_b. In this case, we have
-      f_a = sum_{i=1}^np sum_{j=1}^np p(:, np*(j-1) + i) va_a[i] vb[j],
+       f_a = sum_{i=1}^np sum_{j=1}^np p(:, np*(j-1) + i) va_a[i] vb[j],
    and va_a comes from the GLL polynomial.
 */
 template <typename Cell>
@@ -260,19 +260,18 @@ void calc_isoparametric_sphere_to_ref (
 }
 
 void Remapper::
-init_csl_jacobian (const Mesh& m) {
+init_isl_jacobian (const Mesh& m) {
   const Int dnn = nslices(m.dglln2cglln);
-  resize(density_factor_, dnn);
+  resize(Jdep_, dnn);
   resize(Je_, dnn);
   const Int np = m.np, np2 = square(np);
-  GLL gll;
-  const Real* gll_x, * gll_wt;
-  gll.get_coef(m.np, gll_x, gll_wt);
+  const Real* xnode;
+  m.basis->get_x(m.np, xnode);
   for (Int i = 0; i < dnn; ++i) {
     const Int ti = i / np2;
     const auto& cell = slice(m.cgll_c2n, ti);
     const Int ni = i % np2;
-    const Real ta = gll_x[ni % np], tb = gll_x[ni / np];
+    const Real ta = xnode[ni % np], tb = xnode[ni / np];
     if (md_) {
       // The correct, but increasingly expensive with np, method to compute
       // the density correction. Must use this if property preservation is on.
@@ -467,49 +466,6 @@ adjust_mass_nonbdy (const Real* w_s, const Real* r_s, const Int np2_s,
 #endif
 }
 
-class Source {
-  const Real* dq;
-  pg::PhysgridData::Ptr pg;
-  Int nf, nf2;
-
-  static Int idx (const Int& nf, const Real& ref_coord) {
-    const Real p = (1 + ref_coord)/2;
-    const Int i = std::max(0, std::min(nf-1, Int(nf*p)));
-    return i;
-  }
-
-public:
-  Source (const Real* dq_, pg::PhysgridData::Ptr pg_)
-    : dq(dq_), pg(pg_), nf(pg->ops->nphys), nf2(square(nf))
-  {
-    assert(idx(2,-0.11) == 0); assert(idx(2,0.11) == 1);
-    assert(idx(2,-1.01) == 0); assert(idx(2,1.01) == 1);
-    assert(idx(4,-0.24) == 1); assert(idx(4,0.74) == 3);
-  }
-
-  Real interp (const Int& cell_idx, const Real& a, const Real& b) const {
-    const Int sub_idx = nf*idx(nf, b) + idx(nf, a);
-    return dq[cell_idx*nf2 + sub_idx];
-  }
-
-  Real cell_mass (const int ic) const {
-    assert(0); return 0;
-  }
-
-  Real global_mass () const {
-    Real mass;
-    accumulate_threaded_bfb<1>(
-      [&] (const Int i, Real* a) {
-        *a += pg->ops->mesh_data->fv_metdet[i] * pg->rho[i] * dq[i];
-      },
-      len(pg->ops->mesh_data->fv_metdet), &mass);
-    mass /= square(pg->ops->nphys);
-    return mass;
-  }
-};
-
-typedef std::vector<std::shared_ptr<Source> > Sources;
-
 struct TendencyData {
   typedef std::shared_ptr<TendencyData> Ptr;
   static constexpr Real lcl_cdr_relax_fac = 1e-2;
@@ -518,7 +474,6 @@ struct TendencyData {
   const Real* dq; // mixing ratio
   Int idx_beg, idx_end; // tracers [idx_beg,idx_end) have tendencies
   Array<Real> q_min, q_max;
-  Sources sources;
 
   TendencyData () : cdr_on_src_tends(false) {}
 };
@@ -599,7 +554,6 @@ class Remapper::IslImpl {
       tend->idx_beg = srcterm_beg;
       tend->idx_end = srcterm_end;
       tend->dss_source = ! pg && tend->idx_beg < ntracers;
-      tend->sources.resize(ntracers);
       if (srcterm_tendency || tend->dss_source) {
         const auto ncell = nslices(p->m_t->geo_c2n);
         tend->q_min.optclear_and_resize(ntracers*ncell);
@@ -1120,7 +1074,7 @@ void Remapper::make_c2d_relations (
 }
 
 void Remapper::
-init_csl () {
+init_isl () {
   isl_impl_ = std::make_shared<IslImpl>(pr_);
   const auto& m = *(pr_->experiment == 0 ? pr_->m_f : pr_->m_v);
   const Int cnn = nslices(m.cgll_p);
@@ -1133,7 +1087,7 @@ init_csl () {
   resize(cn_src_cell_, std::max(cnn, cnn_f));
   // independent of mesh
   resize(q_data_, ncell_, 2);
-  init_csl_jacobian(m);
+  init_isl_jacobian(m);
 }
 
 // Always d2c rho because the density factor is discontinuous at the element
@@ -1144,107 +1098,147 @@ static void dss_rho (D2Cer& d2cer, Real* rho, Array<Real>& wrk) {
   d2cer.dss(rho, wrk.data());
 }
 
+static Int
+find_src_cell (const Mesh& m, const RemapData& rd, const Int ne,
+               const Int cni /* target node, continuous ID */,
+               const AVec3s& advected_p) {
+  // Solve for departure point.
+  const auto p = slice(advected_p, cni);
+  // Normalize p for much faster calc_sphere_to_ref.
+  Real p_sphere[3];
+  for (Int d = 0; d < 3; ++d) p_sphere[d] = p[d];
+  siqk::SphereGeometry::normalize(p_sphere);
+  Int ci;
+  if (m.nonuni) {
+    SearchFunctor sf(m, p_sphere);
+    rd.octree().apply(sf.bb, sf); // same regardless of mesh
+    ci = sf.ci_hit;
+  } else {
+    const auto& gr = m.grid_rotation;
+    ci = mesh::get_cell_idx(ne, gr.angle, gr.R, p_sphere[0], p_sphere[1],
+                            p_sphere[2]);
+  }
+  return ci;
+}
+
+static void
+calc_departure_data (const Mesh& m, const RemapData& rd, const Int ne,
+                     const Int cni /* target node, continuous ID */,
+                     const AVec3s& advected_p,
+                     Int& ci, bool calc_ci, Real* va, Real* vb, Real& a, Real& b) {
+  // Solve for departure point.
+  const auto p = slice(advected_p, cni);
+  // Normalize p for much faster calc_sphere_to_ref.
+  Real p_sphere[3];
+  for (Int d = 0; d < 3; ++d) p_sphere[d] = p[d];
+  siqk::SphereGeometry::normalize(p_sphere);
+  if (calc_ci) {
+    if (m.nonuni) {
+      SearchFunctor sf(m, p_sphere);
+      rd.octree().apply(sf.bb, sf); // same regardless of mesh
+      ci = sf.ci_hit;
+    } else {
+      const auto& gr = m.grid_rotation;
+      ci = mesh::get_cell_idx(ne, gr.angle, gr.R, p_sphere[0], p_sphere[1],
+                              p_sphere[2]);
+    }
+  }
+  assert(ci >= 0 && ci < nslices(m.dglln2cglln)/square(m.np));
+  const auto& cell = slice(m.geo_c2n, ci);
+  siqk::sqr::Info info;
+  sphere2ref::calc_sphere_to_ref(m.geo_p, cell, p_sphere, a, b, &info);
+  // Eval basis functions.
+  m.basis->eval(m.np, a, va);
+  m.basis->eval(m.np, b, vb);
+}
+
+static void
+calc_jacobian_departure (const Mesh& m, const bool isoparametric,
+                         const AVec3s& advected_p, Real* Jdep) {
+  const Int dnn = nslices(m.dglln2cglln);
+  const Int np = m.np, np2 = square(np);
+  const Real* xnode;
+  m.basis->get_x(m.np, xnode);
+  ompparfor for (Int i = 0; i < dnn; ++i) {
+    const Int ti = i / np2;
+    const auto& cell = slice(m.cgll_c2n, ti);
+    const Int ni = i % np2;
+    const Real ta = xnode[ni % np], tb = xnode[ni / np];
+    Real Jd;
+    if (isoparametric) {
+      // This is the case of interest and runs when property preservation is
+      // on. For the case of p-refinement, it runs with np=npv=4, not npt.
+      Jd = calc_isoparametric_jacobian(advected_p, cell, np, ta, tb);
+    } else {
+      // This happens only when property preservation is off, in which case
+      // density doesn't couple to the mixing ratios.
+      const Int corners[] = {cell[0], cell[np-1],
+                             cell[np*np-1], cell[np*(np-1)]};
+      Jd = calc_jacobian(advected_p, corners, ta, tb);
+    }
+    Jdep[i] = Jd;
+  }
+}
+
 void Remapper::
 interp (const Mesh& m, const C2DRelations& c2d, const AVec3s& advected_p,
         const Real* const src_tracer, Real* const tgt_tracer, const Int ntracers,
         const Real* const src_density, Real* const tgt_density,
-        const bool rho_csl) {
+        const bool rho_isl) {
+  assert(src_density && tgt_density);
   Timer::start(Timer::ts_isl_interp);
-  static const GLL gll;
   const Int cnn = nslices(c2d.ptr) - 1;
+  assert(cnn == nslices(advected_p));
   const Int dnn = nslices(m.dglln2cglln);
   const Int ne = std::sqrt(ncell_ / 6);
   const Int np = m.np, np2 = square(np);
   const bool footprint = ntracers > 0 && isl_impl_->footprint != nullptr;
-  if (rho_csl) {
-    const Real* gll_x, * gll_wt;
-    gll.get_coef(m.np, gll_x, gll_wt);
-#   pragma omp parallel for
-    for (Int i = 0; i < dnn; ++i) {
-      const Int ti = i / np2;
-      const auto& cell = slice(m.cgll_c2n, ti);
-      const Int ni = i % np2;
-      const Real ta = gll_x[ni % np], tb = gll_x[ni / np];
-      Real Jd;
-      if (md_) {
-        // This is the case of interest and runs when property preservation is
-        // on. For the case of p-refinement, it runs with np=npv=4, not npt.
-        Jd = calc_isoparametric_jacobian(advected_p, cell, np, ta, tb);
-      } else {
-        // This happens only when property preservation is off, in which case
-        // density doesn't couple to the mixing ratios.
-        const Int corners[] = {cell[0], cell[np-1],
-                               cell[np*np-1], cell[np*(np-1)]};
-        Jd = calc_jacobian(advected_p, corners, ta, tb);
-      }
-      density_factor_(i) = Jd/Je_(i);
-    }
-  }
+  const Real* wt;
+  auto& colscl = isl_impl_->wrk;
+  if (rho_isl)
+    calc_jacobian_departure(m, md_ != nullptr, advected_p, Jdep_.data());
   if (footprint) isl_impl_->footprint->start();
-# pragma omp parallel for
-  for (Int cni = 0; cni < cnn; ++cni) {
+  ompparfor for (Int tni = 0; tni < cnn; ++tni) {
     // GLL values in the source cell corresponding to this departure point.
     Int ci;
     Real va[GLL::np_max], vb[GLL::np_max];
     Real a, b;
-    {
-      // Solve for departure point.
-      const auto p = slice(advected_p, cni);
-      // Normalize p for much faster calc_sphere_to_ref.
-      Real p_sphere[3];
-      for (Int d = 0; d < 3; ++d) p_sphere[d] = p[d];
-      siqk::SphereGeometry::normalize(p_sphere);
-      if (m.nonuni) {
-        SearchFunctor sf(m, p_sphere);
-        rd_->octree().apply(sf.bb, sf); // same regardless of mesh
-        ci = sf.ci_hit;
-      } else {
-        ci = mesh::get_cell_idx(ne, p_sphere[0], p_sphere[1], p_sphere[2]);
+    calc_departure_data(m, *rd_, ne, tni, advected_p, ci, true, va, vb, a, b);
+    cn_src_cell_[tni] = ci;
+    if (rho_isl) {
+      const Real* const src_cell = src_density + ci*np2;
+      Real iv = 0;
+      for (Int j = 0, k = 0; j < np; ++j)
+        for (Int i = 0; i < np; ++i, ++k)
+          iv += va[i]*vb[j]*src_cell[k];
+      for (Int j = c2d.ptr[tni]; j < c2d.ptr[tni+1]; ++j) {
+        const Int di = c2d.r[j];
+        tgt_density[di] = (Jdep_(di)/Je_(di))*iv;
       }
-      const auto& cell = slice(m.geo_c2n, ci);
-      siqk::sqr::Info info;
-      sphere2ref::calc_sphere_to_ref(m.geo_p, cell, p_sphere, a, b, &info);
-      cn_src_cell_[cni] = ci;
-      // Eval basis functions.
-      m.basis->eval(np, a, va);
-      m.basis->eval(np, b, vb);
     }
     for (Int tri = 0; tri < ntracers; ++tri) {
       // Interp the source cell to this point.
       const Real* const src_cell = src_tracer + tri*dnn + ci*np2;
+      Real* const tgt = tgt_tracer + tri*dnn;
       Real iv = 0;
       for (Int j = 0, k = 0; j < np; ++j)
         for (Int i = 0; i < np; ++i, ++k)
           iv += va[i]*vb[j]*src_cell[k];
-      if (isl_impl_->tend->sources[tri])
-        iv += isl_impl_->tend->sources[tri]->interp(ci, a, b);
       // Scatter to DGLL target values corresponding to this CGLL node.
-      Real* const tgt = tgt_tracer + tri*dnn;
-      for (Int j = c2d.ptr[cni]; j < c2d.ptr[cni+1]; ++j)
+      for (Int j = c2d.ptr[tni]; j < c2d.ptr[tni+1]; ++j)
         tgt[c2d.r[j]] = iv;
       if (footprint && tri == 0) {
-        for (Int j = c2d.ptr[cni]; j < c2d.ptr[cni+1]; ++j)
+        for (Int j = c2d.ptr[tni]; j < c2d.ptr[tni+1]; ++j)
           isl_impl_->footprint->record(c2d.r[j], ci);
       }
     }
-    if (rho_csl) {
-      const Real* const src_cell = src_density + ci*np2;
-      Real iv = 0;
-      for (Int j = 0, k = 0; j < np; ++j)
-        for (Int i = 0; i < np; ++i, ++k)
-          iv += va[i]*vb[j]*src_cell[k];
-      for (Int j = c2d.ptr[cni]; j < c2d.ptr[cni+1]; ++j) {
-        const Int di = c2d.r[j];
-        tgt_density[di] = density_factor_(di)*iv;
-      }
-    }
   }
   if (footprint) isl_impl_->footprint->finish();
   Timer::stop(Timer::ts_isl_interp);
 }
 
 void Remapper::
-csl_cdr_rho (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
+isl_cdr_rho (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
              const MonoData& md, const spf::MassRedistributor::Ptr& mrd,
              Real* const src_density, const Int np2_src,
              Real* const tgt_density, const Int np2_tgt,
@@ -1296,11 +1290,11 @@ csl_cdr_rho (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
 
 // Constrained density reconstruction for classical semi-Lagrangian.
 void Remapper::
-csl_cdr (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
+isl_cdr (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
          const MonoData& md, const spf::MassRedistributor::Ptr& mrd,
          Real* const src_density, Real* const src_tracer, const Int np_src,
          Real* const tgt_density, Real* const tgt_tracer, const Int ntracers,
-         const bool positive_only, const bool rho_csl, const Filter::Enum cdr_method) {
+         const bool positive_only, const bool rho_isl, const Filter::Enum cdr_method) {
   const auto& td = *isl_impl_->tend;
   const Int np_tgt = m.np, np2_tgt = square(np_tgt);
   const Int np2_src = square(np_src);
@@ -1312,8 +1306,8 @@ csl_cdr (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
     total_mass_redistribution_.resize(ntracers+1);
     total_mass_discrepancy_.resize(ntracers+1);
   }
-  if (rho_csl) {
-    csl_cdr_rho(m, rd_src, rd_tgt, md, mrd,
+  if (rho_isl) {
+    isl_cdr_rho(m, rd_src, rd_tgt, md, mrd,
                 src_density, np2_src,
                 tgt_density, np2_tgt, ntracers);
     assert(nslices(m.dglln2cglln) == d2cer_->get_dnn());
@@ -1329,8 +1323,6 @@ csl_cdr (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
     accumulate_threaded_bfb<1>(
       [&] (const Int i, Real* a) { *a += F_src[i]*src[i]*src_density[i]; },
       len_src, &Q_mass_src);
-    if (isl_impl_->tend->sources[tri])
-      Q_mass_src += isl_impl_->tend->sources[tri]->global_mass();
     accumulate_threaded_bfb<1>(
       [&] (const Int i, Real* a) {
         tgt[i] *= tgt_density[i];
@@ -1356,14 +1348,9 @@ csl_cdr (const Mesh& m, const RemapData& rd_src, const RemapData& rd_tgt,
         }
         if ( ! td.cdr_on_src_tends &&
             (td.dss_source || (tri >= td.idx_beg && tri < td.idx_end))) {
-          if (td.sources[tri]) {
-            q_min = td.q_min[ncell_*tri + ci];
-            q_max = td.q_max[ncell_*tri + ci];
-          } else {
-            // Get the tightest bounds we can.
-            q_min = std::max(q_min, td.q_min[ncell_*tri + ci]);
-            q_max = std::min(q_max, td.q_max[ncell_*tri + ci]);
-          }
+          // Get the tightest bounds we can.
+          q_min = std::max(q_min, td.q_min[ncell_*tri + ci]);
+          q_max = std::min(q_max, td.q_max[ncell_*tri + ci]);
         }
         if (fit_extremum_) {
           assert(np_tgt == np_); // for now
@@ -1488,10 +1475,10 @@ ::init_physgrid (const Real* const rho, const Real* const tracer, const Int ntra
 }
 
 void Remapper
-::csl (const AVec3s& advected_p, Real* const src_tracer,
+::isl (const AVec3s& advected_p, Real* const src_tracer,
        Real* const tgt_tracer, const Int ntracers,
        Real* const src_density, Real* const tgt_density,
-       const bool positive_only, const bool rho_csl,
+       const bool positive_only, const bool rho_isl,
        const Real* srcterm_tendency /* mixing ratio */,
        const Int srcterm_idx_beg, const Int srcterm_idx_end,
        const pg::PhysgridData::Ptr& pg) {
@@ -1526,12 +1513,12 @@ ::csl (const AVec3s& advected_p, Real* const src_tracer,
     // original basic p-refinement.
     interp(*m_, c2d_v_, advected_p,
            src_tracer, tgt_tracer, ntracers,
-           src_density, tgt_density, rho_csl);
+           src_density, tgt_density, rho_isl);
     if (run_cdr)
-      csl_cdr(*m_, *rd_, *rd_, *md_, mrd_,
+      isl_cdr(*m_, *rd_, *rd_, *md_, mrd_,
               src_density, src_tracer, np_,
               tgt_density, tgt_tracer, ntracers,
-              positive_only, rho_csl, cdr_method);
+              positive_only, rho_isl, cdr_method);
     else
       dss_rho(*d2cer_, tgt_density, isl_impl_->wrk);
     if (pg)
@@ -1558,17 +1545,17 @@ ::csl (const AVec3s& advected_p, Real* const src_tracer,
       }
       d2cer_ = std::make_shared<D2Cer>(pr_->m_v->dglln2cglln, pr_->rd_v->dgbfi_mass());
     }
-    if (rho_csl) {
+    if (rho_isl) {
       // Transport rho on v mesh to mimic dycore.
       interp(*pr_->m_v, c2d_v_, advected_p,
              nullptr, nullptr, 0,
-             src_rho_impl, tgt_rho_impl, rho_csl);
+             src_rho_impl, tgt_rho_impl, rho_isl);
       if (run_cdr) {
         // Property-preserve rho on v mesh to mimic dycore.
-        csl_cdr(*pr_->m_v, *pr_->rd_v, *pr_->rd_v, *pr_->md_v, isl_impl_->mrd_c,
+        isl_cdr(*pr_->m_v, *pr_->rd_v, *pr_->rd_v, *pr_->md_v, isl_impl_->mrd_c,
                 src_rho_impl, nullptr, pr_->m_v->np,
                 tgt_rho_impl, nullptr, ntracers,
-                positive_only, rho_csl, cdr_method);
+                positive_only, rho_isl, cdr_method);
       } else {
         dss_rho(*d2cer_, tgt_rho_impl, isl_impl_->wrk);
       }
@@ -1579,9 +1566,11 @@ ::csl (const AVec3s& advected_p, Real* const src_tracer,
     const auto ap_fine = isl_impl_->interp_p(advected_p);
     interp(*m_, c2d_f_, ap_fine,
            src_tracer, tgt_tracer, ntracers,
-           nullptr, nullptr, false);
+           // These are unused unless we're using the local mass conservation
+           // option.
+           src_rho_impl, tgt_rho_impl, false);
     if (run_cdr)
-      csl_cdr(*m_, *rd_, *rd_, *md_, mrd_,
+      isl_cdr(*m_, *rd_, *rd_, *md_, mrd_,
               src_density, src_tracer, np_,
               tgt_density, tgt_tracer, ntracers,
               positive_only, false, cdr_method);
@@ -1593,17 +1582,17 @@ ::csl (const AVec3s& advected_p, Real* const src_tracer,
     Real* const tgt_rho_impl = isl_impl_->rho[isl_impl_->idxs[1]].data();
     Real* const src_tracer_impl = isl_impl_->tracer[isl_impl_->idxs[0]].data();
     Real* const tgt_tracer_impl = isl_impl_->tracer[isl_impl_->idxs[1]].data();
-    if (rho_csl) {
+    if (rho_isl) {
       // Transport rho on v mesh to mimic dycore.
       interp(*pr_->m_v, c2d_v_, advected_p,
              nullptr, nullptr, 0,
-             src_density, tgt_density, rho_csl);
+             src_density, tgt_density, rho_isl);
       if (run_cdr) {
         // Property-preserve rho on v mesh to mimic dycore.
-        csl_cdr(*pr_->m_v, *pr_->rd_v, *pr_->rd_v, *pr_->md_v, mrd_,
+        isl_cdr(*pr_->m_v, *pr_->rd_v, *pr_->rd_v, *pr_->md_v, mrd_,
                 src_density, nullptr, pr_->m_v->np,
                 tgt_density, nullptr, ntracers,
-                positive_only, rho_csl, cdr_method);
+                positive_only, rho_isl, cdr_method);
       } else {
         dss_rho(*d2cer_, tgt_density, isl_impl_->wrk);
       }
@@ -1633,13 +1622,15 @@ ::csl (const AVec3s& advected_p, Real* const src_tracer,
     const auto ap_fine = isl_impl_->interp_p(advected_p);
     interp(*pr_->m_f, c2d_f_, ap_fine,
            src_tracer_impl, tgt_tracer_impl, ntracers,
-           nullptr, nullptr, false);
+           // These are unused unless we're using the local mass conservation
+           // option.
+           src_rho_impl, tgt_rho_impl, false);
     // Interpolate current rho to f mesh.
     isl_impl_->interp_rho(tgt_density, tgt_rho_impl, rho_continuous);
     // Global property preservation on t mesh, local on both.
     if (run_cdr) {
       // Property preserve q on f mesh.
-      csl_cdr(*pr_->m_f, *pr_->rd_f, *pr_->rd_f, *pr_->md_f, isl_impl_->mrd_f,
+      isl_cdr(*pr_->m_f, *pr_->rd_f, *pr_->rd_f, *pr_->md_f, isl_impl_->mrd_f,
               src_rho_impl, src_tracer_impl, pr_->m_f->np,
               tgt_rho_impl, tgt_tracer_impl, ntracers,
               positive_only, false /* don't apply cdr to rho */, cdr_method);
diff --git a/methods/slmm/slmmir_test.cpp b/methods/slmm/slmmir_test.cpp
index 53f0ddf..7d27d97 100644
--- a/methods/slmm/slmmir_test.cpp
+++ b/methods/slmm/slmmir_test.cpp
@@ -1,59 +1,33 @@
 #include <stdexcept>
 
-#include "slmm_defs.hpp"
-#include "slmm_mesh.hpp"
-#include "slmm_debug.hpp"
-using namespace slmm;
+#include "slmmir.hpp"
 
 struct Command {
-  enum Enum { test_make_cubedsphere };
+  enum Enum { test_csl_density };
   static Enum from_string (const std::string& s) {
-    if (s == "test_make_cubedsphere") return test_make_cubedsphere;
+    if (s == "test_csl_density") return test_csl_density;
     throw std::runtime_error(s + " is not a command.");
   }
 };
 
 struct Input {
   Command::Enum command;
-  Int n;
-  Real angle;
-  bool write_matlab;
 
   Input(Int argc, char** argv);
   void print(std::ostream& os) const;
 };
 
-static Int test_make_cubedsphere (const Input& in) {
-  AVec3s cp;
-  AIdxs ce;
-  mesh::make_cubedsphere(cp, ce, in.n);
+static Int test_csl_density (const Input& in) {
   Int nerr = 0;
-  {
-    const Int ne = mesh::check_elem_normal_against_sphere(cp, ce);
-    if (ne) std::cerr << "FAIL: check_elem_normal_against_sphere\n";
-    nerr += ne;
-  }
-  if (in.write_matlab)
-    write_matlab("cm", cp, ce);
   return nerr;
 }
 
-inline bool
-eq (const std::string& a, const char* const b1, const char* const b2 = 0) {
-  return (a == std::string(b1) || (b2 && a == std::string(b2)) ||
-          a == std::string("-") + std::string(b1));
-}
-
 Input::Input (Int argc, char** argv)
-  : command(Command::test_make_cubedsphere), n(10), angle(M_PI*1e-1),
-    write_matlab(false)
+  : command(Command::test_csl_density)
 {
   for (Int i = 1; i < argc; ++i) {
     const std::string& token = argv[i];
     if (eq(token, "-c", "--comand")) command = Command::from_string(argv[++i]);
-    else if (eq(token, "-n")) n = atoi(argv[++i]);
-    else if (eq(token, "-m", "--write-matlab")) write_matlab = true;
-    else if (eq(token, "--angle")) angle = atof(argv[++i]);
   }
 
   print(std::cout);
@@ -61,9 +35,7 @@ Input::Input (Int argc, char** argv)
 
 void Input::print (std::ostream& os) const {
   os << "command " << command << "\n"
-     << "n (-n): " << n << "\n"
-     << "write matlab (-m): " << write_matlab << "\n"
-     << "angle (--angle): " << angle << "\n";
+     << "\n";
 }
 
 int main (int argc, char** argv) {
@@ -72,8 +44,8 @@ int main (int argc, char** argv) {
     Input in(argc, argv);
     Int nerr = 0;
     switch (in.command) {
-    case Command::test_make_cubedsphere:
-      nerr = test_make_cubedsphere(in);
+    case Command::test_csl_density:
+      nerr = test_csl_density(in);
       break;
     }
     std::cerr << (nerr ? "FAIL" : "PASS") << "ED\n";
diff --git a/methods/slmm/slmmir_time_int.cpp b/methods/slmm/slmmir_time_int.cpp
index f3dfc9a..f43008e 100644
--- a/methods/slmm/slmmir_time_int.cpp
+++ b/methods/slmm/slmmir_time_int.cpp
@@ -116,6 +116,9 @@ ::create (const Enum& ode, const bool use_xyz_form,
   case MovingVortices:
     return std::make_shared<MeshIntegratorWithOdeFn<
       gallery::MovingVortices> >(p, use_xyz_form);
+  case TestWindField:
+    return std::make_shared<MeshIntegratorWithOdeFn<
+      gallery::TestWindField> >(p, use_xyz_form);
   default:
     assert(0);
     return nullptr;
diff --git a/siqk/siqk_geometry.hpp b/siqk/siqk_geometry.hpp
index 9a19d6b..b3fba94 100644
--- a/siqk/siqk_geometry.hpp
+++ b/siqk/siqk_geometry.hpp
@@ -122,8 +122,8 @@ struct SphereGeometry {
     c[1] = a[2]*b[0] - a[0]*b[2];
     c[2] = a[0]*b[1] - a[1]*b[0];
   }
-  template <typename CV> KOKKOS_INLINE_FUNCTION
-  static Real dot (const CV a, const CV b) {
+  template <typename CVA, typename CVB> KOKKOS_INLINE_FUNCTION
+  static Real dot (const CVA a, const CVB b) {
     return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
   }
   template <typename CV> KOKKOS_INLINE_FUNCTION