diff --git a/code/MOR b/code/MOR
index c69a8383..94eacb92 100644
--- a/code/MOR
+++ b/code/MOR
@@ -12,7 +12,7 @@ c     parameters and common blocks for MOR
      $                  wb(lx1*ly1*lz1*lelm,0:lub),
      $                  uvwb(lx1*ly1*lz1*lelm,ldim,0:lub),
      $                  pb(lx2*ly2*lz2*lelm,0:lpb),
-     $                  tb(lx1*ly1*lz1*lelm,0:ltb)
+     $                  tb(lx1*ly1*lz1*lelm,0:ltb,ldimt)
 
       common /moric/ uic(lx1*ly1*lz1*lelm),
      $               vic(lx1*ly1*lz1*lelm),
@@ -25,7 +25,7 @@ c     parameters and common blocks for MOR
 
       common /morusnap/ us0(lx1*ly1*lz1*lelm,ldim,lsu)
       common /morpsnap/ prs(lx2*ly2*lz2*lelm,lsp)
-      common /mortsnap/ ts0(lx1*ly1*lz1*lelm,lst)
+      common /mortsnap/ ts0(lx1*ly1*lz1*lelm,lst,ldimt)
 
       common /morsnaptmp/ snapt(lx1*ly1*lz1*lelm,ldim,lsu)
 
@@ -166,8 +166,9 @@ c     parameters and common blocks for MOR
       ! Flags
 
       logical ifl2,ifvort,ifstep,ifpart,ifavg0,ifplay,
-     $        ifcdrag,ifrms,ifread,ifreads(3),ifforce,ifpod(0:ldimt1),
-     $        ifbuoy,ifrom(0:ldimt1),ifcintp,iffastc,ifcnuss,ifctke,
+     $        ifcdrag,ifrms,ifread,ifreads(0:ldimt1),ifforce,
+     $        ifpod(0:ldimt1),ifbuoy,ifrom(0:ldimt1),
+     $        ifcintp,iffastc,ifcnuss,ifctke,
      $        iffasth,ifavisc,ifsource,ifei,ifsub0,ifrecon,iftneu,
      $        ifdecpl,ifcore,ifquad,ifhelm,ifsrct,iftflux,ifcomb,ifpb,
      $        ifcp,ifcflow
diff --git a/code/ana.f b/code/ana.f
index acec0790..5e94edc8 100644
--- a/code/ana.f
+++ b/code/ana.f
@@ -420,8 +420,8 @@ subroutine snap_analysis
 
 
          if (ifpod(2)) then
-            call ps2b(ut,ts0(1,j),tb) 
-            call add2(ts0(1,j),tb,n)
+            call ps2b(ut,ts0(1,j,1),tb)
+            call add2(ts0(1,j,1),tb,n)
          endif
 
          nio = nid
@@ -432,7 +432,7 @@ subroutine snap_analysis
          call opcopy(t1,t2,t3,us0(1,1,j),us0(1,2,j),us0(1,ldim,j))
          energy=op_glsc2_wt(t1,t2,t3,t1,t2,t3,bm1)
 
-         if (ifpod(2)) call copy(t4,ts0(1,j),n)
+         if (ifpod(2)) call copy(t4,ts0(1,j,1),n)
 
          ttmp = time
          itmp = istep
@@ -457,7 +457,7 @@ subroutine snap_analysis
             s=-ut(i)
             ss = 0
             if (ifpod(2)) then
-               call add2s2(t4,tb(1,i),s,n)
+               call add2s2(t4,tb(1,i,1),s,n)
                do ii=1,n
                   ss=ss+bm1(ii,1,1,1)*(t4(ii)*t4(ii))
                enddo
@@ -503,7 +503,7 @@ subroutine tj_analysis
       npts=128
 
       do i=0,nb
