diff --git a/methods/islet/figures/figs-adv-diag.hy b/methods/islet/figures/figs-adv-diag.hy
index 9ceab27..fb25f8e 100644
--- a/methods/islet/figures/figs-adv-diag.hy
+++ b/methods/islet/figures/figs-adv-diag.hy
@@ -1,61 +1,13 @@
 (require [amb3 [*]])
 (import amb3 [amb3 [*]]
         [figsutils [*]]
-        math glob struct os)
+        math glob struct os matplotlib.ticker)
 
 (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]
-  (sv txt (.split (readall fname) "\n")
-      bo None d {})
-  (for [ln txt]
-    (sv ln2 (cut ln 0 2))
-    (cond [(and (= "cm" ln2) (in "cmd>" ln))
-           (sv cmd ln c (parse-cmd ln :map-nstepfac map-nstepfac))]
-          [(= ln2 "M ")
-           (sv mp (parse-midpoint-check ln))
-           (unless (none? (geton d #*(+ (cmd->key-base c)
-                                        (, (:ic mp) cyc (:ne c) :M))))
-             ;; handle repeated ic used for src-term OOA testing.
-             (assoc mp :ic (+ (:ic mp) "-src")))
-           (assoc-nested d (+ (cmd->key-base c) (, (:ic mp) cyc (:ne c) :M))
-                         {:l1 (:l1 mp) :l2 (:l2 mp)})]
-          [(= ln2 "C ")
-           (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") (!= (:mono c) "none") (<= cyc 10))
-                    (prf "FAIL {}" cmd))]
-                 [:else
-                  (sv cl (parse-C ln))
-                  (unless (none? (geton d #*(+ (cmd->key-base c)
-                                               (, (:ic cl) cyc (:ne c) :C))))
-                    ;; handle repeated ic used for src-term OOA testing.
-                    (assoc cl :ic (+ (:ic cl) "-src")))
-                  (assoc-nested d (+ (cmd->key-base c) (, (:ic cl) cyc (:ne c) :C))
-                                cl)])]
-          [(= ln2 "L ")
-           (cond [(in "L file" ln)
-                  (assoc-nested d (+ (cmd->key-base c)
-                                     (, "cos" cyc (:ne c) :L :mixing-file))
-                                (last (.split ln)))]
-                 [(in "phimin" ln)
-                  (sv (, - ic - l1 - l2 - li - phimin - phimax)
-                      (sscanf ln "s,s,s,f,s,f,s,f,s,f,s,f"))
-                  (assoc-nested d (+ (cmd->key-base c) (, ic cyc (:ne c) :Lerr))
-                                {:l1 l1 :l2 l2 :li li :phimin phimin :phimax phimax})]
-                 [:else
-                  (sv bo (parse-bakeoff-diag bo ln (:timeint c)))
-                  (when (:done bo)
-                    (assoc-nested d (+ (cmd->key-base c) (, "cos" cyc (:ne c) :L)) bo)
-                    (sv bo None))])]))
-  d)
-
 ;;; print a long table of all accuracy results
 
 (defn acc-print-txt-table [c d &optional]
@@ -194,6 +146,15 @@
 
 ;;; accuracy figs
 
+(defn show-logy-ticks [&optional ax]
+  (svifn ax (pl.gca))
+  (sv ma (matplotlib.ticker.LogLocator :base 10, :numticks 100))
+  (ax.yaxis.set-major-locator ma)
+  (sv mi (matplotlib.ticker.LogLocator :base 10 :subs (npy.linspace 0.1 1 10)
+                                       :numticks 100))
+  (ax.yaxis.set-minor-locator mi)
+  (ax.yaxis.set-minor-formatter (matplotlib.ticker.NullFormatter)))
+
 (defclass AccFig []
   (defn --init-- [me])
 
@@ -205,7 +166,7 @@
      :nps c.nps :nes c.nes :pat-line "-" :pat-clr (.copy c.npclrs) :C-line :C
      :pat-mark (.copy c.npmarks) :fs 11 :lw 1 :markersize 4 :yform :log2
      :xticks :degrees :ooa-text False :filter-floor None :figsize (, 4 4)
-     :ref-ooa-2 False :ref-cos-033 True :pg None})
+     :ref-ooa-2 False :ref-cos-033 True :pg None :jcp False})
 
   (defn plot [me ax d-all &optional [o None]]
     (svifn o (me.get-defaults))
@@ -215,6 +176,7 @@
         gray (* 0.5 (npy.ones 3))
         y-ext [1e10 -1e10])
     (when (none? d1) (return))
+    (sv x-min 1000 x-max 0)
     (for [np (:nps o)]
       (sv keys (if (none? (:pg o))
                    (, np (:ic o) (:cyc o))
@@ -224,12 +186,20 @@
           d2 (geton d1 #*keys))
       (when (none? d2) (continue))
       (sv x [] y [])
-      (for [ne (:nes o)]
-        (sv val (geton d2 ne (:C-line o) (:measure o)))
+      (for [tne (:nes o)]
+        (sv ne-key (if (:jcp o)
+                     (jcp-tne2ne tne np)
+                     tne)
+            tne-actual (if (:jcp o)
+                         (jcp-effective-tne tne np)
+                         tne)
+            val (geton d2 ne-key (:C-line o) (:measure o)))
         (when (none? val) (continue))
         (unless (none? (:filter-floor o))
           (when (< val (:filter-floor o)) (continue)))
-        (.append x ne)
+        (sv x-min (min x-min tne-actual)
+            x-max (max x-max tne-actual))
+        (.append x tne-actual)
         (.append y val)
         (sv (get y-ext 0) (min (get y-ext 0) val)
             (get y-ext 1) (max (get y-ext 1) val)))
@@ -258,21 +228,24 @@
                    :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))) "-."
-                   :zorder -10 :lw (:lw o) :color gray)))
-      (when (:ooa-text o)
+               (pl.yticks (list (range -40 10)) :fontsize (:fs o))]))
+      (when (and (:ooa-text o) (> (len x) 1))
         (sv i (- (len x) 2))
-        (pl.text (+ (get xtform (inc i)) -0.08) (* (get ytform (inc i)) 2)
+        (pl.text (+ (get xtform (inc i)) -0.08)
+                 (* (get ytform (inc i)) (if (= (:measure o) :l2) 0.25 2))
                  (.format "{:1.1f}"
                           (first (calc-ooa (cut y i (+ i 2))
                                            :x (cut x i (+ i 2)))))))
       (when (= np (last (:nps o)))
+        (when (and (:ref-cos-033 o) (= (:ic o) "cos") (= (:ode o) "nondivergent"))
+          (pl-plot (npy.log2 (npa [x-min x-max]))
+                   (* 0.033 (npy.ones 2)) "-."
+                   :zorder -10 :lw (:lw o) :color gray))
         (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)
+          (show-logy-ticks)
           (sv (, y ys) (pl.yticks))
           (when (> (/ (last y) (first y)) 1e3)
             (for [i (range (len ys))]
@@ -504,6 +477,139 @@
             (pl.ylim (, (/ (first ye) 4) (nextpow10 (second ye)))))
           :ref-ooa-2 True)))
 
+(defn figs-acc-jcp [c d]
+  (sv prefix ""
+      ref-ooa-2 False
+      legend True
+      general-timeint "exact"
+      general-prefine 0
+      pp False
+      jcp True)
+  (sv p (AccFig)
+      o (p.get-defaults c))
+  (assoc o :ode "nondivergent" :ic "gau" :nstepfac 5
+         :yform :semilogy :cdrglb "none" :cdrlcl "none"
+         :timeint general-timeint :prefine general-prefine
+         :filter-floor None :jcp True :ref-cos-033 False)
+  (defn plot [spi o plot-fn &optional [ref-ooa-2 False] xlabel]
+    (svifn xlabel "Resolution")
+    (sv fname (+ prefix "acc-" (:ode o) "-" (:ic o)
+                 "-" (:timeint o) "-" (if (= (:cdrglb o) "none")
+                                          "nopp" "pp")
+                 "-fac" (str (:nstepfac o))
+                 "-jcp"))
+    (print fname)
+    (do
+      (sv ax (pl.subplot 1 2 (inc spi)))
+      (plot-fn ax d o)
+      (my-grid)
+      (sv legs [(, "k-" "$l_2$")])
+      (when (= ic "gau")
+        (.append legs (, "k--" "$l_{\infty}$")))
+      (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)
+      (sv flowname (flow-short2long (:ode o))
+          flowname (+ (.upper (first flowname)) (cut flowname 1))
+          f (pl.gcf))
+      (f.text (+ (* 0.5 spi) 0.05) 0.05
+              (if (zero? spi) "(a)" "(b)") :fontsize (:fs o))
+      (p.title (+ flowname ", "
+                  (ic-short2long (cut (:ic o) 0 3)) ", "
+                  (nstepfac->word (:nstepfac o)) " steps") o)))
+  (for [nstepfac (, 1 5)]
+    (with [(pl-plot (, 8 4) (+ c.fig-dir "acc-jcp-nstepfac" (str nstepfac)))]
+      (for [(, iic ic) (enumerate (, "gau" "cos"))
+            ode (, "nondivergent")]
+        (assoc o :nstepfac nstepfac :ooa-text None :ode ode
+               :ic ic)
+        (plot iic o
+              (fn [ax d o]
+                (sv ye [1e10 -1e10])
+                (defn update-ye [ye1]
+                  (when (none? ye1) (return))
+                  (sv (get ye 0) (min (first ye) (first ye1))
+                      (get ye 1) (max (last ye) (last ye1))))
+                (for [norm (, :l2 :li)]
+                  (when (and (= norm :li) (= ic "cos")) (continue))
+                  (assoc o :cdrglb (if pp "caas" "none") :timeint "exact" :prefine 0
+                         :measure norm :pat-line (if (= norm :l2) "-" "--")
+                         :ref-ooa-2 (and (!= ic "gau") (= norm :l2))
+                         :ooa-text (and (= ic "gau") (= norm :l2) (= nstepfac 1)))
+                  (sv ye1 (p.plot ax d o))
+                  (update-ye ye1))
+                (case/eq ic
+                         ["cos"
+                          (pl.ylim (, (/ (first ye) 3) 0.4))]
+                         ["gau"
+                          (pl.ylim (, (/ (first ye) 60) 0.5))]))
+              :ref-ooa-2 (and ref-ooa-2 (!= ic "slo")))))))
+
+(defn figs-vortex-acc-jcp [c d]
+  (sv prefix ""
+      ref-ooa-2 False
+      legend True
+      general-timeint "exact"
+      general-prefine 0
+      show-linf True
+      pp False
+      jcp True)
+  (sv p (AccFig)
+      o (p.get-defaults c))
+  (assoc o :ode "movingvortices" :ic "vor"
+         :yform :semilogy :cdrglb "none" :cdrlcl "none"
+         :timeint general-timeint :prefine general-prefine
+         :filter-floor None :jcp True)
+  (defn plot [spi o plot-fn &optional [ref-ooa-2 False] xlabel]
+    (svifn xlabel "Resolution")
+    (sv fname (+ prefix "acc-" (:ode o) "-" (:ic o)
+                 "-" (:timeint o) "-" (if (= (:cdrglb o) "none")
+                                          "nopp" "pp")
+                 "-fac" (str (:nstepfac o))
+                 "-jcp"))
+    (print fname)
+    (do
+      (sv ax (pl.subplot 1 2 (inc spi)))
+      (plot-fn ax d o)
+      (my-grid)
+      (sv legs [(, "k-" "$l_2$")])
+      (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)
+      (sv flowname (flow-short2long (:ode o))
+          flowname (+ (.upper (first flowname)) (cut flowname 1))
+          f (pl.gcf))
+      (f.text (+ (* 0.5 spi) 0.05) 0.05
+              (if (zero? spi) "(a)" "(b)") :fontsize (:fs o))
+      (p.title (+ flowname ", "
+                  (nstepfac->word (:nstepfac o)) " steps, "
+                  "cycle " (str (:cyc o))) o)))
+  (for [nstepfac (, 1 5)]
+    (with [(pl-plot (, 8 4)
+                    (+ c.fig-dir "vortex-acc-jcp-nstepfac" (str nstepfac)))]
+      (for [(, icyc cycle) (enumerate (, 1 2))
+            ode (, "moving vortices")]
+        (assoc o :nstepfac nstepfac :ooa-text None :cyc cycle)
+        (plot icyc o
+              (fn [ax d o]
+                (sv ye [1e10 -1e10])
+                (defn update-ye [ye1]
+                  (when (none? ye1) (return))
+                  (sv (get ye 0) (min (first ye) (first ye1))
+                      (get ye 1) (max (last ye) (last ye1))))
+                (for [norm (, :l2)]
+                  (when (and (= norm :li) (not show-linf)) (continue))
+                  (assoc o :cdrglb (if pp "caas" "none") :timeint "exact" :prefine 0
+                         :measure norm :pat-line (if (= norm :l2) "-" "--")
+                         :ref-ooa-2 False :ooa-text (and (= norm :l2) (zero? icyc)))
+                  (sv ye1 (p.plot ax d o))
+                  (update-ye ye1))
+                (pl.ylim (, (/ (first ye) (if (zero? icyc) 60 6))
+                            (* (second ye) 1.5))))
+              :ref-ooa-2 False)))))
+
 ;;; filament diagnostic
 
 (defn fig-filament [c d-all]
@@ -550,16 +656,58 @@
                  (.format (if 0 "$n_p$ {}{}" "$n_p$ {} {}") np extra)
                  :fontsize (+ (:fs o) 2))))))
 