-         call sol_intp_xline_qoi(ub(1,i),vb(1,i),wb(1,i),tb(1,i),
+         call sol_intp_xline_qoi(ub(1,i),vb(1,i),wb(1,i),tb(1,i,1),
      $                           2.5,0.,128,i)
       enddo
 
diff --git a/code/aux.f b/code/aux.f
index 9954083e..b57e821a 100644
--- a/code/aux.f
+++ b/code/aux.f
@@ -28,7 +28,7 @@ subroutine recont(tt,coef)
       call rzero(tt,n)
 
       do i=0,nb
-         call add2s2(tt,tb(1,i),coef(i),n)
+         call add2s2(tt,tb(1,i,1),coef(i),n)
       enddo
 
       return
@@ -51,7 +51,7 @@ subroutine recont_rms(tt)
 
       do j=0,nb
       do i=0,nb
-         call col3(tbt,tb(1,i),tb(1,j),n)
+         call col3(tbt,tb(1,i,1),tb(1,j,1),n)
          call add2s2(tt,tbt,ut2a(1+i+(nb+1)*j),n)
       enddo
       enddo
@@ -371,7 +371,7 @@ subroutine add_sol(vx,vy,vz,pr,t,ux,uy,uz,pp,tt)
       call opadd2(vx,vy,vz,ux,uy,uz)
       call add2(pr,pp,lx2*ly2*lz2*nelv)
 
-      do idim=1,1
+      do idim=1,ldimt
          call add2(t(1,idim),tt(1,idim),lx1*ly1*lz1*nelt)
       enddo
 
@@ -874,10 +874,10 @@ subroutine dump_sfld
          call opzero(ux1,uy1,uz1)
          do j=0,nb
          do i=0,nb
-            call admcol3(ux1,ub(1,i),tb(1,j),uuta(1+i+(nb+1)*j),n)
-            call admcol3(uy1,vb(1,i),tb(1,j),uuta(1+i+(nb+1)*j),n)
+            call admcol3(ux1,ub(1,i),tb(1,j,1),uuta(1+i+(nb+1)*j),n)
+            call admcol3(uy1,vb(1,i),tb(1,j,1),uuta(1+i+(nb+1)*j),n)
             if (ldim.eq.3)
-     $         call admcol3(uz1,wb(1,i),tb(1,j),uuta(1+i+(nb+1)*j),n)
+     $         call admcol3(uz1,wb(1,i),tb(1,j,1),uuta(1+i+(nb+1)*j),n)
          enddo
          enddo
          call outpost(ux1,uy1,uz1,pr,tt,'tmn')
@@ -1418,7 +1418,7 @@ subroutine k_mean(k,nsu,nsp,nst,fn,seed)
                label(i) = j
             endif
          enddo
-         call copy(cent_fld(1,i),ts0(1,label(i)),n)
+         call copy(cent_fld(1,i),ts0(1,label(i),1),n)
       enddo
 
       ! minimize distortion measure
@@ -1428,7 +1428,7 @@ subroutine k_mean(k,nsu,nsp,nst,fn,seed)
          call rzero(rnk,ls*k)
          do i=1,ls
             do j=1,k
-              call sub3(tmp(1,j),ts0(1,i),cent_fld(1,j),n)
+              call sub3(tmp(1,j),ts0(1,i,1),cent_fld(1,j),n)
               dist(j) = glsc2(tmp(1,j),tmp(1,j),n)
             enddo
             write(6,*)ls,minloc(dist),sample(i)
@@ -1459,7 +1459,7 @@ subroutine k_mean(k,nsu,nsp,nst,fn,seed)
          do i=1,k
             call rzero(cent_fld,n*k)
             do j=1,ls
-               call add2s2(cent_fld(1,i),ts0(1,j),rnk(j,i),n)
+               call add2s2(cent_fld(1,i),ts0(1,j,1),rnk(j,i),n)
             enddo
             call cmult(cent_fld(1,i),1./num_sc(i),n)
             centroid(i) = glsc2(sample,rnk(1,i),ls)/num_sc(i)
@@ -1476,7 +1476,7 @@ subroutine k_mean(k,nsu,nsp,nst,fn,seed)
                   label(i) = j
                endif
             enddo
-            call copy(cent_fld(1,i),ts0(1,label(i)),n)
+            call copy(cent_fld(1,i),ts0(1,label(i),1),n)
          enddo
          call c_distortion_measure(obj_f,cent_fld,rnk,k)
          write(6,*)kk,obj_f,'distortion measure M'
@@ -1534,7 +1534,7 @@ subroutine c_distortion_measure(obj_f,cent_fld,rnk,k)
       do i=1,ls
          do j=1,k
             if (abs(rnk(i,j)-1).le.1e-8) then
-               call sub3(tmp(1,j),ts0(1,i),cent_fld(1,j),n)
+               call sub3(tmp(1,j),ts0(1,i,1),cent_fld(1,j),n)
                dist(j) = glsc2(tmp(1,j),tmp(1,j),n)
                obj_f = obj_f + dist(j)
             endif
@@ -1581,8 +1581,8 @@ subroutine projtoprerb(nocp,uuk,ttk)
       do j=1,ns
       ttk(0,j) = 1.
       do i=1,nocp
-         ww=sip(tb(1,i),tb(1,i))
-         vv=sip(tb(1,i),ts0(1,j))
+         ww=sip(tb(1,i,1),tb(1,i,1))
+         vv=sip(tb(1,i,1),ts0(1,j,1))
          ttk(i,j) = vv/ww
       enddo
       enddo
@@ -1597,7 +1597,7 @@ subroutine projtoprerb(nocp,uuk,ttk)
             call sub2(us0(1,2,i),vy,n)
             if (ldim.eq.3) call sub2(us0(1,ldim,i),vz,n)
          endif
-         if (ifrom(2)) call sub2(ts0(1,i),t,n)
+         if (ifrom(2)) call sub2(ts0(1,i,1),t,n)
       enddo
 
       if (nio.eq.0) write (6,*) 'exiting projtoprerb'
@@ -1755,7 +1755,7 @@ subroutine recont_wo0(tt,coef,nocp)
       call rzero(tt,n)
 
       do i=1,nocp
-         call add2s2(tt,tb(1,i),coef(i),n)
+         call add2s2(tt,tb(1,i,1),coef(i),n)
       enddo
 
       return
diff --git a/code/dump.f b/code/dump.f
index b0f19e05..4d677bc9 100644
--- a/code/dump.f
+++ b/code/dump.f
@@ -252,7 +252,7 @@ subroutine dump_all
          time=i
          itmp=i
          ifxyo=(i.eq.0)
-         call outpost(ub(1,i),vb(1,i),wb(1,i),pb(1,i),tb(1,i),'bas')
+         call outpost(ub(1,i),vb(1,i),wb(1,i),pb(1,i),tb(1,i,1),'bas')
       enddo
 
       istep=itmp
@@ -325,12 +325,14 @@ subroutine dump_bas
       parameter (lt=lx1*ly1*lz1*lelt)
 
       common /dumpglobal/ wk1(lcloc),wk2(lcloc)
+      common /romdbas/ tmp(lt,ldimt)
 
       logical iftmp1,iftmp2,iftmp3
 
       call nekgsync
       dbas_time=dnekclock()
 
+      n=lx1*ly1*lz1*lelt
       ttmp=time
       itmp=istep
 
@@ -348,7 +350,10 @@ subroutine dump_bas
          time=i
          itmp=i
          ifxyo=(i.eq.0)
-         call outpost(ub(1,i),vb(1,i),wb(1,i),pb(1,i),tb(1,i),'bas')
+         do j=1,ldimt
+            call copy(tmp(1,j),tb(1,i,j),n)
+         enddo
+         call outpost2(ub(1,i),vb(1,i),wb(1,i),pb(1,i),tmp,ldimt,'bas')
       enddo
 
       istep=itmp
@@ -451,7 +456,7 @@ subroutine dump_snaps
 
       do i=1,ns
          call outpost(us0(1,1,i),us0(1,2,i),us0(1,ldim,i),
-     $      pr,ts0(1,i),'sna')
+     $      pr,ts0(1,i,1),'sna')
          ifxyo=.false.
       enddo
 
diff --git a/code/ei.f b/code/ei.f
index 3a3978f4..2948428b 100644
--- a/code/ei.f
+++ b/code/ei.f
@@ -502,7 +502,7 @@ subroutine setr_poisson(rhs)
 
       do i=1,nb
 c        rhs(i)=wl2sip(qq,tb(1,i))
-         rhs(i)=glsc2(qq,tb(1,i),n)
+         rhs(i)=glsc2(qq,tb(1,i,1),n)
       enddo
 
       return
@@ -610,7 +610,7 @@ subroutine set_residual
       do i=0,nb
          ! setup rhs for temperature representator
          ifield=2
-         call copy(wk4,tb(1,i),n)
+         call copy(wk4,tb(1,i,1),n)
          call axhelm(riesz_rt(1,l2),wk4,ones,zeros,1,1)
          call cmult(riesz_rt(1,l2),param(8),n)
          call chsign(riesz_rt(1,l2),n) 
@@ -622,7 +622,7 @@ subroutine set_residual
          ifield=1
          call opcopy(riesz_ru(1,1,l1),riesz_ru(1,2,l1)
      $               ,riesz_ru(1,ldim,l1)
-     $               ,tb(1,i),zeros,zeros)
+     $               ,tb(1,i,1),zeros,zeros)
          call opcolv(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
      $               riesz_ru(1,ldim,l1),bm1)
          call opchsgn(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
@@ -635,7 +635,7 @@ subroutine set_residual
          ifield=1
          call opcopy(riesz_ru(1,1,l1),riesz_ru(1,2,l1)
      $               ,riesz_ru(1,ldim,l1)
-     $               ,zeros,tb(1,i),zeros)
+     $               ,zeros,tb(1,i,1),zeros)
          call opcolv(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
      $               riesz_ru(1,ldim,l1),bm1)
          call opchsgn(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
@@ -659,7 +659,7 @@ subroutine set_residual
      $                   riesz_ru(1,ldim,l1))
 
          ifield=2