+(defn fig-filament-jcp [c d-all]
+  (defn get-pat [nstepfac ne]
+    (+ (get {20 "r." 40 "k"} ne)
+       (get {1 "-" 5 "--"} nstepfac)))
+  (sv p (AccFig) o (p.get-defaults c)
+      tnes (, 20 40)
+      nc (len c.nps))
+  (with [(pl-plot (, 10 4) (+ c.fig-dir "filament-jcp") :tight False)]
+    (for [(, itne tne) (enumerate tnes)]
+      (for [(, inp np) (enumerate c.nps)]
+        (sv spi (inc (+ (* itne (len c.nps)) inp))
+            ;;ax (pl.subplot 2 (len c.nps) spi)
+            ax (pl.axes [(/ inp nc) (/ (% (inc itne) 2) 2)
+                         (* 0.95 (/ 1 nc)) (* 0.95 (/ 1 2))])
+            ne (jcp-tne2ne tne np))
+        (for [nstepfac (, 1 5)]
+          (sv d (get d-all "exact" "nondivergent" nstepfac
+                     "pcsl" "none" "none" 0
+                     np "cos" 1))
+          (when (none? (geton d ne)) (continue))
+          (sv thr (get d ne :L :thr)
+              fil (get d ne :L :fil))
+          (pl.plot thr fil (get-pat nstepfac tne)
+                   :label (.format "{} step" (nstepfac->word nstepfac)))
+          (pl.xticks [0 0.2 0.4 0.6 0.8 1] :fontsize (:fs o))
+          (pl.xlim (, 0.05 1.05))
+          (pl.ylim (, -5 125))
+          (sv yl (pl.ylim))
+          (my-grid)
+          (when (= spi (* 2 (len c.nps)))
+            (pl.legend :loc "center" :fontsize (:fs o) :frameon False))
+          (sv yt [0 20 40 60 80 100 120]
+              xt [0.2 0.4 0.6 0.8 1.0])
+          (if (zero? inp)
+            (pl.yticks yt)
+            (pl.yticks yt (* [""] (len yt))))
+          (if (zero? itne)
+            (pl.xticks xt (* [""] (len xt)))
+            (pl.xticks xt)))
+        (pl.text 0.1 (+ 0 (* 0.05 (npy.diff yl)))
+                 (.format "$n_e \, {}$, $n_p \, {}$" ne np)
+                 :fontsize (+ (:fs o) 2))))))
+
 ;;; mixing diagnostic
 
 (defn triplot-read-dat [fname]
-    (sv raw (with [f (open fname "rb")]
-            (.read f))
+  (sv raw (with [f (open fname "rb")] (.read f))
       n (first (struct.unpack "i" (get raw (slice 0 4))))
       data (npy.zeros (* 2 n)))
   (for [i (range (* 2 n))]
     (sv os (+ 4 (* 4 i))
-        (get data i) (first (struct.unpack "f" (cut raw os (+ os 4))))))
+        (get data i) (first (struct.unpack "f" (get raw (slice os (+ os 4)))))))
   (sv cb (cut data 0 n)
       ccb (cut data n (* 2 n)))
   (, cb ccb))
@@ -597,6 +745,38 @@
     (text x 0.32 "$l_r$" 'lr 'me-lr)
     (text x 0.11 "$l_u$" 'lu 'me-lu)))
 
+(defn triplot-jcp [cb ccb &optional [data None]]
+  (defn triplot-curve [x]
+    (+ 0.9 (* -0.8 x x)))
+  (sv x (npy.linspace 0.1 1 100)
+      lw 2)
+  (pl.plot x (triplot-curve x) "k-" :lw lw)
+  (sv tl (triplot-curve (last x))
+      br (triplot-curve (first x)))
+  (pl.plot x (* (triplot-curve (first x)) (npy.ones (len x))) "k-" :lw lw)
+  (pl.plot (npy.ones (len x)) (npy.linspace tl br (len x)) "k-" :lw lw)
+  (pl.plot x (npy.linspace br tl (len x)) "k-" :lw lw)
+  (pl.plot cb ccb "r." :markersize 1)
+  (pl.xlim -0.02 1.05)
+  (pl.ylim  0.02 0.97)
+  (sv t [0.0 0.2 0.4 0.6 0.8 1])
+  (pl.xticks t []) (pl.yticks t [])
+  (my-grid)
+  ;(pl.axis "off")
+  (defn text [x y txt sym]
+    (pl.text x (+ y 0.1) txt)
+    (sv dx 0.1)
+    (pl.text (+ x dx) (+ y 0.1) (.format "{:1.2e}" (get data sym))))
+  (unless (none? data)
+    (sv x 0.03)
+    (pl.text x 0.04 (.format "$n_e \, {}$, $n_p \, {}$, {} step"
+                             (get data 'ne) (get data 'np)
+                             (nstepfac->word (get data 'nstepfac)))
+             :fontsize 11)
+    (text x 0.32 "$l_r$" 'lr)
+    (text x 0.21 "$l_u$" 'lu)
+    (text x 0.1 "$l_o$" 'lo)))
+
 (defn figs-mixing [c d]
   (sv p (AccFig) o (p.get-defaults c))
   (for [ne (, 20 40)]
@@ -618,6 +798,32 @@
                (get e me-mixing :lo)))
         (triplot cb ccb :data data)))))
 
+(defn figs-mixing-jcp [c d]
+  (sv p (AccFig) o (p.get-defaults c)
+      nc (len (:nps o)))
+  (for [tne (, 20 40)]
+    (with [(pl-plot (, 8.3 3.5)
+                    (+ c.fig-dir "mixing-tne" (str tne) "-jcp")
+                    :format "png" :tight False)]
+      (for [(, instepfac nstepfac) (enumerate (, 1 5))
+            (, inp np) (enumerate (:nps o))]
+        (sv spi (inc (+ (* instepfac (len (:nps o))) inp))
+            ax (pl.axes [(/ inp nc) (/ (% (inc instepfac) 2) 2)
+                         (* 0.95 (/ 1 nc)) (* 0.95 (/ 1 2))])
+            ne (jcp-tne2ne tne np)
+            - (print np ne)
+            e (get d "exact" "nondivergent" nstepfac "pcsl"
+                   "none" "none" 0 np
+                   "cos" 1 ne :L)
+            (, cb ccb) (triplot-read-dat (+ c.data-dir (:mixing-file e)))
+            data {'ne ne 'np np 'nstepfac nstepfac
+                  'lr (get e :mixing :lr)
+                  'lu (get e :mixing :lu)
+                  'lo (get e :mixing :lo)})
+        (unless (and (zero? (get e :mixing :lo)) (zero? (get e me-mixing :lo)))
+          (prf "{}: lo {}" (:mixing-file e) (get e :mixing :lo)))
+        (triplot-jcp cb ccb :data data)))))
+
 ;;; slotted cylinders images
 
 (defn img-slo-filament [c d direc img-idx outname &optional nps nps-right]
@@ -1025,7 +1231,7 @@
     (sv d.data (stabgrind-parse s ne :d d.data)))
   d)
 
-(defn stabgrind-fig [c ds outname &optional ic]
+(defn stabgrind-fig [c ds outname &optional ic map-ne-equivalent]
   (svifn ic "gau")
   (sv nes (, 5 10 20 40 60 80)
       yl (, 1e-11 1)
@@ -1033,8 +1239,8 @@
       ;;fsz (, 4 8) sbs (, 2 1)
       fs 12
       np (. (first ds) s np))
-  (defn xticks-ne []
-    (sv xa (nes->degstrs nes))
+  (defn xticks-ne [nes]
+    (sv xa (nes->degstrs nes :convert-all True))
     (pl.xticks (npy.log2 (:ne xa)) (:deg xa) :fontsize fs))
   (defn xticks-cnt []
     (sv xa (npy.array (list (, 1 10 100))))
@@ -1051,18 +1257,21 @@
   (for [d ds]
     (assert (= d.type 'stabgrind-read))
     (sv s d.s keys (stabgrind-keys s)
-        x [] cls [])
+        x [] cls [] ne-keys [])
     (for [ne nes]
       (sv cl (geton d.data #*(+ keys (, ic ne))))
       (when (none? cl) (continue))
-      (.append x ne)
+      (.append ne-keys ne)
+      (.append x (if (none? map-ne-equivalent)
+                   ne
+                   (map-ne-equivalent ne np)))
       (.append cls cl))
-    (sv d.nes x d.cls cls))
+    (sv d.nes x d.cls cls d.ne-keys ne-keys))
   (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))))))
+          ncyc (len (get d.data #*(+ keys (, ic (first d.ne-keys))))))
       (for [cyc (, 1 10 100)]
         (when (> cyc ncyc) (break))
         (sv y [])
@@ -1072,6 +1281,7 @@
         (pl.semilogy (npy.log2 d.nes) y
                      (make-pat s.basis s.nstepfac :marker True)
                      :fillstyle "none")))
+    (show-logy-ticks)
     (do
      (sv xl (pl.xlim))
      (for [basis (, "Gll" "GllNodal") nstepfac (, 1 5)]
@@ -1084,8 +1294,8 @@
     (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)
+    (xticks-ne d.nes) (yticks) (pl.ylim yl) (my-grid)
+    (pl.xlabel "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)
@@ -1096,6 +1306,7 @@
         (sv y [])
         (for [cl cls] (.append y (:l2 cl)))
         (pl.loglog (range 1 (inc (len y))) y (make-pat s.basis s.nstepfac))))
+    (show-logy-ticks)
     (f.text 0.52 0.05 "(b)" :fontsize fs)
     (pl.xlabel "Cycle" :fontsize fs)
     (pl.ylabel "log$_{10}$ $l_2$ relative error" :fontsize fs)
@@ -1226,9 +1437,10 @@
   (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)))))
+                   :map-nstepfac
+                   (if 0
+                     (fn [n &rest r] (get {0.5 1 1 1 2.5 5 5 5} n))
+                     (fn [n &rest r] (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))
@@ -1244,19 +1456,56 @@
                      (first p) ".bin")))
   (img-instab c pairs (+ c.fig-dir "img-instab")))
 
-(when-inp ["fig-stabgrind"]
+(when-inp ["resolution-pairs"]
+  ;; Good np based on (1) matching the np=4 number of grid points (* ne 3) and
+  ;; (2) OOA is a jump over previous np: 4, 6, 11, 13. We also want to include 8
+  ;; or 9, probably 8.
+  (for [ne (, 20 40)]
+    (prf ">>> {}" ne)
+    (sv res-tgt (* ne 3))
+    (for [np (range 4 14)]
+      (sv ene (math.floor (/ res-tgt (dec np))))
+      (prf "{:2d}: {:2d} {:3d}" np ene (* (dec np) ene)))))
+
+(when-inp ["fig-stabgrind-jcp"]
   (sv c (get-context)
       ode "nondivergent"
       np 11
       nes (, 5 10 20 40 80)
       ds [])
   (for [basis (, "Gll" "GllNodal")
-        nstepfac (, 1 5)
-        ]
+        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")))
+  (stabgrind-fig c ds (+ c.fig-dir "fig-stabgrind")
+                 :map-ne-equivalent jcp-np4-ne))
+
+(when-inp ["figs-acc-jcp"]
+  (sv c (jcp-context (get-context))
+      c.nes [10 20 40 80 160]
+      d (acc-parse (+ c.data-dir "acc-jcp.txt")
+                   :map-nstepfac jcp-nenp2nstepfac))
+  (figs-acc-jcp c d))
+
+(when-inp ["figs-vortex-acc-jcp"]
+  (sv c (jcp-context (get-context))
+      c.nes [10 20 40 80 160]
+      d (acc-parse (+ c.data-dir "vortex-jcp.txt")
+                   :map-nstepfac jcp-nenp2nstepfac))
+  (figs-vortex-acc-jcp c d))
+
+(when-inp ["figs-mixing-jcp"]
+  (sv c (jcp-context (get-context))
+      d (acc-parse (+ c.data-dir "mixing-jcp.txt")
+                   :map-nstepfac jcp-nenp2nstepfac))
+  (figs-mixing-jcp c d))
+
+(when-inp ["fig-filament-jcp"]
+  (sv c (jcp-context (get-context))
+      d (acc-parse (+ c.data-dir "mixing-jcp.txt")
+                   :map-nstepfac jcp-nenp2nstepfac))
+  (fig-filament-jcp c d))
diff --git a/methods/islet/figures/figs-methods.hy b/methods/islet/figures/figs-methods.hy
index b6753d5..167bf0b 100644
--- a/methods/islet/figures/figs-methods.hy
+++ b/methods/islet/figures/figs-methods.hy
@@ -440,7 +440,8 @@
 (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]
+(defn plot1-slmmir-vs-heuristic [c d nps prop-preserve ic norm pum-thr lebesgue
+                                 &optional [jcp False]]
   (sv npa npy.array fs 11
       plot (if lebesgue pl.semilogy pl.loglog))
   (for [np nps]
@@ -466,7 +467,10 @@
                    :fontsize fs)))
   (my-grid)
   (pl.title (.format (+ "Test problem vs.~heuristic:\n"
-                        "nondivergent flow, 1.5$^\circ$, {}, long steps\n"
+                        "nondivergent flow, "
+                        (if jcp "$n_e = 20$, " "1.5$^\circ$, ")
+                        "{}, "
+                        (if jcp "120 steps\n" "long steps\n")
                         "$p$-refinement, {}")
                      (futils.ic-short2long ic)
                      (+ (if prop-preserve "" "no ") "property preservation"))
@@ -498,7 +502,7 @@
       (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)
+                                 lebesgue :jcp True)
       (sv f (pl.gcf))
       (f.text (nth (, 0.02 0.52) i) 0.05 (nth (, "(a)" "(b)") i)
               :fontsize fs))))
@@ -644,7 +648,7 @@
 ;;; drivers
 
 (when-inp ["dev-parse"]
-  (sv fname "data/search-0.txt"
+  (sv fname "data/mar21/search-0.txt"
       blns (parse-search-list fname))
   (print (len blns))
   (for [b blns]
@@ -664,7 +668,7 @@
                      "search-findnodal_given_bestosn-7.txt")
       blns [])
   (for [fname data-fnames]
-    (.extend blns (parse-search-list (+ "data/" fname))))
+    (.extend blns (parse-search-list (+ "data/mar21/" fname))))
   (sv blns (uniquify-search-list blns))
   (write-slmmir-script blns script-fname))
 