-            call convect_new(riesz_rt(1,l2),tb(1,i),.false.,
+            call convect_new(riesz_rt(1,l2),tb(1,i,1),.false.,
      $                       ub(1,j),vb(1,j),wb(1,j),.false.)
             call chsign(riesz_rt(1,l2),n) 
             l1=l1+1
@@ -1610,7 +1610,7 @@ subroutine set_residual_unsteady
          ifield=1
          call opcopy(riesz_ru(1,1,l1),riesz_ru(1,2,l1)
      $               ,riesz_ru(1,ldim,l1)
-     $               ,tb(1,i),zeros,zeros)
+     $               ,tb(1,i,1),zeros,zeros)
          call opcolv(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
      $               riesz_ru(1,ldim,l1),bm1)
          call opchsgn(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
@@ -1623,7 +1623,7 @@ subroutine set_residual_unsteady
          ifield=1
          call opcopy(riesz_ru(1,1,l1),riesz_ru(1,2,l1)
      $               ,riesz_ru(1,ldim,l1)
-     $               ,zeros,tb(1,i),zeros)
+     $               ,zeros,tb(1,i,1),zeros)
          call opcolv(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
      $               riesz_ru(1,ldim,l1),bm1)
          call opchsgn(riesz_ru(1,1,l1),riesz_ru(1,2,l1),
@@ -1654,7 +1654,7 @@ subroutine set_residual_unsteady
       l2=1
       do i=0,nb
          ifield=2
-         call copy(riesz_rt(1,l2),tb(1,i),n)
+         call copy(riesz_rt(1,l2),tb(1,i,1),n)
          call col2(riesz_rt(1,l2),bm1,n)
          call chsign(riesz_rt(1,l2),n)
          l2=l2+1
@@ -1662,7 +1662,7 @@ subroutine set_residual_unsteady
       if (nid.eq.0) write(6,*)l2,'lres_t_1'
       do i=0,nb
          ifield=2
-         call copy(wk4,tb(1,i),n)
+         call copy(wk4,tb(1,i,1),n)
          call axhelm(riesz_rt(1,l2),wk4,ones,zeros,1,1)
          call cmult(riesz_rt(1,l2),param(8),n)
          call chsign(riesz_rt(1,l2),n) 
@@ -1672,7 +1672,7 @@ subroutine set_residual_unsteady
       do j=0,nb
          do i=0,nb
          ifield=2
-            call convect_new(riesz_rt(1,l2),tb(1,i),.false.,
+            call convect_new(riesz_rt(1,l2),tb(1,i,1),.false.,
      $                       ub(1,j),vb(1,j),wb(1,j),.false.)
             call chsign(riesz_rt(1,l2),n) 
             l2=l2+1
@@ -1917,7 +1917,7 @@ subroutine c_rieszrd_uns
       call mxm(utj,(nb+1),alphaj,6,coef,1)
       do i=0,nb
          ifield=1
-         call opcopy(wk1,wk2,wk3,tb(1,i),zeros,zeros)
+         call opcopy(wk1,wk2,wk3,tb(1,i,1),zeros,zeros)
          call opcolv(wk1,wk2,wk3,bm1)
          call opchsgn(wk1,wk2,wk3)
 
@@ -1936,7 +1936,7 @@ subroutine c_rieszrd_uns
       call mxm(utj,(nb+1),alphaj,6,coef,1)
       do i=0,nb
          ifield=1
-         call opcopy(wk1,wk2,wk3,zeros,tb(1,i),zeros)
+         call opcopy(wk1,wk2,wk3,zeros,tb(1,i,1),zeros)
          call opcolv(wk1,wk2,wk3,bm1)
          call opchsgn(wk1,wk2,wk3)
 
@@ -1985,7 +1985,7 @@ subroutine c_rieszrd_uns
       call mxm(utj,nb+1,betaj,6,coef,1)
       do i=0,nb
          ifield=2
-         call copy(wk1,tb(1,i),n)
+         call copy(wk1,tb(1,i,1),n)
          call col2(wk1,bm1,n)
          call chsign(wk1,n)
 
@@ -1996,7 +1996,7 @@ subroutine c_rieszrd_uns
       if (nid.eq.0) write(6,*)l2,'lres_t_1'
       do i=0,nb
          ifield=2
-         call copy(wk2,tb(1,i),n)
+         call copy(wk2,tb(1,i,1),n)
          call axhelm(wk1,wk2,ones,zeros,1,1)
          call cmult(wk1,param(8),n)
          call chsign(wk1,n) 
@@ -2012,7 +2012,7 @@ subroutine c_rieszrd_uns
       do j=0,nb
          do i=0,nb
             ifield=2
-            call convect_new(wk1,tb(1,i),.false.,
+            call convect_new(wk1,tb(1,i,1),.false.,
      $                       ub(1,j),vb(1,j),wb(1,j),.false.)
             call chsign(wk1,n) 
 
@@ -2312,7 +2312,7 @@ subroutine resid_in_time(msg)
          call mxm(utj,nb+1,betaj,8,coef,1)
          do i=0,nb
             ifield=2
-            call copy(wk1,tb(1,i),n)
+            call copy(wk1,tb(1,i,1),n)
             call col2(wk1,bm1,n)
             call chsign(wk1,n)
 
@@ -2369,7 +2369,7 @@ subroutine resid_in_diffusion(msg)
       elseif (msg.eq.'tmp') then
          do i=0,nb
             ifield=2
-            call copy(wk2,tb(1,i),n)
+            call copy(wk2,tb(1,i,1),n)
             call axhelm(wk1,wk2,ones,zeros,1,1)
             call cmult(wk1,param(8),n)
             call chsign(wk1,n)
@@ -2383,7 +2383,7 @@ subroutine resid_in_diffusion(msg)
             if(nid.eq.0) write(6,*)'helm is on'
             do i=0,nb
                ifield=2
-               call copy(wk1,tb(1,i),n)
+               call copy(wk1,tb(1,i,1),n)
                call col2(wk1,bm1,n)
                call chsign(wk1,n)
 
@@ -2449,7 +2449,7 @@ subroutine resid_in_advec(msg)
          do j=0,nb
             do i=0,nb
                ifield=2
-               call convect_new(wk1,tb(1,i),.false.,
+               call convect_new(wk1,tb(1,i,1),.false.,
      $                          ub(1,j),vb(1,j),wb(1,j),.false.)
                call chsign(wk1,n)
 
@@ -2489,7 +2489,7 @@ subroutine resid_in_buoy
       call mxm(utj,(nb+1),alphaj,8,coef,1)
       do i=0,nb
          ifield=1
-         call opcopy(wk1,wk2,wk3,tb(1,i),zeros,zeros)
+         call opcopy(wk1,wk2,wk3,tb(1,i,1),zeros,zeros)
          call opcolv(wk1,wk2,wk3,bm1)
          call opchsgn(wk1,wk2,wk3)
 
@@ -2506,7 +2506,7 @@ subroutine resid_in_buoy
       call mxm(utj,(nb+1),alphaj,8,coef,1)
       do i=0,nb
          ifield=1
-         call opcopy(wk1,wk2,wk3,zeros,tb(1,i),zeros)
+         call opcopy(wk1,wk2,wk3,zeros,tb(1,i,1),zeros)
          call opcolv(wk1,wk2,wk3,bm1)
          call opchsgn(wk1,wk2,wk3)
 
diff --git a/code/filter.f b/code/filter.f
index 086e0580..49a9728c 100644
--- a/code/filter.f
+++ b/code/filter.f
@@ -125,7 +125,7 @@ subroutine set_les_imp(fles1,fles2)
          endif
 
          if (ifrom(2)) then
-            call copy(t,tb(1,i),nt)
+            call copy(t,tb(1,i,1),nt)
             call add2(t,tb,nt)
          endif
 
diff --git a/code/mpar.f b/code/mpar.f
index 429c7e48..35ea29d4 100644
--- a/code/mpar.f
+++ b/code/mpar.f
@@ -533,7 +533,7 @@ subroutine bcastmpar
 
       call bcast(ifrecon,lsize)
 
-      do i=1,2
+      do i=0,ldimt1
          call bcast(ifpod(i),lsize)
          call bcast(ifrom(i),lsize)
       enddo
diff --git a/code/pod.f b/code/pod.f
index 4dccd1ad..1afc6a1f 100644
--- a/code/pod.f
+++ b/code/pod.f
@@ -37,7 +37,8 @@ subroutine setbases
             call opcopy(ub,vb,wb,uic,vic,wic)
          endif
          if (ifrom(2)) then
-            call pod(tb(1,1),eval,ug,ts0,1,ips,nb,ns,ifpb,'ops/gt  ')
+            call pod(tb(1,1,1),eval,ug,ts0(1,1,1),1,ips,
+     $               nb,ns,ifpb,'ops/gt  ')
             if (.not.ifcomb.and.ifpb) call snorm(tb)
          endif
 
@@ -98,10 +99,10 @@ subroutine setbases
             nt=lx1*ly1*lz1*nelt
             do i=0,nb
                call rzero(upup,nv)
-               call rzero(tb(1,i+nb+1),nt)
+               call rzero(tb(1,i+nb+1,1),nt)
 
                call evalcflds(
-     $            upup,uvwb(1,1,0),tb(1,i),1,1,.true.)
+     $            upup,uvwb(1,1,0),tb(1,i,1),1,1,.true.)
 
                call col2(upup,tmask,nt)
                call dssum(upup,lx1,ly1,lz1)
@@ -110,7 +111,7 @@ subroutine setbases
                sc=1./sqrt(glsc3(upup,upup,bm1,nv))
 
                call cmult(upup,sc,nv)
-               call copy(tb(1,i+nb+1),upup,nv)
+               call copy(tb(1,i+nb+1,1),upup,nv)
             enddo
          endif
 
@@ -172,10 +173,10 @@ subroutine setbases
             do i=0,nb
 
                call rzero(upup,nv)
-               call rzero(tb(1,i+nb+1),nt)
+               call rzero(tb(1,i+nb+1,1),nt)
 
                call evalcflds(
-     $            upup,uvwb(1,1,i),tb(1,i),1,1,.true.)
+     $            upup,uvwb(1,1,i),tb(1,i,1),1,1,.true.)
 
                call col2(upup,tmask,nt)
                call dssum(upup,lx1,ly1,lz1)
@@ -184,7 +185,7 @@ subroutine setbases
                sc=1./sqrt(glsc3(upup,upup,bm1,nv))
 
                call cmult(upup,sc,nv)
-               call copy(tb(1,i+nb+1),upup,nv)
+               call copy(tb(1,i+nb+1,1),upup,nv)
             enddo
          endif
 
@@ -270,10 +271,10 @@ subroutine setbases
             nt=lx1*ly1*lz1*nelt
             do i=0,nb
                call rzero(upup,nt)
-               call rzero(tb(1,i+nb+1),nt)
+               call rzero(tb(1,i+nb+1,1),nt)
 
                call evalcflds(
-     $            upup,uvwb(1,1,0),tb(1,i),1,1,.true.)
+     $            upup,uvwb(1,1,0),tb(1,i,1),1,1,.true.)
 
                call col2(upup,tmask,nt)
                call dssum(upup,lx1,ly1,lz1)
@@ -282,15 +283,15 @@ subroutine setbases
                sc=1./sqrt(glsc3(upup,upup,bm1,nt))
 
                call cmult(upup,sc,nt)
-               call copy(tb(1,i+nb+1),upup,nt)
+               call copy(tb(1,i+nb+1,1),upup,nt)
             enddo
             do i=1,nb
 
                call rzero(upup,nt)
-               call rzero(tb(1,i+2*nb+1),nt)
+               call rzero(tb(1,i+2*nb+1,1),nt)
 
                call evalcflds(
-     $            upup,uvwb(1,1,i),tb(1,i),1,1,.true.)
+     $            upup,uvwb(1,1,i),tb(1,i,1),1,1,.true.)
 
                call col2(upup,tmask,nt)
                call dssum(upup,lx1,ly1,lz1)