@@ -674,7 +678,7 @@
       lebesgue False
       d {})
   (for [fname fnames]
-    (sv d (parse-slmmir-output (+ "data/" fname) :d d :lebesgue lebesgue)))
+    (sv d (parse-slmmir-output (+ "data/mar21/" fname) :d d :lebesgue lebesgue)))
   (for [(, norm pp ic) (, (, :l2 False "gau") (, :l2 True "cos"))]
     (plot-slmmir-vs-heuristic
       c d (.format (+ "{}slmmir-vs-heuristic-{}-{}-{}" (if lebesgue "-leb" ""))
@@ -688,7 +692,7 @@
       lebesgue False
       d {})
   (for [fname fnames]
-    (sv d (parse-slmmir-output (+ "data/" fname) :d d :lebesgue lebesgue)))
+    (sv d (parse-slmmir-output (+ "data/mar21/" 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"))
@@ -702,13 +706,13 @@
 
 (when-inp ["pum-vs-perturb"]
   (sv fname "pum_perturb_plot-041021.txt"
-      d (parse-pum-vs-perturb (+ "data/" fname))
+      d (parse-pum-vs-perturb (+ "data/mar21/" fname))
       c (futils.get-context))
   (plot-pum-vs-perturb c d (+ c.fig-dir "pum-vs-perturb")))
 
 (when-inp ["meam1-and-pum-vs-dx"]
   (sv c (futils.get-context)
-      data-dir "data/"
+      data-dir "data/mar21/"
       meam1-fname (+ data-dir "run_meam1_sweep-np8.txt")
       method-fnames (zip (, "gll_best" "gll_natural" "uniform_offset_nodal_subset")
                          (lfor fname (, "pum_sweep-np8-gll_best.txt"
diff --git a/methods/islet/figures/figs-vortex.hy b/methods/islet/figures/figs-vortex.hy
new file mode 100644
index 0000000..43f0746
--- /dev/null
+++ b/methods/islet/figures/figs-vortex.hy
@@ -0,0 +1,120 @@
+(require [amb3 [*]])
+(import [amb3 [*]]
+        [figsutils [*]]
+        math colorsys)
+
+(assoc matplotlib.rcParams "savefig.dpi" 300)
+(do (pl-require-type1-fonts))
+
+(defn get-errs [d ne np cycle]
+  (get d "exact" "movingvortices" 1 "pcsl" "none" "none" 0 np "vor" cycle ne :C))
+
+(defn make-vortex-cmap []
+  (sv lsc (fn [colors lo hi n]
+            ((matplotlib.colors.LinearSegmentedColormap.from-list
+              "filament" colors 11)
+             (npy.linspace lo hi n)))
+      h2r colorsys.hsv-to-rgb
+      s 0.2
+      red (lsc [(h2r 0 1 1) (h2r 0 s 1)] 0 1 11)
+      b (/ 2 3)
+      blue (lsc [(h2r b 1 1) (h2r b s 1)] 1 0 11)
+      cmap (matplotlib.colors.ListedColormap
+            (npy.concatenate (, blue red) :axis 0)))
+  cmap)
+
+(defn fig-vortex [c d configs]
+  (print "TODO config labels")
+  (defn draw-circle []
+    (sv th (npy.linspace 0 (* 2 math.pi) 512)
+        r 0.995
+        x (* r (npy.cos th)) y (* r (npy.sin th))
+        g 0.8)
+    (pl.plot x y "-" :lw 0.75 :color [g g g]))
+  (sv ncol (len configs)
+      fs 11
+      cmap (make-vortex-cmap)
+      ecmap (matplotlib.colors.LinearSegmentedColormap.from-list
+             "filament" [[1 1 1] [0 0 0]]))
+  (with [(pl-plot (, (* 2 ncol) 4) (+ c.fig-dir "vortex")
+                  :format "png"
+                  :tight False)]
+    (for [(, icon config) (enumerate configs)]
+      (sv (, ne np nstep cycle recon) config
+          e (get-errs d ne np cycle)
+          arrs (read-slmmir-io-arrays
+                (+ c.data-dir f"vortex-imgs/ne{ne}np{np}nstep{nstep}{recon}.bin")))
+      (sv img (get arrs (+ 1 (* 3 cycle)))
+          diff (/ (npy.abs (- img (get arrs (+ 2 (* 3 cycle)))))
+                  (npy.max (npy.abs (get arrs (+ 2 (* 3 cycle))))))
+          imgs [img diff])
+      (for [i (range 2)]
+        (sv img (get imgs i)
+            mask (!= img 0)
+            ext (, (npy.min (get img mask)) (npy.max img)))
+        (when (zero? i)
+          (sv nmask (npy.logical-not mask)))
+        (sv (get img nmask) npy.nan)
+        (print "corner" (get img (, 0 0)))
+        (print "ext" ext)
+        (if (zero? i)
+          (sv vmin 0.45 vmax 1.55)
+          (sv vmax (get e :li)
+              vmin 0))
+        (sv rect [(/ icon ncol) (- 1 (* i 0.5))
+                  (/ 0.97 ncol) (* 0.97 0.5)]
+            ax (pl.axes rect))
+        (pl.imshow img :cmap (if (zero? i) cmap ecmap)
+                   :vmin vmin :vmax vmax
+                   :origin "upper"
+                   :extent [-1 1 -1 1])
+        (pl.axis "off")
+        (pl.text -1 -1
+                 (.format "({})" (if (zero? i)
+                                   (get "abcd" icon)
+                                   (get "efgh" icon)))
+                 :ha "left" :va "bottom" :fontsize fs)
+        (draw-circle)
+        (unless (zero? i)
+          (pl.text 0 0.5
+                   (.format "$n_e \, {}$, $n_p \, {}$, cycle {}, $n_{{step}} \, {}$\n"
+                            ne np cycle (* cycle nstep))
+                   :ha "center" :fontsize fs)
+          (unless (= recon "constant")
+            (for [k (range 3)]
+              (pl.text -0.4 (+ 0.0 (* 0.2 (- 2 k)))
+                       (.format (+ (get ["$l_1$" "$l_2$" "$l_{{\infty}}$"] k)
+                                   " {:1.2e}")
+                                (get e (get (, :l1 :l2 :li) k)))
+                       :ha "right" :fontsize fs)))
+          (when (= recon "constant")
+            (pl.text 0 0.9 "constant in subcell" :ha "center" :fontsize fs)))
+        (cond [(and (zero? i) (= icon 0))
+               (sv f (pl.gcf)
+                   cax (f.add-axes [(- (get rect 0) 0.04)
+                                    (+ (get rect 1) (* 0.1 (get rect 3)))
+                                    (* 0.1 (get rect 2))
+                                    (* 0.8 (get rect 3))]))
+               (pl.colorbar :cax cax :ticks [0.5 0.75 1 1.25 1.5]
+                            :aspect 10)
+               (cax.yaxis.set-ticks-position "left")]
+              [(= i 1)
+               (sv f (pl.gcf)
+                   cax (f.add-axes [(+ (get rect 0) (* 0.1 (get rect 2)))
+                                    (+ (get rect 1) 0.1)
+                                    (* 0.8 (get rect 2))
+                                    (* 0.05 (get rect 3))])
+                   cb (pl.colorbar :cax cax :orientation "horizontal"
+                                   :ticks [0 vmax]
+                                   :shrink 0.6 :aspect 10))
+               (cb.ax.set-xticklabels ["0" (.format "{:1.2e}" vmax)])])))))
+
+(when-inp ["fig-vortex"]
+  (sv c (jcp-context (get-context))
+      d (acc-parse (+ c.data-dir "vortex-imgs.txt")
+                   :map-nstepfac jcp-nenp2nstepfac)
+      configs [(, 3 13 72 1 "constant")
+               (, 3 13 72 1 "bilin")
+               (, 6 13 144 1 "bilin")
+               (, 12 13 288 2 "bilin")])
+  (fig-vortex c d configs))
diff --git a/methods/islet/figures/figsutils.hy b/methods/islet/figures/figsutils.hy
index f359731..fefcb4f 100644
--- a/methods/islet/figures/figsutils.hy
+++ b/methods/islet/figures/figsutils.hy
@@ -3,8 +3,8 @@
 
 (defn get-context []
   (sv c (Box)
-      c.data-dir "data/"
-      c.fig-dir "figs/"
+      c.data-dir "data/mar21/"
+      c.fig-dir "figs/drafts/"
       c.odes (, "rotate" "divergent" "nondivergent")
       c.cdrs (, (, "none" "none") (, "caas-node" "caas") (, "caas" "caas"))
       c.nstepfacs (, 1 5)
@@ -14,32 +14,88 @@
       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 "m" 13 "m" 16 "g"}
+      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.npmarks {4 "o" 5 "x" 6 "s" 7 "x" 8 "p" 9 "+" 10 "." 11 "^" 12 "." 13 "*" 16 "."})
   c)
 
 (defn flow-short2long [flow]
   (get {"divergent" "divergent flow"
         "nondivergent" "nondivergent flow"
-        "rotate" "solid-body rotation"} flow))
+        "rotate" "solid-body rotation"
+        "movingvortices" "moving vortices"} flow))
 
 (defn ic-short2long [ic]
   (get {"gau" "Gaussian hills"
         "cos" "cosine bells"
         "slo" "slotted cylinders"} ic))
 
-(defn nes->degstrs [nes]
+(defn nes->degstrs [nes &optional [convert-all False]]
   (sv x [] xstr [])
   (for [ne nes]
-    (sv deg (geton {5 6 10 3 20 "1.5" 40 "0.75" 80 "0.375" 160 "0.1875"} ne))
-    (when (none? deg) (continue))
+    (sv deg (geton {2 15 5 6 10 3 20 "1.5" 40 "0.75" 80 "0.375" 160 "0.1875"} ne))
     (.append x ne)
-    (.append xstr (.format "${}^{{\circ}}$" deg)))
+    (if (none? deg)
+      (do (sv deg (/ 30 ne))
+          (.append xstr (.format "${:1.4f}^{{\circ}}$" deg)))
+      (.append xstr (.format "${}^{{\circ}}$" deg))))
   {:ne x :deg xstr})
 
 (defn cdr-name [short]
   (get {"caas" "CAAS-CAAS" "caas-node" "CAAS-point"} short))
 
+;;; parse cmd, L, M, C slmmir output
+
+(defn acc-parse [fname &optional map-nstepfac]
+  (sv txt (.split (readall fname) "\n")
+      bo None d {})
+  (for [ln txt]
+    (sv ln2 (cut ln 0 2))
+    (cond [(and (= "cm" ln2) (in "cmd>" ln))
+           (sv cmd ln c (parse-cmd ln :map-nstepfac map-nstepfac))]
+          [(= ln2 "M ")
+           (sv mp (parse-midpoint-check ln))
+           (unless (none? (geton d #*(+ (cmd->key-base c)
+                                        (, (:ic mp) cyc (:ne c) :M))))
+             ;; handle repeated ic used for src-term OOA testing.
+             (assoc mp :ic (+ (:ic mp) "-src")))
+           (assoc-nested d (+ (cmd->key-base c) (, (:ic mp) cyc (:ne c) :M))
+                         {:l1 (:l1 mp) :l2 (:l2 mp)})]
+          [(= ln2 "C ")
+           (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") (!= (:mono c) "none") (<= cyc 10))
+                    (prf "FAIL {}" cmd))]
+                 [:else
+                  (sv cl (parse-C ln))
+                  (unless (none? (geton d #*(+ (cmd->key-base c)
+                                               (, (:ic cl) cyc (:ne c) :C))))
+                    ;; handle repeated ic used for src-term OOA testing.
+                    (assoc cl :ic (+ (:ic cl) "-src")))
+                  #_(print "found:" (+ (cmd->key-base c) (, (:ic cl) cyc (:ne c) :C)))
+                  (assoc-nested d (+ (cmd->key-base c) (, (:ic cl) cyc (:ne c) :C))
+                                cl)])]
+          [(= ln2 "L ")
+           (cond [(in "L file" ln)
+                  (assoc-nested d (+ (cmd->key-base c)
+                                     (, "cos" cyc (:ne c) :L :mixing-file))
+                                (last (.split ln)))
+                  ;; If nothing else has marked bo as done, do it now.
+                  (unless (none? bo)
+                    (assoc bo :done True))]
+                 [(in "phimin" ln)
+                  (sv (, - ic - l1 - l2 - li - phimin - phimax)
+                      (sscanf ln "s,s,s,f,s,f,s,f,s,f,s,f"))
+                  (assoc-nested d (+ (cmd->key-base c) (, ic cyc (:ne c) :Lerr))
+                                {:l1 l1 :l2 l2 :li li :phimin phimin :phimax phimax})]
+                 [:else
+                  (sv bo (parse-bakeoff-diag bo ln (:timeint c)))])
+           (when (and (not (none? bo)) (:done bo))
+             (for [(, k v) (.items bo)]
+               (assoc-nested d (+ (cmd->key-base c) (, "cos" cyc (:ne c) :L k)) v))
+             (sv bo None))]))
+  d)
+
 ;;; slmmir I/O
 
 (defn read-slmmir-io-arrays [fname &optional beg end stride]
@@ -91,6 +147,7 @@
 
 ;;; parse slmmir text output
 
+;; (map-nstepfac nstepfac-initial ne np)
 (defn parse-cmd [cmd &optional map-nstepfac]
   (sv toks (.split cmd))
   (defn int-or-none [x]
@@ -105,7 +162,8 @@
   (for [e (.items keys)]
     (assoc d (keyword (first e)) ((second e) (get-key-val (first e)))))
   (sv nstepfac (/ (:nsteps d) (:ne d) 6))
-  (unless (none? map-nstepfac) (sv nstepfac (map-nstepfac nstepfac)))
+  (unless (none? map-nstepfac)
+    (sv nstepfac (map-nstepfac nstepfac (:ne d) (:np d))))
   (assoc d :nstepfac (int nstepfac))
   (when (= (:timeint d) "exact") (assoc d :prefine 0))
   d)
@@ -224,3 +282,32 @@
   (if (in " offst " ln)
       (parse-search-offset-nodal-subset-line ln)
       (parse-search-nodal-subset-line ln)))
+
+;;; JCP manuscript revision utilities
+
+(defn jcp-context [c]
+  (sv c.data-dir "data/feb24/"
+      c.nps (, 4 6 8 11 13))
+  c)
+
+;; Compute ne for np given tne for np=4 such that (ne,np) has the maximal #DOF
+;; <= $DOF for (tne,4).
+(defn jcp-tne2ne [tne np]
+  (// (* tne 3) (dec np)))
+
+;; The floating point actual value of resolution converted to np=4.
+(defn jcp-effective-tne [tne np]
+  (/ (* (jcp-tne2ne tne np) (dec np)) 3))
+
+(defn jcp-np4-ne [ne np]
+  (/ (* ne (dec np)) 3))
+
+;; Deduce nstepfac 1 or 5 from parse-cmd's internal values.
+(defn jcp-nenp2nstepfac [nstepfac ne np]
+  (sv nsteps (* nstepfac ne 6)
+      tne (/ (* ne (dec np)) 3)
+      tnstepfac (/ nsteps tne 6))
+  (cond [(< tnstepfac 3) 1]
+        [(< tnstepfac 8) 5]
+        [:else (raisefmt "jcp-nenp2nstepfac failed: {} {} {} => {}"
+                         nstepfac ne np tnstepfac)]))
diff --git a/methods/islet/figures/run-accuracy-jcp.sh b/methods/islet/figures/run-accuracy-jcp.sh
new file mode 100644
index 0000000..682c201
--- /dev/null
+++ b/methods/islet/figures/run-accuracy-jcp.sh
@@ -0,0 +1,28 @@
+cat $0
+
+exe=slmmir
+
+ctr=0
+function run {
+    ctr=$(expr $ctr + 1)
+    cmd="OMP_NUM_THREADS=64 KMP_AFFINITY=balanced ../$exe -method pcsl -ic gaussianhills -ic cosinebells -ic correlatedcosinebells -ic slottedcylinders -we 0 -rit -T 12 -d2c -ode $ode -ne $ne -np $np -nsteps $nstep -prefine 0 -mono none -lim none -timeint exact -rotate-grid"
+    grepcmd='grep "^C \|^L \|^M "'
+    echo "cmd> $ctr $cmd"
+    eval "$cmd | $grepcmd"
+}
+
+for tne in 5 10 20 40 80 160; do
+    for nstepfac in 1 5; do
+        for ode in nondivergent; do
+            tgtres=$(($tne * 3))
+            for np in $(seq 4 13) 16; do
+                ne=$(($tgtres / ($np - 1)))
+                nstep=$(($nstepfac * 2 * ($np - 1) * $ne))
+                if [[ $ne == 1 ]]; then
+                    continue
+                fi
+                run
+            done
+        done
+    done
+done
diff --git a/methods/islet/figures/run-mixing-jcp.sh b/methods/islet/figures/run-mixing-jcp.sh
new file mode 100644
index 0000000..8e66b40
--- /dev/null
+++ b/methods/islet/figures/run-mixing-jcp.sh
@@ -0,0 +1,33 @@
+cat $0
+
+exe=slmmir
+
+ctr=0
+function run {
+    ctr=$(expr $ctr + 1)
+    cmd="OMP_NUM_THREADS=1 KMP_AFFINITY=balanced ../$exe -method pcsl -ic gaussianhills -ic cosinebells -ic correlatedcosinebells -ic slottedcylinders -we 0 -rit -T 12 -d2c -lauritzen -lauritzen-io -ode $ode -ne $ne -np $np -nsteps $nstep -prefine $prefine -mono $cdrglb -lim $cdrlcl -timeint $timeint -rotate-grid -o mixing-jcp/${ode}-$timeint-nsteps${nstep}-prefine${prefine}-$cdrglb-$cdrlcl-ne$ne-np$np"
+    grepcmd='grep "^C \|^L \|^M "'
+    echo "cmd> $ctr $cmd"
+    eval "$cmd | $grepcmd"
+}
+
+cdrglbs=(none caas-node);
+cdrlcls=(none caas);
+for tne in 20 40; do
+    tgtres=$(($tne * 3))
+    for nstepfac in 1 5; do
+        for ode in nondivergent; do
+            icdr=0
+            cdrglb=${cdrglbs[$icdr]}
+            cdrlcl=${cdrlcls[$icdr]}
+            timeint=exact
+            prefine=0
+            for np in 4 6 8 11 13; do
+                ne=$(($tgtres / ($np - 1)))
+                nstep=$(($nstepfac * 2 * ($np - 1) * $ne))
+                echo ">>> $ne $np $nstep"
+                run
+            done
+        done
+    done
+done
diff --git a/methods/islet/figures/run-vortex-imgs.sh b/methods/islet/figures/run-vortex-imgs.sh
new file mode 100644
index 0000000..8206c6e
--- /dev/null
+++ b/methods/islet/figures/run-vortex-imgs.sh
@@ -0,0 +1,36 @@
+cat $0
+
+exe=slmmir
+
+ctr=0
+function run {
+    ctr=$(expr $ctr + 1)
+    cmd="OMP_NUM_THREADS=1 KMP_AFFINITY=balanced ../$exe -method pcsl -ic vortextracer -rit -T 24 -d2c -ode movingvortices -ne $ne -np $np -nsteps $nstep -prefine $prefine -mono $cdrglb -lim $cdrlcl -timeint $timeint -rotate-grid -we $nstep -io-type internal -io-recon $iorecon -res 1024 -o vortex-imgs/ne${ne}np${np}nstep${nstep}${iorecon}"
+    grepcmd='grep "^C \|^L \|^M "'
+    echo "cmd> $ctr $cmd"
+    eval "$cmd | $grepcmd"
+}
+
+cdrglb=none
+cdrlcl=none
+timeint=exact
+prefine=0
+iorecon=bilin
+
+ne=3
+np=13
+nstep=$((12 * 6))
+run
+iorecon=constant
+run
+iorecon=bilin
+
+ne=6
+np=13
+nstep=$((24 * 6))
+run
+
+ne=12
+np=13
+nstep=$((48 * 6))
+run
diff --git a/methods/islet/figures/run-vortex-jcp.sh b/methods/islet/figures/run-vortex-jcp.sh
new file mode 100644
index 0000000..4a57f4a
--- /dev/null
+++ b/methods/islet/figures/run-vortex-jcp.sh
@@ -0,0 +1,31 @@
+cat $0
+
+exe=slmmir
+
+ctr=0
+function run {
+    ctr=$(expr $ctr + 1)
+    cmd="OMP_NUM_THREADS=64 KMP_AFFINITY=balanced ../$exe -method pcsl -ic vortextracer -we 0 -rit -T 24 -d2c -ode movingvortices -ne $ne -np $np -nsteps $nstep -prefine $prefine -mono $cdrglb -lim $cdrlcl -timeint $timeint -rotate-grid"
+    grepcmd='grep "^C \|^L \|^M "'
+    echo "cmd> $ctr $cmd"
+    eval "$cmd | $grepcmd"
+}
+
+cdrglb=none
+cdrlcl=none
+timeint=exact
+prefine=0
+for nstepfac in 1 5; do
+    for tne in 5 10 20 40 80 160; do
+        tgtres=$(($tne * 3))
+        for np in $(seq 4 13) 16; do
+            ne=$(($tgtres / ($np - 1)))
+            nstep=$(($nstepfac * 2 * ($np - 1) * $ne))
+            if [[ $ne == 1 ]]; then
+                continue
+            fi
+            echo ">>> $ne $np $nstep"
+            run
+        done
+    done
+done
diff --git a/methods/slmm/slmm_gallery.cpp b/methods/slmm/slmm_gallery.cpp
index 8b65d00..a41eef2 100644
--- a/methods/slmm/slmm_gallery.cpp
+++ b/methods/slmm/slmm_gallery.cpp
@@ -138,12 +138,7 @@ void InitialCondition::init (const Shape shape, const Size n, const Real* const
     }      
   } break;
   case VortexTracer: {
-    for (Size i = 0; i < n; ++i) {
-      const Real lambda_p = std::atan2(-std::cos(lon[i]), std::tan(lat[i]));
-      const Real rr = 3*std::sqrt((1 - slmm::square(std::cos(lat[i])) *
-                                   slmm::square(std::sin(lon[i]))));
-      u[i] = 0.5*(1 - std::tanh(0.2*rr*std::sin(lambda_p)));
-    }
+    MovingVortices::calc_tracer(0, n, lat, lon, u);
   } break;
   case CorrelatedCosineBells: {
     const Real a = -0.8, b = 0.9;
@@ -396,10 +391,45 @@ ::eval (const Real t, const Real* const d, Real* const f) const {
 //     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.
+// The code uses the formulas in
+//     Bosler, P. A., Kent, J., Krasny, R., & Jablonowski, C. (2017). A
+//     Lagrangian particle method with remeshing for tracer transport on the
+//     sphere. Journal of Computational Physics, 340, 639-654.
+
+const Real MovingVortices::rho0 = 3;
+const Real MovingVortices::gamma = 5;
+
+Real MovingVortices::get_Omega () {
+  return 2*M_PI/slmm::day2sec(12);
+}
+
+Real MovingVortices::calc_rho (const Real theta, const Real lambda) {
+  return rho0*std::sqrt(1 - slmm::square(std::cos(theta)*std::sin(lambda)));
+}
+
+Real MovingVortices::calc_omega (const Real Omega, const Real rho) {
+  return (rho == 0 ?
+          0 :
+          (Omega*slmm::consts::earth_radius_m*
+           1.5*std::sqrt(3)*std::tanh(rho) /
+           (rho*slmm::square(std::cosh(rho)))));
+}
+
+void MovingVortices
+::calc_tracer (const Real time, const Size n, const Real* const lat,
+               const Real* const lon, Real* const u) {
+  const auto Omega = get_Omega(), R = slmm::consts::earth_radius_m;
+  for (Size i = 0; i < n; ++i) {
+    const Real
+      lon_d = lon[i] - Omega*time,
+      lambda_p = std::atan2(-std::cos(lon_d), std::tan(lat[i])),
+      rho = calc_rho(lat[i], lon_d),
+      omega = calc_omega(Omega, rho);
+    u[i] = 1 - std::tanh((rho/MovingVortices::gamma)*
+                         std::sin(lambda_p - (omega/R)*time));
+  }
+}
+
 bool MovingVortices
 ::eval (const Real t, const Real* const d, Real* const f) const {
   Real theta, lambda;
@@ -410,25 +440,20 @@ ::eval (const Real t, const Real* const d, Real* const f) const {
     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),
+    R = slmm::consts::earth_radius_m,
+    Omega = get_Omega(),
     // 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))));
+    lambda_p = lambda - fac*Omega*t,
+    costh = std::cos(theta),
+    rho = calc_rho(theta, lambda_p),
+    omega = calc_omega(Omega, rho);
   // v
-  f[0] = omg*std::cos(lambda_p)*cost;
+  f[0] = omega*std::cos(lambda_p);
   // u
-  f[1] = omg*std::sin(lambda_p)*std::sin(theta)*cost + u0*costh*fac;
+  f[1] = omega*std::sin(lambda_p)*std::sin(theta) + R*Omega*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 {
diff --git a/methods/slmm/slmm_gallery.hpp b/methods/slmm/slmm_gallery.hpp
index 2991d16..0ca5ff1 100644
--- a/methods/slmm/slmm_gallery.hpp
+++ b/methods/slmm/slmm_gallery.hpp
@@ -60,6 +60,15 @@ struct DivergentWindField : public OdeFn {
 
 struct MovingVortices : public OdeFn {
   bool eval(const Real t, const Real* const d, Real* const f) const override;
+  static const Real rho0, gamma;
+  // Omega is the solid-body rotation rate.
+  static Real get_Omega();
+  static Real calc_rho(const Real theta, const Real lambda);
+  // omega is the vortex rotation.
+  static Real calc_omega(const Real Omega, const Real rho);
+  // Time is in seconds.
+  static void calc_tracer(const Real time, const Size n, const Real* const lat,
+                          const Real* const lon, Real* const u);
 };
 
 struct NonDivergentWindFieldHack : public OdeFn {
diff --git a/methods/slmm/slmm_test.cpp b/methods/slmm/slmm_test.cpp
index 6c79afa..0f283c1 100644
--- a/methods/slmm/slmm_test.cpp
+++ b/methods/slmm/slmm_test.cpp
@@ -17,8 +17,7 @@ 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;
diff --git a/methods/slmm/slmm_vis.cpp b/methods/slmm/slmm_vis.cpp
index d5fd59d..98eaee5 100644
--- a/methods/slmm/slmm_vis.cpp
+++ b/methods/slmm/slmm_vis.cpp
@@ -10,6 +10,17 @@
 namespace slmm {
 namespace vis {
 
+Reconstruction::Enum Reconstruction::convert (const std::string& s) {
+  if (s == "constant") return Enum::constant;
+  return Enum::bilin;
+}
+
+std::string Reconstruction::convert (const Enum e) {
+  if (e == Enum::constant) return "constant";
+  return "bilin";
+}
+
+
 struct SourceType {
   enum Enum { gll, physgrid };
 };
@@ -57,8 +68,12 @@ struct SparseMat {
 // y = A*x
 static void csr_mvp (const SparseMat& A, const Real* x, Real* y) {
   const Int nr = (Int) A.r.size() - 1;
-  Int nc = 0;
+# pragma omp parallel for
   for (Int r = 0; r < nr; ++r) {
+    if (A.r[r+1] == A.r[r]) {
+      y[r] = 0;
+      continue;
+    }
     Real accum = 0;
     Real min_x = 1e10, max_x = -1e10;
     for (Int j = A.r[r]; j < A.r[r+1]; ++j) {
@@ -72,7 +87,6 @@ static void csr_mvp (const SparseMat& A, const Real* x, Real* y) {
     // diagnostic. So limit it.
     accum = std::min(max_x, std::max(accum, min_x));
     y[r] = accum;
-    for (Int j = A.r[r]; j < A.r[r+1]; ++j) nc = std::max(nc, A.c[j]);
   }
 }
 
@@ -86,10 +100,39 @@ static void cclinspace (const Int n, const Real lo, const Real hi,
 }
 
 static void
-init_interp_map (const AVec3s& p, const AIdxs& c2n,
-                 const Octree& ot, const std::vector<Real>& lats,
-                 const std::vector<Real>& lons, const SourceType::Enum src_type,
-                 SparseMat& op) {
+fill_row (const AVec3s& p, const AIdxs& c2n, const Int ci, const Int r,
+          const Real* s, const SourceType::Enum src_type,
+          const Reconstruction::Enum recon, SparseMat& op) {
+  const auto& cell = slice(c2n, ci);
+  Real* ls = &op.s[op.r[r]];
+  if (src_type == SourceType::gll) {
+    // Compute reference coordinates.
+    Real a, b;
+    siqk::sqr::Info info;
+    siqk::sqr::calc_sphere_to_ref(p, cell, s, a, b, &info);
+    // Fill matrix row.
+    for (Int i = 0; i < 4; ++i) op.c[op.r[r]+i] = cell[i];
+    if (recon == Reconstruction::Enum::bilin) {
+      static const Real f = 0.25;
+      ls[0] = f*(1-a)*(1-b); ls[1] = f*(1+a)*(1-b);
+      ls[2] = f*(1+a)*(1+b); ls[3] = f*(1-a)*(1+b);
+    } else {
+      ls[0] = ls[1] = ls[2] = ls[3] = 0.25;
+    }
+  } else {
+    op.c[op.r[r]] = ci;
+    *ls = 1;
+  }
+}
+
+// Map for const or bilin interp for lat-lon tensor-grid output. c2n should be
+// the I/O cell->node map, i.e., subcells of the full element.
+static void init_latlon_ioc2n_map (
+  const AVec3s& p, const AIdxs& c2n, const Octree& ot,
+  const std::vector<Real>& lats, const std::vector<Real>& lons,
+  const SourceType::Enum src_type, const Reconstruction::Enum recon,
+  SparseMat& op)
+{
   // Compute edge normals for octree search.
   AVec3s nml;
   AIdxs c2nml;
@@ -112,24 +155,55 @@ init_interp_map (const AVec3s& p, const AIdxs& c2n,
     ot.apply(sf.bb, sf);
     const Int ci = sf.ci_hit;
     assert(ci >= 0 && ci < nslices(c2n));
-    Real* ls = &op.s[op.r[r]];
-    const auto& cell = slice(c2n, ci);
-    if (src_type == SourceType::gll) {
-      // Compute reference coordinates.
-      Real a, b;
-      siqk::sqr::Info info;
-      siqk::sqr::calc_sphere_to_ref(p, cell, s, a, b, &info);
-      // Fill matrix row.
-      for (Int i = 0; i < 4; ++i) op.c[op.r[r]+i] = cell[i];
-      {
-        static const Real f = 0.25;
-        ls[0] = f*(1-a)*(1-b); ls[1] = f*(1+a)*(1-b);
-        ls[2] = f*(1+a)*(1+b); ls[3] = f*(1-a)*(1+b);
-      }
+    fill_row(p, c2n, ci, r, s, src_type, recon, op);
+  }
+}
+
+// Map for const or bilin interp for orthographic output. c2n should be the I/O
+// cell->node map, i.e., subcells of the full element.
+static void init_orthographic_ioc2n_map (
+  const AVec3s& p, const AIdxs& c2n, const Octree& ot,
+  const std::vector<Real>& xs, const std::vector<Real>& ys,
+  const Real xhat[3], const Real yhat[3], const Real zhat[3],
+  const Reconstruction::Enum recon, SparseMat& op)
+{
+  using geo = siqk::SphereGeometry;
+  // Compute edge normals for octree search.
+  AVec3s nml;
+  AIdxs c2nml;
+  fill_normals(p, c2n, nml, c2nml);
+  // Make a sparse CSR matrix implementing bilinear interpolation. We know most
+  // of the structure of this matrix up front.
+  const Int nx = (Int) xs.size(), nr = (Int) ys.size() * nx;
+  const Int nnz_per_row = 4;
+  op.r.resize(nr+1);
+  op.r[0] = 0;
+  for (Int r = 0; r < nr; ++r) {
+    const Int ix = r % nx, iy = r / nx;
+    const Real xy_mag2 = square(xs[ix]) + square(ys[iy]);
+    op.r[r+1] = op.r[r] + (xy_mag2 < 1 ? nnz_per_row : 0);
+  }
+  op.c.resize(op.r[nr], 0);
+  op.s.resize(op.r[nr], 0);
+# pragma omp parallel for
+  for (Int r = 0; r < nr; ++r) {
+    // Find the cell containing this point.
+    const Int ix = r % nx, iy = r / nx;
+    const Real xy_mag2 = square(xs[ix]) + square(ys[iy]);
+    Real s[3] = {0};
+    if (xy_mag2 < 1) {
+      const Real z = std::sqrt(1 - xy_mag2);
+      geo::axpy(xs[ix], xhat, s);
+      geo::axpy(ys[iy], yhat, s);
+      geo::axpy(z, zhat, s);
     } else {
-      op.c[op.r[r]] = ci;
-      *ls = 1;
+      continue;
     }
+    SearchFunctor sf(p, c2n, nml, c2nml, s);
+    ot.apply(sf.bb, sf);
+    const Int ci = sf.ci_hit;
+    assert(ci >= 0 && ci < nslices(c2n));
+    fill_row(p, c2n, ci, r, s, SourceType::gll, recon, op);
   }
 }
 
@@ -139,11 +213,13 @@ struct MapToLatLon::Impl {
   SparseMat op;
 
   Impl (const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
-        Int nlat, Int nlon, const SourceType::Enum src_type) {
+        Int nlat, Int nlon, const SourceType::Enum src_type,
+        const Reconstruction::Enum recon) {
     cclinspace(nlat, -M_PI/2, M_PI/2, lats);
     cclinspace(nlon, -M_PI, M_PI, lons);
     ot.init(a2ConstVec3s(cgll_p), a2ConstIdxs(cgll_io_c2n));
-    init_interp_map(cgll_p, cgll_io_c2n, ot, lats, lons, src_type, op);
+    init_latlon_ioc2n_map(cgll_p, cgll_io_c2n, ot, lats, lons, src_type, recon,
+                          op);
   }
 
   // Input is CGLL data, e.g., from D2Cer::d2c. Output is lon-lat data in a
@@ -155,17 +231,21 @@ struct MapToLatLon::Impl {
 
 const std::vector<Real>& MapToLatLon::get_lons () const { return p->lons; }
 const std::vector<Real>& MapToLatLon::get_lats () const { return p->lats; }
+size_t MapToLatLon::get_x_size () const { return get_lons().size(); }
+size_t MapToLatLon::get_y_size () const { return get_lats().size(); }
 
 BilinGLLToLatLon
 ::BilinGLLToLatLon (const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
-                    Int nlat, Int nlon) {
-  p = std::make_shared<Impl>(cgll_p, cgll_io_c2n, nlat, nlon, SourceType::gll);
+                    Int nlat, Int nlon, Reconstruction::Enum recon) {
+  p = std::make_shared<Impl>(cgll_p, cgll_io_c2n, nlat, nlon, SourceType::gll,
+                             recon);
 }
 
 PhysgridToLatLon
 ::PhysgridToLatLon (const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
                     Int nlat, Int nlon) {
-  p = std::make_shared<Impl>(cgll_p, cgll_io_c2n, nlat, nlon, SourceType::physgrid);
+  p = std::make_shared<Impl>(cgll_p, cgll_io_c2n, nlat, nlon, SourceType::physgrid,
+                             Reconstruction::Enum::constant);
 }
 
 // Input is CGLL data, e.g., from D2Cer::d2c. Output is lon-lat data in a
@@ -178,29 +258,71 @@ void PhysgridToLatLon::remap (const Real* field_pg, Real* field_ll) const {
   p->remap(field_pg, field_ll);
 }
 
+struct BilinGLLToOrthographic::Impl {
+  Real zhat[3];
+  std::vector<Real> xs, ys;
+  Octree ot;
+  SparseMat op;
+
+  Impl (const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
+        const Real xhat[3], const Real yhat[3],
+        Int nx, Int ny, Reconstruction::Enum recon)
+  {
+    using geo = siqk::SphereGeometry;
+    assert(std::abs(geo::norm2(xhat) - 1) < 1e-10);
+    assert(std::abs(geo::norm2(yhat) - 1) < 1e-10);
+    geo::cross(xhat, yhat, zhat);
+    assert(std::abs(geo::norm2(zhat) - 1) < 1e-10);
+    const auto r = 1;
+    cclinspace(nx, -r, r, xs);
+    cclinspace(ny, -r, r, ys);
+    ot.init(a2ConstVec3s(cgll_p), a2ConstIdxs(cgll_io_c2n));
+    init_orthographic_ioc2n_map(cgll_p, cgll_io_c2n, ot, xs, ys,
+                                xhat, yhat, zhat, recon, op);
+  }
+
+  void remap (const Real* field_cgll, Real* field_xy) const {
+    csr_mvp(op, field_cgll, field_xy);
+  }
+};
+
+BilinGLLToOrthographic
+::BilinGLLToOrthographic (const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
+                          const Real xhat[3], const Real yhat[3],
+                          Int nx, Int ny, Reconstruction::Enum recon)
+{
+  p = std::make_shared<Impl>(cgll_p, cgll_io_c2n, xhat, yhat, nx, ny, recon);
+}
+
+void BilinGLLToOrthographic::remap (const Real* field_cgll, Real* field_xy) const {
+  p->remap(field_cgll, field_xy);
+}
+
+size_t BilinGLLToOrthographic::get_x_size () const { return p->xs.size(); }
+size_t BilinGLLToOrthographic::get_y_size () const { return p->ys.size(); }
+
 struct VisWriter::Impl {
-  MapToLatLon::Ptr op;
+  MapToArray::Ptr op;
   io::InternalWriter::Ptr iw;
   std::vector<Real> wrk;
-  Impl (const MapToLatLon::Ptr& op_, const std::string& filename)
+  Impl (const MapToArray::Ptr& op_, const std::string& filename)
     : op(op_), iw(std::make_shared<io::InternalWriter>(filename)),
-      wrk(op->get_lats().size() * op->get_lons().size())
+      wrk(op->get_x_size() * op->get_y_size())
   {}
 };
 
-VisWriter
-::VisWriter (const MapToLatLon::Ptr& op, const std::string& filename) {
+VisWriter::VisWriter (const MapToArray::Ptr& op, const std::string& filename) {
   p = std::make_shared<Impl>(op, filename);
 }
 
 void VisWriter::write (const Real* field_cgll) {
   p->op->remap(field_cgll, p->wrk.data());
-  const Int dims[] = {(Int) p->op->get_lats().size(),
-                      (Int) p->op->get_lons().size()};
+  const Int dims[] = {(Int) p->op->get_y_size(),
+                      (Int) p->op->get_x_size()};
   p->iw->write_array(2, dims, p->wrk.data());
 }
 
-MapToLatLon::Ptr VisWriter::get_map () { return p->op; }
+MapToArray::Ptr VisWriter::get_map () { return p->op; }
 
 } // namespace vis
 } // namespace slmm
diff --git a/methods/slmm/slmm_vis.hpp b/methods/slmm/slmm_vis.hpp
index a4cb50e..7da22e9 100644
--- a/methods/slmm/slmm_vis.hpp
+++ b/methods/slmm/slmm_vis.hpp
@@ -8,7 +8,22 @@
 namespace slmm {
 namespace vis {
 
-class MapToLatLon {
+struct Reconstruction {
+  enum class Enum { bilin, constant };
+  static Enum convert(const std::string& s);
+  static std::string convert(const Enum e);
+};
+
+struct MapToArray {
+  typedef std::shared_ptr<MapToArray> Ptr;
+
+  virtual void remap(const Real* field_from, Real* field_ll) const = 0;
+
+  virtual size_t get_x_size() const = 0;
+  virtual size_t get_y_size() const = 0;
+};
+
+class MapToLatLon : public MapToArray {
 protected:
   class Impl;
   std::shared_ptr<Impl> p;
@@ -19,7 +34,8 @@ class MapToLatLon {
   const std::vector<Real>& get_lons() const;
   const std::vector<Real>& get_lats() const;
 
-  virtual void remap(const Real* field_from, Real* field_ll) const = 0;
+  size_t get_x_size() const override;
+  size_t get_y_size() const override;
 };
 
 struct BilinGLLToLatLon : public MapToLatLon {
@@ -27,7 +43,8 @@ struct BilinGLLToLatLon : public MapToLatLon {
 
   BilinGLLToLatLon(const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
                    Int nlat,  // cc(linspace(-pi/2, pi/2, nlat+1))
-                   Int nlon); // cc(linspace(-pi  , pi  , nlon+1))
+                   Int nlon,  // cc(linspace(-pi  , pi  , nlon+1))
+                   Reconstruction::Enum recon = Reconstruction::Enum::bilin);
 
   // Input is CGLL data, e.g., from D2Cer::d2c. Output is lon-lat data in a
   // rectangle, with longitude the faster index.
@@ -45,6 +62,27 @@ struct PhysgridToLatLon : public MapToLatLon {
   void remap(const Real* field_pg, Real* field_ll) const override;
 };
 
+struct BilinGLLToOrthographic : public MapToArray {
+  typedef std::shared_ptr<BilinGLLToOrthographic> Ptr;
+
+  BilinGLLToOrthographic(const AVec3s& cgll_p, const AIdxs& cgll_io_c2n,
+                         // Image frame.
+                         const Real xhat[3], const Real yhat[3],
+                         Int nx, Int ny,
+                         Reconstruction::Enum recon = Reconstruction::Enum::bilin);
+
+  // Input is CGLL data, e.g., from D2Cer::d2c. Output is xhat-yhat data in a
+  // rectangle, with x the faster index.
+  void remap(const Real* field_cgll, Real* field_xy) const override;
+
+  size_t get_x_size() const override;
+  size_t get_y_size() const override;
+
+private:
+  struct Impl;
+  std::shared_ptr<Impl> p;
+};
+
 class VisWriter {
   class Impl;
   std::shared_ptr<Impl> p;
@@ -52,10 +90,10 @@ class VisWriter {
 public:
   typedef std::shared_ptr<VisWriter> Ptr;
 
-  VisWriter(const MapToLatLon::Ptr& bilin, const std::string& filename);
+  VisWriter(const MapToArray::Ptr& bilin, const std::string& filename);
   void write(const Real* field_cgll);
 
-  MapToLatLon::Ptr get_map();
+  MapToArray::Ptr get_map();
 };
 
 } // namespace vis
diff --git a/methods/slmm/slmmir.cpp b/methods/slmm/slmmir.cpp
index 9030a09..9523215 100644
--- a/methods/slmm/slmmir.cpp
+++ b/methods/slmm/slmmir.cpp
@@ -55,6 +55,9 @@
      according to the DSS.
    -io-start-cycle: Wait to start output, except ICs, until this cycle. Cycle
     count starts at 1.
+   -io-recon: For the internal I/O type, choose the reconstruction from {bilin,
+    const}. 'bilin' and 'const' act on subcells. 'bilin' is bilinear interp;
+    'const' is a simple average of nodal values. Default: bilin.
    -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
@@ -300,12 +303,13 @@ struct IntegrateOptions {
   bool track_footprint;
   // Netcdf or custom.
   IOType::Enum io_type;
+  vis::Reconstruction::Enum io_recon;
   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;
+  bool rhot0, vortex_problem;
 };
 
 struct Input {
@@ -325,7 +329,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;
+  bool rotate_grid, ode_tight;
 
   // Super-duper debug/analysis options.
   struct Debug {
@@ -963,6 +967,13 @@ class Observer {
                (f->ms.back().max - f->ms.front().max) / phi0);
       }
   }
+
+  // Hack! Don't use this for anything other than the vortex problem. That
+  // problem is run w/o returning to its IC, so we need to hack our
+  // infrastructure to make the "IC" the solution at the desired time.
+  Real* get_ic_field_for_vortex_problem_to_overwrite (const int idx) {
+    return field_[idx]->ic.data();
+  }
 };
 
 static void write_binary (const Mesh& m, const AVec3s& departure_p,
@@ -974,6 +985,24 @@ static void write_binary (const Mesh& m, const AVec3s& departure_p,
                m.dglln2cglln.data(), "binary.dat");
 }
 
+static void vortex_problem_hook (
+  const Mesh& m, const Observer::Ptr& so, const D2Cer::Ptr& d2cer,
+  const IntegrateOptions& opts, const vis::VisWriter::Ptr& iw,
+  const int field_idx, const Real time)
+{
+  const auto dnn = len(m.dglln2cglln), cnn = nslices(m.cgll_p);
+  std::vector<Real> wrk(cnn);
+  for (Int i = 0; i < cnn; ++i) {
+    const auto n = slice(m.cgll_p, i);
+    Real lat, lon;
+    xyz2ll(n[0], n[1], n[2], lat, lon);
+    gallery::MovingVortices::calc_tracer(time, 1, &lat, &lon, &wrk[i]);
+  }
+  Real* const soln = so->get_ic_field_for_vortex_problem_to_overwrite(field_idx);
+  d2cer->c2d(wrk.data(), soln);
+  if (iw) iw->write(wrk.data());
+}
+
 static pg::PhysgridData::Ptr
 init_physgrid (const PRefineData& pr, const Int nphys,
                const Real* rho, const Real* q,
@@ -1050,9 +1079,17 @@ static void integrate (
       ncw->end_definition();
     } else if (opts.io_type == IOType::internal) {
       const auto res = opts.output_resolution;
-      const auto bilin = std::make_shared<vis::BilinGLLToLatLon>(
-        m.cgll_p, m.cgll_io_c2n, 2*res+1, 4*res+1);
-      iw = std::make_shared<vis::VisWriter>(bilin, out_fn + ".bin");
+      vis::MapToArray::Ptr map;
+      if (opts.vortex_problem) {
+        Real xhat[] = {-1,0,0}, yhat[] = {0,0,1};
+        siqk::SphereGeometry::normalize(xhat);
+        map = std::make_shared<vis::BilinGLLToOrthographic>(
+          m.cgll_p, m.cgll_io_c2n, xhat, yhat, 2*res+1, 2*res+1, opts.io_recon);
+      } else {
+        map = std::make_shared<vis::BilinGLLToLatLon>(
+          m.cgll_p, m.cgll_io_c2n, 2*res+1, 4*res+1, opts.io_recon);
+      }
+      iw = std::make_shared<vis::VisWriter>(map, out_fn + ".bin");
     }
   }
 
@@ -1105,7 +1142,14 @@ static void integrate (
       if (iic == 0)
         copy(error_data.data(), tracer_p[0]->data(), len);
       if (ncw) ncw->write_field(tracer_name(ics[iic], iic), data);
-      if (iw) iw->write(data);
+      if (iw) {
+        iw->write(data);
+        if (opts.vortex_problem &&
+            ics[iic] == gallery::InitialCondition::VortexTracer) {
+          // Write this field again to match vortex_problem_hook writes.
+          iw->write(data);
+        }
+      }
       if (so) {
         so->add_field(tracer_name(ics[iic], iic),
                       tracer_p[0]->data() + iic*len);
@@ -1324,12 +1368,21 @@ static void integrate (
       }
     }
 
+    const bool report = (((step + 1) % nsteps_per_12days == 0) ||
+                         step == last_step);
+
+    if (opts.vortex_problem && (do_io || report)) {
+      for (Int iic = 0; iic < nics; ++iic) {
+        if (ics[iic] != gallery::InitialCondition::VortexTracer) continue;
+        vortex_problem_hook(m, so, d2cer, opts, iw, iic+1, tf);
+        break;
+      }
+    }
+
     // Observer.
     if (so) {
       if (step % nsteps_per_12days == 0)
         std::cout << "\nC cycle " << cycle << "\n";
-      const bool report = (((step + 1) % nsteps_per_12days == 0) ||
-                           step == last_step);
       so->set_time(tf);
       so->add_obs(0, m, rd->dgbfi_gll(), rd->dgbfi(), nullptr,
                   density_p[1]->data(),
@@ -1506,10 +1559,25 @@ static void run (const Input& in) {
     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;
+      if (in.integrate_options.vortex_problem) {
+        // Follow Nair & Jablonowski (see comment above MovingVortices::eval for
+        // full citation) and orient the grid so that the vortices move over 4
+        // of the 8 cubed-sphere corners, thus probing the effects of the
+        // corners. The grid corners are sometimes considered to be troublesome.
+        gr.axis[0] = 1; gr.axis[1] = 0; gr.axis[2] = 0;
+        // Fixed but random number a little under 1. See below.
+        gr.angle = 0.97654321*M_PI/4;
+      } else {
+        // Rotate by a fixed but random amount to keep the center of the
+        // background solid-body rotation off an element corner. An unmoving
+        // solid-body rotation exactly on a collocation point excites a mild
+        // instability, and since the center never moves, the fact that the
+        // instability is mild doesn't help: after enough cycles, it will show
+        // up. This is true for even np=3.
+        gr.axis[0] = 0.11111; gr.axis[1] = -0.051515; gr.axis[2] = 1;
+        gr.angle = 0.142314*M_PI;
+      }
       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);
@@ -1532,6 +1600,7 @@ static void run (const Input& in) {
     switch (in.timeint) {
     case TimeInt::exact:
       mi = MeshIntegratorFactory::create(in.ode, in.xyz_form, p);
+      if (in.ode_tight) mi->config_for_best_possible_accuracy();
       break;
     case TimeInt::line: {
       StudyTimeIntegratorOptions o;
@@ -1645,13 +1714,16 @@ Input::Input (int argc, char** argv)
   iopts.fitext = false;
   iopts.track_footprint = false;
   iopts.io_type = IOType::netcdf;
+  iopts.io_recon = vis::Reconstruction::Enum::bilin;
   iopts.io_start_cycle = 0;
   iopts.io_no_dss = false;
   iopts.output_resolution = 64;
   iopts.pg = 0;
   iopts.rhot0 = false;
+  iopts.vortex_problem = false;
   p_refine_experiment = 0;
   rotate_grid = false;
+  ode_tight = false;
   bool tq_specified = false, method_specified = false;
   for (int i = 1; i < argc; ++i) {
     const std::string& token = argv[i];
@@ -1706,6 +1778,8 @@ Input::Input (int argc, char** argv)
       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-recon"))
+      iopts.io_recon = vis::Reconstruction::convert(argv[++i]);
     else if (eq(token, "-io-nodss"))
       iopts.io_no_dss = true;
     else if (eq(token, "-res"))
@@ -1732,6 +1806,8 @@ Input::Input (int argc, char** argv)
       iopts.perturb_rho = atof(argv[++i]);
     else if (eq(token, "-rotate-grid"))
       rotate_grid = true;
+    else if (eq(token, "-ode-tight"))
+      ode_tight = true;
     else if (eq(token, "-rhot0"))
       iopts.rhot0 = true;
     else
@@ -1783,6 +1859,8 @@ Input::Input (int argc, char** argv)
   std::cout << "Warning: Turning off Lauritzen diagnostics b/c not built.\n";
   iopts.lauritzen_diag = false;
 #endif
+  iopts.vortex_problem = (gallery::WindFieldType::from_string(ode) ==
+                          gallery::WindFieldType::MovingVortices);
   print(std::cout);
   SIQK_THROW_IF(iopts.subcell_bounds && ! Method::is_isl(iopts.method),
                 "-subcellbounds and not ISL is not supported.");
@@ -1827,7 +1905,10 @@ void Input::print (std::ostream& os) const {
      << "write every (-we): " << write_every << "\n"
      << "d2c (-d2c): " << integrate_options.d2c << "\n"
      << "record_in_time (-rit): " << integrate_options.record_in_time << "\n"
-     << "p-refine exp (-prefine): " << p_refine_experiment << "\n";
+     << "p-refine exp (-prefine): " << p_refine_experiment << "\n"
+     << "rotate grid (-rotate-grid): " << rotate_grid << "\n"
+     << "vortex problem: " << integrate_options.vortex_problem << "\n"
+     << "ODE tightest tol: " << ode_tight << "\n";
   if (integrate_options.fitext)
     os << "fitext\n";
   for (auto& ic : initial_conditions)
diff --git a/methods/slmm/slmmir_lauritzen_diag.cpp b/methods/slmm/slmmir_lauritzen_diag.cpp
index 3240687..82cabb6 100644
--- a/methods/slmm/slmmir_lauritzen_diag.cpp
+++ b/methods/slmm/slmmir_lauritzen_diag.cpp
@@ -18,14 +18,20 @@ static void write_mixing_plot_data (const Real* cb, const Real* ccb, const Int n
   auto* const ccbd = cgdata.data() + cnn;
   d2cer.d2c(cb, cbd);
   d2cer.d2c(ccb, ccbd);
-  FILE* fid = fopen((out_prefix + "-lauritzen-diag.dat").c_str(), "w");
+  const auto fname = out_prefix + "-lauritzen-diag.dat";
+  FILE* fid = fopen(fname.c_str(), "wb");
+  if ( ! fid) {
+    printf("Could not open %s; exiting.\n", fname.c_str());
+    exit(-1);
+  }
   std::vector<float> cgsp(2*cnn);
   for (int i = 0; i < 2*cnn; ++i) cgsp[i] = cgdata[i];
+  for (int i = 0; i < 2*cnn; ++i) if (std::isnan(cgsp[i])) prc(i);
   auto* const cbdsp = cgsp.data();
   auto* const ccbdsp = cgsp.data() + cnn;
   fwrite(&cnn, sizeof(Int), 1, fid);
-  fwrite(cbd, sizeof(*cbdsp), cnn, fid);
-  fwrite(ccbd, sizeof(*ccbdsp), cnn, fid);
+  fwrite(cbdsp, sizeof(*cbdsp), cnn, fid);
+  fwrite(ccbdsp, sizeof(*ccbdsp), cnn, fid);
   fclose(fid);
 }