@@ -299,7 +300,7 @@ subroutine setbases
                sc=1./sqrt(glsc3(upup,upup,bm1,nt))
 
                call cmult(upup,sc,nv)
-               call copy(tb(1,i+2*nb+1),upup,nt)
+               call copy(tb(1,i+2*nb+1,1),upup,nt)
             enddo
          endif
 
diff --git a/code/qoi.f b/code/qoi.f
index 637bff37..fae302a2 100644
--- a/code/qoi.f
+++ b/code/qoi.f
@@ -272,11 +272,11 @@ subroutine cnuss_setup
          else
             do j=0,nb
                do i=0,nb
-                  call ctbulk_num(
-     $               tbn(i+j*(nb+1),0),ub(1,i),vb(1,i),wb(1,i),tb(1,j))
+                  call ctbulk_num(tbn(i+j*(nb+1),0),
+     $               ub(1,i),vb(1,i),wb(1,i),tb(1,j,1))
                enddo
                call ctbulk_den(tbd(j),ub(1,j),vb(1,j),wb(1,j))
-               call ctsurf(tsa(j),tb(1,j))
+               call ctsurf(tsa(j),tb(1,j,1))
             enddo
 
             call dump_serial(tsa,nb+1,'qoi/tsa ',nid)
@@ -301,7 +301,7 @@ subroutine cnuss_setup
          endif
       else if (inus.eq.2) then
          do i=0,nb
-            call gradm1(tx,ty,tz,tb(1,i))
+            call gradm1(tx,ty,tz,tb(1,i,1))
 
             eps=1.e-6
             ta=0.
@@ -339,7 +339,7 @@ subroutine cnuss_setup
       else if (inus.eq.3) then
          tbn(0,0)=2.
          do j=0,nb
-            call ctsurf3(tsa(j),tb(1,j))
+            call ctsurf3(tsa(j),tb(1,j,1))
          enddo
 
          do i=0,nb
@@ -349,7 +349,7 @@ subroutine cnuss_setup
          call rone(ones,lx1*ly1*lz1*nelt)
          if (rmode.ne.'ON ') then
          do i=0,nb
-            call gradm1(tx,ty,tz,tb(1,i))
+            call gradm1(tx,ty,tz,tb(1,i,1))
 
             eps=1.e-3
             ta=0.
@@ -387,7 +387,7 @@ subroutine cnuss_setup
          call rone(ones,lx1*ly1*lz1*nelt)
          if (rmode.ne.'ON ') then
          do i=0,nb
-            call gradm1(tx,ty,tz,tb(1,i))
+            call gradm1(tx,ty,tz,tb(1,i,1))
 
             eps=1.e-3
             ta=0.
@@ -426,7 +426,7 @@ subroutine cnuss_setup
          iobj=1
 
          do i=0,nb
-            call gradm1(tx,ty,tz,tb(1,i))
+            call gradm1(tx,ty,tz,tb(1,i,1))
             a=0.
             s=0.
 
diff --git a/code/read.f b/code/read.f
index e67b4dfd..ccacec9b 100644
--- a/code/read.f
+++ b/code/read.f
@@ -135,9 +135,9 @@ subroutine loadbases
 
       if (ifexist) then
          nn=nb+1
-         ifreads(1)=ifrom(1)
-         ifreads(2)=ifrom(0)
-         ifreads(3)=ifrom(2)
+         do i=0,ldimt1
+            ifreads(i)=ifrom(i)
+         enddo
 
          call read_fields(
      $      us0,prs,ts0,nn,0,ifreads,tk,'bas.list ',.false.)
@@ -146,7 +146,7 @@ subroutine loadbases
             if (ifrom(0)) call copy(pb(1,i),prs(1,i+1),n2)
             if (ifrom(1)) call opcopy(ub(1,i),vb(1,i),wb(1,i),
      $                        us0(1,1,i+1),us0(1,2,i+1),us0(1,ldim,i+1))
-            if (ifrom(2)) call copy(tb(1,i),ts0(1,i+1),n)
+            if (ifrom(2)) call copy(tb(1,i,1),ts0(1,i+1,1),n)
          enddo
          if (nn.lt.nb) call exitti(
      $   'number of files in bas.list fewer than nb$',nb-nn)
@@ -166,7 +166,7 @@ subroutine loadbases
             call restart_filen(fname,11+len)
             if (ifrom(0)) call copy(pb(1,i),pr,n2)
             if (ifrom(1)) call opcopy(ub(1,i),vb(1,i),wb(1,i),vx,vy,vz)
-            if (ifrom(2)) call copy(tb(1,i),t,n)
+            if (ifrom(2)) call copy(tb(1,i,1),t,n)
          enddo
       endif
 
@@ -199,13 +199,13 @@ subroutine read_fields(usave,psave,tsave,ns,nskp,ifread,tk,fn,ifa)
       parameter (lt=lx1*ly1*lz1*lelt)
       parameter (lt2=lx2*ly2*lz2*lelt)
 
-      real usave(lt,ldim,1),psave(lt2,1),tsave(lt,1)
+      real usave(lt,ldim,1),psave(lt2,1),tsave(lt,1,ldimt)
       real tk(1)
 
       character*128 fn
       character*128 fnlint
 
-      logical ifa,ifread(3)
+      logical ifa,ifread(0:ldimt1)
 
       common /scrk2/ t4(lt),t5(lt),t6(lt)
 
@@ -263,12 +263,15 @@ subroutine read_fields(usave,psave,tsave,ns,nskp,ifread,tk,fn,ifa)
                if (ldim.eq.3) call add2col2(uvms,vz,t,n)
             endif
 
+            if (ifread(0)) call copy(psave(1,ip),pr,n2)
             if (ifread(1))
      $         call opcopy(usave(1,1,ip),usave(1,2,ip),usave(1,ldim,ip),
      $                    vx,vy,vz)
 
-            if (ifread(2)) call copy(psave(1,ip),pr,n2)
-            if (ifread(3)) call copy(tsave(1,ip),t,n)
+            do j=1,ldimt
+               idx=j+1
+               if (ifread(idx)) call copy(tsave(1,ip,j),t(1,1,1,1,j),n)
+            enddo
          else
             goto 999
          endif
@@ -870,7 +873,7 @@ subroutine loadpbases(nsave)
             call restart_filen(fname,11+len)
             if (ifrom(0)) call copy(pb(1,i),pr,n2)
             if (ifrom(1)) call opcopy(ub(1,i),vb(1,i),wb(1,i),vx,vy,vz)
-            if (ifrom(2)) call copy(tb(1,i),t,n)
+            if (ifrom(2)) call copy(tb(1,i,1),t,n)
          enddo
       endif
 
@@ -946,7 +949,7 @@ subroutine get_p_rb(nocp,nsu,nsp,nst,ttk,fn)
      $                    vx,vy,vz)
 
             if (icount.le.nsp) call copy(pb(1,ip),pr,n2)
-            if (icount.le.nst) call copy(tb(1,ip),t,n)
+            if (icount.le.nst) call copy(tb(1,ip,1),t,n)
          else
             goto 999
          endif
diff --git a/code/riesz.f b/code/riesz.f
index 1216115a..24d00101 100644
--- a/code/riesz.f
+++ b/code/riesz.f
@@ -16,12 +16,12 @@ subroutine set_xi_poisson
         call exitti('(set_xi_poisson) ifield.eq.1 not supported...$',nb)
       else
          if (ips.eq.'L2 ') then
-            call set_xi_a(xi(1,l),tb(1,1),ones,1,nb,2)
+            call set_xi_a(xi(1,l),tb(1,1,1),ones,1,nb,2)
             l=l+nb
             call set_xi_b(xi(1,l),qq,1,1,2)
             l=l+1
             do i=1,nb
-               call set_gradn(wk,tb(1,i))
+               call set_gradn(wk,tb(1,i,1))
                call set_surf(xi(1,l),wk)
                call binv1_nom(xi(1,l))
                l=l+1
@@ -90,7 +90,7 @@ subroutine set_xi_ad
                   do j=0,nb
                      call opcopy(vx,vy,vz,ub(1,j),vb(1,j),wb(1,j))
                      do i=0,nb
-                        call conv1d(xi(1,l),tb(1,i))
+                        call conv1d(xi(1,l),tb(1,i,1))
                         l=l+1
                      enddo
                   enddo
@@ -178,7 +178,7 @@ subroutine set_xi_ns
                call exitti('Buoyancy in EI disabled for now...l',1)
                do i=0,nb
                   call opcopy(wk1,wk2,wk3,gx,gy,gz)
-                  call opcolv(wk1,wk2,wk3,tb(1,i))
+                  call opcolv(wk1,wk2,wk3,tb(1,i,1))
                   call invcol2(wk1,bm1,n)
                   call invcol2(wk2,bm1,n)
                   if (ldim.eq.3) call invcol2(wk3,bm1,n)
diff --git a/code/rom.f b/code/rom.f
index 55a7639e..6c82893f 100644
--- a/code/rom.f
+++ b/code/rom.f
@@ -425,7 +425,8 @@ subroutine setqoi
          m=0
          if (nid.eq.0) m=nintp
          do i=0,nb
-            call gfldi(tbintp(1+nintp*i),tb(1,i),xintp,yintp,zintp,m,1)
+            call gfldi(tbintp(1+nintp*i),tb(1,i,1),
+     $                 xintp,yintp,zintp,m,1)
          enddo
          call dump_serial(xintp,nintp*(nb+1),'qoi/xintp',nid)
          call dump_serial(yintp,nintp*(nb+1),'qoi/yintp',nid)
@@ -1052,9 +1053,9 @@ subroutine mor_init_fields
       if (rmode.eq.'ALL'.or.rmode.eq.'OFF'.or.rmode.eq.'AEQ') then
          fname1='file.list '
 
-         ifreads(1)=ifrom(1)
-         ifreads(2)=ifrom(0)
-         ifreads(3)=ifrom(2)
+         do i=0,ldimt1
+            ifreads(i)=ifrom(i)
+         enddo
 
          call read_fields(
      $      us0,prs,ts0,ns,nskip,ifreads,timek,fname1,.true.)
@@ -1124,7 +1125,7 @@ subroutine mor_init_fields
                   call sub2(us0(1,2,i),vb,n)
                   if (ldim.eq.3) call sub2(us0(1,ldim,i),wb,n)
                endif
-               if (ifrom(2)) call sub2(ts0(1,i),tb,n)
+               if (ifrom(2)) call sub2(ts0(1,i,1),tb,n)
             enddo
             call sub2(uavg,ub,n)
             call sub2(vavg,vb,n)
@@ -1281,7 +1282,7 @@ subroutine setc(cl,fname)
                call copy(u2va(1,i),u2v,nd)
                if (ldim.eq.3) call copy(u3va(1,i),u3v,nd)
             else
-               call setcnv_u1(tb(1,i))
+               call setcnv_u1(tb(1,i,1))
                call copy(u1va(1,i),u1v,nd)
             endif
          enddo
@@ -1295,7 +1296,7 @@ subroutine setc(cl,fname)
      $                  ub(1,j),vb(1,j),wb(1,j))
                      call convect_axis(cu,ldim,ux,uy,uz,wku)
                   else
-                     call convect_axis(cu,1,ux,uy,uz,tb(1,j))
+                     call convect_axis(cu,1,ux,uy,uz,tb(1,j,1))
                   endif
                else
                   if (ifield.eq.1) then
@@ -1315,7 +1316,7 @@ subroutine setc(cl,fname)
                      cel=op_glsc2_wt(ub(1,i),vb(1,i),wb(1,i),
      $                  cu(1,1),cu(1,2),cu(1,ldim),ones)
                   else
-                     cel=glsc2(tb(1,i),cu,n)
+                     cel=glsc2(tb(1,i,1),cu,n)
                   endif
                   call setc_local(cl,cel,ic1,ic2,jc1,jc2,kc1,kc2,i,j,k)
                   if (nid.eq.0) write (100,*) cel
@@ -1376,9 +1377,9 @@ subroutine seta(a,a0,fname)
      $                           ub(1,j),vb(1,j),wb(1,j))
                else
                   if (nelgt.ne.nelvt) then
-                     a0(i,j)=h10sip_vd(tb(1,i),tb(1,j),vdm1)
+                     a0(i,j)=h10sip_vd(tb(1,i,1),tb(1,j,1),vdm1)
                   else
-                     a0(i,j)=h10sip(tb(1,i),tb(1,j))
+                     a0(i,j)=h10sip(tb(1,i,1),tb(1,j,1))
                   endif
                endif
             enddo
@@ -1494,9 +1495,9 @@ subroutine setb(b,b0,fname)
      $                           ub(1,j),vb(1,j),wb(1,j))
                else
                   if (nelgt.eq.nelvt) then
-                     b0(i,j)=wl2sip(tb(1,i),tb(1,j))
+                     b0(i,j)=wl2sip(tb(1,i,1),tb(1,j,1))
                   else
-                     b0(i,j)=wl2sip_vd(tb(1,i),tb(1,j),brhom1)
+                     b0(i,j)=wl2sip_vd(tb(1,i,1),tb(1,j,1),brhom1)
                   endif
                endif
             enddo
@@ -1548,7 +1549,7 @@ subroutine setae(a,fname)
                   a(i,j,k)=h10vip_vd(ub(1,i),vb(1,i),wb(1,i),
      $                             ub(1,j),vb(1,j),wb(1,j),udfld(1,1,k))
                else
-                  a(i,j,k)=h10sip_vd(tb(1,i),tb(1,j),tdfld(1,k))
+                  a(i,j,k)=h10sip_vd(tb(1,i,1),tb(1,j,1),tdfld(1,k))
                endif
             enddo
             nio=nid
@@ -1614,8 +1615,8 @@ subroutine setu
 
          if (ifrom(2)) then
             ifield=2
-            call sub2(tic,tb,n)
-            call ps2b(ut,tic,tb)
+            call sub2(tic,tb(1,0,1),n)
+            call ps2b(ut,tic,tb(1,0,1))
             do i=0,nb
                if (nio.eq.0) write (6,*) 'ut',ut(i)
             enddo
@@ -1703,7 +1704,7 @@ subroutine setf
 
       if (ifsource.and.ifrom(2)) then ! assume qq has mass
          do i=1,nb
-            rq(i)=glsc2(qq,tb(1,i),n)
+            rq(i)=glsc2(qq,tb(1,i,1),n)
             if (nio.eq.0) write (6,*) rq(i),i,'rq'
          enddo
          call copy(wk1,qq,n)
@@ -1711,7 +1712,7 @@ subroutine setf
          call outpost(vx,vy,vz,pavg,wk1,'qqq')
          if (ifsrct) then
             do i=1,nb
-               rqt(i)=glsc2(qqxyz,tb(1,i),n)
+               rqt(i)=glsc2(qqxyz,tb(1,i,1),n)
                if (nio.eq.0) write (6,*) rqt(i),i,'rqt'
             enddo
             call copy(wk1,qqxyz,n)
@@ -1774,7 +1775,7 @@ subroutine setfluc(fv,ftt,fname)
       do i=1,nbavg
          fv(i,j)=wl2vip(ub(1,i),vb(1,i),wb(1,i),
      $                  flucv(1,1,j),flucv(1,2,j),flucv(1,ldim,j))
-         ftt(i,j)=wl2sip(tb(1,i),fluct(1,j))
+         ftt(i,j)=wl2sip(tb(1,i,1),fluct(1,j))
       enddo
       enddo
 
@@ -1871,9 +1872,9 @@ subroutine setbut(bx,by,bz)
       if (rmode.eq.'ALL'.or.rmode.eq.'OFF'.or.rmode.eq.'AEQ') then
          do j=0,nb
          do i=0,nb
-            bx(i,j)=glsc3(ub(1,i),tb(1,j),bm1,n)
-            by(i,j)=glsc3(vb(1,i),tb(1,j),bm1,n)
-            if (ldim.eq.3) bz(i,j)=glsc3(wb(1,i),tb(1,j),bm1,n)
+            bx(i,j)=glsc3(ub(1,i),tb(1,j,1),bm1,n)
+            by(i,j)=glsc3(vb(1,i),tb(1,j,1),bm1,n)
+            if (ldim.eq.3) bz(i,j)=glsc3(wb(1,i),tb(1,j,1),bm1,n)
          enddo
          enddo
          call dump_serial(bx,(nb+1)**2,'ops/buxt ',nid)