diff --git a/actions/actiontsc.h b/actions/actiontsc.h
index 899d720..1a071bd 100644
--- a/actions/actiontsc.h
+++ b/actions/actiontsc.h
@@ -24,10 +24,20 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
  auto [tindexes, err] = GetTIndexes(ds, args, pars);
  if(err.Exist()) return err;
 
- if(!args.contains("var")) return "Variable not specified";
- MString vname = args.at("var");
- if(!ds.CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
- pars.SetParameter("variable", vname);
+ if(!args.contains("var")) return "Variables not specified";
+ auto                 vlist = michlib::Split_on_words(args.at("var"), " ,", false);
+ std::vector<MString> vnames{std::make_move_iterator(std::begin(vlist)), std::make_move_iterator(std::end(vlist))};
+ {
+  std::set<MString> used;
+  for(const auto& vname: vnames)
+  {
+   if(used.contains(vname)) return "Duplicate variable " + vname + " in list " + args.at("var");
+   if(!ds.CheckVar(vname)) return "Variable " + vname + " not exists in this dataset";
+   used.insert(vname);
+  }
+ }
+
+ pars.SetParameter("variables", args.at("var"));
 
  std::unique_ptr<const BaseParameters> sourcepars;
  if constexpr(ParametersSupported<D>)
@@ -47,11 +57,15 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
  }
  auto p = sourcepars.get();
 
- auto data = Read(ds, vname, p, tindexes);
- if(!data) return "Can't read data";
- if(!data.Unit().Exist()) michlib::errmessage("Unknown measurement unit!");
- if(!data.StandartName().Exist()) data.SetStandartName(NCFuncs::Name2StName(vname));
- if(!data.LongName().Exist()) data.SetStandartName(NCFuncs::Name2LongName(vname));
+ auto data = Read(ds, vnames, p, tindexes);
+ if(data.size() != vnames.size()) return "Can't read data";
+
+ for(size_t i = 0; i < vnames.size(); i++)
+ {
+  if(!data[i].Unit().Exist()) michlib::errmessage("Unknown measurement unit for variable " + vnames[i] + "!");
+  if(!data[i].StandartName().Exist()) data[i].SetStandartName(NCFuncs::Name2StName(vnames[i]));
+  if(!data[i].LongName().Exist()) data[i].SetStandartName(NCFuncs::Name2LongName(vnames[i]));
+ }
 
  MString name   = args.contains("out") ? args.at("out") : "out.bin";
  MString outfmt = args.contains("format") ? args.at("format") : (GetExt(name) == "nc" ? "nc" : "bin");
@@ -60,16 +74,17 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
  {
   BFileW fw;
 
-  fw.Create(name, 3);
-  fw.SetColumnName(1, "Longitude");
-  fw.SetColumnName(2, "Latitude");
-  fw.SetColumnName(3, vname + ", " + (data.Unit().Exist() ? data.Unit() : "unknown"));
+  size_t ic = 0;
+  fw.Create(name, 2 + data.size());
+  fw.SetColumnName(++ic, "Longitude");
+  fw.SetColumnName(++ic, "Latitude");
+  for(size_t i = 0; i < data.size(); i++) fw.SetColumnName(++ic, vnames[i] + ", " + (data[i].Unit().Exist() ? data[i].Unit() : "unknown"));
   fw.SetParameters(pars);
-  for(size_t i = 0; i < data.N(); i++)
+  for(size_t r = 0; r < data[0].N(); r++)
   {
-   fw.Write(data.Lon(i));
-   fw.Write(data.Lat(i));
-   fw.Write(data.IsFill(i) ? NAN : data(i));
+   fw.Write(data[0].Lon(r));
+   fw.Write(data[0].Lat(r));
+   for(size_t i = 0; i < data.size(); i++) fw.Write(data[i].IsFill(r) ? NAN : data[i](r));
   }
   fw.Finalize();
   fw.Close();
@@ -84,10 +99,12 @@ template<class D> MString ActionTSC::DoAction(const CLArgs& args, D& ds)
 
   if(args.contains("compress")) compress = args.at("compress").ToInt();
 
-  if(!err.Exist()) err = fw.Create(data, name, args.at("_cmdline"), compress);
-  if(!err.Exist()) err = fw.AddVariable(vname, data.StandartName(), data.LongName(), data.Unit(), data.Comment());
-  if(!err.Exist()) err = fw.WriteGrid(data);
-  if(!err.Exist()) err = fw.WriteVariable(data, vname);
+  if(!err.Exist()) err = fw.Create(data[0], name, args.at("_cmdline"), compress);
+  for(size_t i = 0; i < data.size(); i++)
+   if(!err.Exist()) err = fw.AddVariable(vnames[i], data[i].StandartName(), data[i].LongName(), data[i].Unit(), data[i].Comment());
+  if(!err.Exist()) err = fw.WriteGrid(data[0]);
+  for(size_t i = 0; i < data.size(); i++)
+   if(!err.Exist()) err = fw.WriteVariable(data[i], vnames[i]);
   if(err.Exist()) return err;
   fw.Close();
 
diff --git a/include/actiondep.h b/include/actiondep.h
index 14dd54d..93b1be4 100644
--- a/include/actiondep.h
+++ b/include/actiondep.h
@@ -105,41 +105,50 @@ template<class D> std::pair<TIndex, MString> GetTIndexes(const D& data, const CL
  return {tindexes, ""};
 }
 
-template<class D> ReadType<D> Read(const D& data, const MString& vname, const BaseParameters* p, const TIndex& tindex)
+template<class D> std::vector<ReadType<D>> Read(const D& data, const std::vector<MString>& vnames, const BaseParameters* p, const TIndex& tindex)
 {
  using RT = ReadType<D>;
- size_t ind;
+ size_t          ind;
+ std::vector<RT> out;
  if(tindex.size() == 1)
  {
   ind = tindex[0];
   michlib::message("Time: " + data.Time(ind).ToTString());
-  if constexpr(ReadPSupported<D>)
-   return data.Read(vname, p, ind);
-  else if constexpr(ReadSupported<D>)
-   return data.Read(vname, ind);
+  std::map<MString, RT> cache;
+  for(const auto& vname: vnames)
+  {
+   bool res;
+   if constexpr(ReadPSupported<D>)
+    res = data.Read(vname, cache, p, ind);
+   else if constexpr(ReadSupported<D>)
+    res = data.Read(vname, cache, ind);
+   if(!res) return out;
+  }
+  for(size_t i = 0; i < vnames.size(); i++) out.emplace_back(std::move(cache[vnames[i]]));
  }
  else
  {
-  Averager<RT> out;
-  bool         ok = true;
+  std::vector<Averager<RT>> aver(vnames.size());
   for(size_t i = 0; i < tindex.size(); i++)
   {
-   if(!ok) break;
    ind = tindex[i];
    michlib::message("Time: " + data.Time(ind).ToTString());
-   RT dat;
-   if constexpr(ReadPSupported<D>)
-    dat = data.Read(vname, p, ind);
-   else if constexpr(ReadSupported<D>)
-    dat = data.Read(vname, ind);
-   if(dat)
-    out.Add(std::move(dat));
-   else
-    ok = false;
+   std::map<MString, RT> cache;
+   for(const auto& vname: vnames)
+   {
+    bool res;
+    if constexpr(ReadPSupported<D>)
+     res = data.Read(vname, cache, p, ind);
+    else if constexpr(ReadSupported<D>)
+     res = data.Read(vname, cache, ind);
+    if(!res) return out;
+   }
+   for(size_t i = 0; i < vnames.size(); i++) aver[i].Add(std::move(cache[vnames[i]]));
   }
-  if(ok) return std::move(out.Div());
+  for(size_t i = 0; i < vnames.size(); i++) out.emplace_back(std::move(aver[i].Div()));
+  return out;
  }
- return RT();
+ return out;
 }
 
 template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, size_t ind)
@@ -147,18 +156,15 @@ template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters
  using RT = ReadType<D>;
  using UV = UVData<RT>;
  michlib::message("Time: " + data.Time(ind).ToTString());
- RT u, v;
+ std::map<MString, RT> cache;
+ bool                  res = false;
  if constexpr(ReadPSupported<D>)
- {
-  u = data.Read("u", p, ind);
-  v = data.Read("v", p, ind);
- }
+  res = data.Read("u", cache, p, ind) && data.Read("v", cache, p, ind);
  else if constexpr(ReadSupported<D>)
- {
-  u = data.Read("u", ind);
-  v = data.Read("v", ind);
- }
- return UV(u, v);
+  res = data.Read("u", cache, ind) && data.Read("v", cache, ind);
+ if(!res) return UV();
+
+ return UV(cache.at("u"), cache.at("v"));
 }
 
 template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters* p, const TIndex& tindex)
@@ -178,18 +184,15 @@ template<class D> UVData<ReadType<D>> ReadUV(const D& data, const BaseParameters
    if(!ok) break;
    ind = tindex[i];
    michlib::message("Time: " + data.Time(ind).ToTString());
-   RT u, v;
+   std::map<MString, RT> cache;
+   bool                  res = false;
    if constexpr(ReadPSupported<D>)
-   {
-    u = data.Read("u", p, ind);
-    v = data.Read("v", p, ind);
-   }
+    res = data.Read("u", cache, p, ind) && data.Read("v", cache, p, ind);
    else if constexpr(ReadSupported<D>)
-   {
-    u = data.Read("u", ind);
-    v = data.Read("v", ind);
-   }
-   UV dat(u, v);
+    res = data.Read("u", cache, ind) && data.Read("v", cache, ind);
+   if(!res) return UV();
+
+   UV dat(cache.at("u"), cache.at("v"));
    if(dat)
     out.Add(std::move(dat));
    else
diff --git a/include/layereddata.h b/include/layereddata.h
index c88e272..9e5b984 100644
--- a/include/layereddata.h
+++ b/include/layereddata.h
@@ -135,7 +135,7 @@ class LayeredData: public NCFuncs
 
  std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;
 
- Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
+ bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t i) const;
 
  bool isOk() const { return nc.size() > 0; }
 
diff --git a/include/simple2ddata.h b/include/simple2ddata.h
index 66b9924..527e139 100644
--- a/include/simple2ddata.h
+++ b/include/simple2ddata.h
@@ -18,6 +18,8 @@ class Simple2DData: public BaseData
  {
  }
 
+ Simple2DData CopyGrid() const { return {nx, ny, x0, y0, xstep, ystep}; }
+
  const real& V(size_t ix, size_t iy) const { return V(iy * nx + ix); }
  real&       V(size_t ix, size_t iy) { return V(iy * nx + ix); }
 
diff --git a/include/traits.h b/include/traits.h
index b85dd5e..f19ae71 100644
--- a/include/traits.h
+++ b/include/traits.h
@@ -31,16 +31,16 @@ template<class T>
 concept ParametersSupported = ParametersRequiredRegion<T> || ParametersNotRequiredRegion<T>;
 
 template<class T>
-concept ReadPSupported = requires(T t, const MString& vname, const BaseParameters* ip, size_t i) {
+concept ReadPSupported = requires(T t, const MString& vname, std::map<MString, typename T::Data>& cache, const BaseParameters* ip, size_t i) {
  {
-  t.Read(vname, ip, i)(0)
- } -> std::convertible_to<real>;
+  t.Read(vname, cache, ip, i)
+ } -> std::same_as<bool>;
 };
 template<class T>
-concept ReadSupported = requires(T t, const MString& vname, size_t i) {
+concept ReadSupported = requires(T t, const MString& vname, std::map<MString, typename T::Data>& cache, size_t i) {
  {
-  t.Read(vname, i)(0)
- } -> std::convertible_to<real>;
+  t.Read(vname, cache, i)
+ } -> std::same_as<bool>;
 };
 
 template<class T>
@@ -70,20 +70,14 @@ concept HaveXYStep = requires(T t) {
 
 template<typename...> using void_ = void;
 
-template<class D, bool RP, bool R> struct GetReadType_;
-
-template<class D, bool R> struct GetReadType_<D, true, R>
-{
- using p    = BaseParameters*;
- using type = decltype(D().Read(MString(), p(), 0));
-};
+template<class D, bool R> struct GetReadType_;
 
-template<class D> struct GetReadType_<D, false, true>
+template<class D> struct GetReadType_<D, true>
 {
- using type = decltype(D().Read(MString(), 0));
+ using type = typename D::Data;
 };
 
-template<class D> struct GetReadType_<D, false, false>
+template<class D> struct GetReadType_<D, false>
 {
  using type = std::decay_t<D>;
 };
@@ -125,7 +119,7 @@ template<class Act, class S> consteval bool IsDisabled()
  }
 }
 
-template<class D> using ReadType = internal::GetReadType_<D, ReadPSupported<D>, ReadSupported<D>>::type;
+template<class D> using ReadType = internal::GetReadType_<D, ReadPSupported<D> || ReadSupported<D>>::type;
 
 template<class T>
 concept ReadIsGrid = requires {
diff --git a/sources/AVISOLOCAL.cpp b/sources/AVISOLOCAL.cpp
index 093ab3f..1847d58 100644
--- a/sources/AVISOLOCAL.cpp
+++ b/sources/AVISOLOCAL.cpp
@@ -86,9 +86,10 @@ std::pair<const BaseParameters*, MString> AVISOLOCALData::Parameters(michlib_int
  return {ppar.release(), ""};
 }
 
-AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
+bool AVISOLOCALData::Read(const MString& vname, std::map<MString, AVISOLOCALData::Data>& cache, const BaseParameters* ip, size_t i) const
 {
- if(!isOk()) return Data();
+ if(cache.contains(vname)) return true;
+ if(!isOk()) return false;
 
  auto    p = dynamic_cast<const struct Parameters*>(ip);
  NCFileA nc;
@@ -96,7 +97,7 @@ AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParame
  bool    isfloat = false, isint2 = false, isint4 = false;
 
  nc.Reset(datapath + "/uv-" + times[i].ToString() + ".nc");
- if(!nc) return Data();
+ if(!nc) return false;
  {
   auto head = nc.Header();
   for(const auto& v: head.Variables())
@@ -119,27 +120,36 @@ AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParame
   if(vname == "U" || vname == "U2")
   {
    bool square = vname == "U2";
-   auto u      = Read("u", ip, i);
-   auto v      = Read("v", ip, i);
-   if(!(u && v)) return Data();
-   if(square && u.Unit().Exist()) u.SetUnit("(" + u.Unit() + ")2");
-   for(size_t ind = 0; ind < u.N(); ind++)
+   if(!(Read("u", cache, ip, i) && Read("v", cache, ip, i))) return false;
+   cache[vname]  = cache.at("u").CopyGrid();
+   auto&       U = cache.at(vname);
+   const auto& u = cache.at("u");
+   const auto& v = cache.at("v");
+   if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
+   for(size_t ind = 0; ind < U.N(); ind++)
    {
     if(u.IsFill(ind) || v.IsFill(ind))
-     u.V(ind) = u.Fillval();
+     U.V(ind) = U.Fillval();
     else
-     u.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
+     U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
    }
-   return u;
+   return true;
   }
-  return Data();
+  return false;
  }
 
  // Direct read
- if(isint2) return ReadVarRaw<int2>(nc, name, p);
- if(isint4) return ReadVarRaw<int4>(nc, name, p);
- if(isfloat) return ReadVarRaw<float>(nc, name, p);
- return Data();
+ Data data;
+ if(isint2) data = ReadVarRaw<int2>(nc, name, p);
+ if(isint4) data = ReadVarRaw<int4>(nc, name, p);
+ if(isfloat) data = ReadVarRaw<float>(nc, name, p);
+ if(data)
+ {
+  cache[vname] = std::move(data);
+  return true;
+ }
+
+ return false;
 }
 
 template<class DataType> AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const NCFileA& nc, const MString& name, const struct AVISOLOCALData::Parameters* p) const
diff --git a/sources/AVISOLOCAL.h b/sources/AVISOLOCAL.h
index 1365d28..cd9f285 100644
--- a/sources/AVISOLOCAL.h
+++ b/sources/AVISOLOCAL.h
@@ -60,7 +60,7 @@ class AVISOLOCALData: public NCFuncs
   return NCFuncs::CheckVar(vname, [&nc = std::as_const(nc)](const MString& vn) { return HaveVar(nc, vn); });
  }
 
- Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
+ bool Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t i) const;
 
  template<class DataType> Data ReadVarRaw(const NCFileA& nc, const MString& name, const struct Parameters* p) const;
 };
diff --git a/sources/BINFILE.cpp b/sources/BINFILE.cpp
index 596d6e2..7a4f111 100644
--- a/sources/BINFILE.cpp
+++ b/sources/BINFILE.cpp
@@ -53,9 +53,10 @@ MString BINFILEData::Open(const CLArgs& args)
  return "";
 }
 
-BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
+bool BINFILEData::Read(const MString& vname, std::map<MString, BINFILEData::Data>& cache, size_t i) const
 {
- if(!isOk()) return Data();
+ if(cache.contains(vname)) return true;
+ if(!isOk()) return false;
 
  // Only rectangular grids are supported
  real xs = (data->Lon(data->Nx() - 1, data->Ny() - 1) - data->Lon(0, 0)) / (data->Nx() - 1);
@@ -66,18 +67,20 @@ BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
  if(vname == "U" || vname == "U2")
  {
   bool square = vname == "U2";
-  auto u      = Read("u", i);
-  auto v      = Read("v", i);
-  if(!(u && v)) return Data();
-  if(square && u.Unit().Exist()) u.SetUnit("(" + u.Unit() + ")2");
-  for(size_t ind = 0; ind < u.N(); ind++)
+  if(!(Read("u", cache, i) && Read("v", cache, i))) return false;
+  cache[vname]  = cache.at("u").CopyGrid();
+  auto&       U = cache.at(vname);
+  const auto& u = cache.at("u");
+  const auto& v = cache.at("v");
+  if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
+  for(size_t ind = 0; ind < U.N(); ind++)
   {
    if(u.IsFill(ind) || v.IsFill(ind))
-    u.V(ind) = u.Fillval();
+    U.V(ind) = U.Fillval();
    else
-    u.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
+    U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
   }
-  return u;
+  return true;
  }
 
  if(vname == "u" || vname == "v")
@@ -91,7 +94,13 @@ BINFILEData::Data BINFILEData::Read(const MString& vname, size_t i) const
     else
      out.V(ix, iy) = out.Fillval();
    }
+
+  if(out)
+  {
+   cache[vname] = std::move(out);
+   return true;
+  }
  }
 
- return out;
+ return false;
 }
diff --git a/sources/BINFILE.h b/sources/BINFILE.h
index 00018cc..ad9d7b9 100644
--- a/sources/BINFILE.h
+++ b/sources/BINFILE.h
@@ -44,5 +44,5 @@ class BINFILEData
 
  bool CheckVar(const MString& vname) const { return vname == "u" || vname == "v" || vname == "U" || vname == "U2"; }
 
- Data Read(const MString& vname, size_t i) const;
+ bool Read(const MString& vname, std::map<MString, Data>& cache, size_t i) const;
 };
diff --git a/sources/MODISBINLOCAL.cpp b/sources/MODISBINLOCAL.cpp
index 36e5151..bf5da18 100644
--- a/sources/MODISBINLOCAL.cpp
+++ b/sources/MODISBINLOCAL.cpp
@@ -83,8 +83,9 @@ std::pair<const BaseParameters*, MString> MODISBINLOCALData::Parameters(michlib_
  return {ppar.release(), ""};
 }
 
-MODISBINLOCALData::Data MODISBINLOCALData::Read(const MString& var, const BaseParameters* ip, size_t tind) const
+bool MODISBINLOCALData::Read(const MString& vname, std::map<MString, MODISBINLOCALData::Data>& cache, const BaseParameters* ip, size_t tind) const
 {
+ if(cache.contains(vname)) return true;
  auto p = pointer_cast<const struct Parameters*>(ip);
 
  Averager<Data> out;
@@ -94,7 +95,12 @@ MODISBINLOCALData::Data MODISBINLOCALData::Read(const MString& var, const BasePa
   Data dat = ReadFile(suf, p, tind);
   if(dat) out.Add(std::move(dat));
  }
- return std::move(out.Div());
+ if(out)
+ {
+  cache[vname] = std::move(out.Div());
+  return true;
+ }
+ return false;
 }
 
 MODISBINLOCALData::Data MODISBINLOCALData::ReadFile(const MString& suf, const struct MODISBINLOCALData::Parameters* p, size_t tind) const
diff --git a/sources/MODISBINLOCAL.h b/sources/MODISBINLOCAL.h
index b6348f8..663647f 100644
--- a/sources/MODISBINLOCAL.h
+++ b/sources/MODISBINLOCAL.h
@@ -40,7 +40,7 @@ class MODISBINLOCALData
 
  MString Info() const;
  MString Open(const CLArgs& args);
- Data    Read(const MString& var, const BaseParameters* ip, size_t tind) const;
+ bool    Read(const MString& vname, std::map<MString, Data>& cache, const BaseParameters* ip, size_t tind) const;
  Data    ReadFile(const MString& suf, const struct Parameters* p, size_t tind) const;
 
  std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args, const struct Region& reg) const;
diff --git a/src/layereddata.cpp b/src/layereddata.cpp
index 98fbc63..f15ab7b 100644
--- a/src/layereddata.cpp
+++ b/src/layereddata.cpp
@@ -178,10 +178,10 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
  return {ppar.release(), ""};
 }
 
-LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
+bool LayeredData::Read(const MString& vname, std::map<MString, LayeredData::Data>& cache, const BaseParameters* ip, size_t i) const
 {
- if(!isOk()) return Data();
- bool nodepth = false;
+ if(cache.contains(vname)) return true;
+ if(!isOk()) return false;
 
  auto p               = dynamic_cast<const struct Parameters*>(ip);
  auto [name, id, tid] = VarNameLoc(vname, times[i]);
@@ -190,86 +190,101 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters*
   // ptemp from temp and sal
   if(vname == "ptemp")
   {
-   auto temp = Read("temp", ip, i);
-   auto sal  = Read("sal", ip, i);
-   if(!(temp && sal)) return Data();
-   temp.SetUnit("degrees_C");
-   for(size_t ind = 0; ind < temp.N(); ind++)
+   if(!(Read("temp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
+   cache[vname]      = cache.at("temp").CopyGrid();
+   auto&       ptemp = cache.at(vname);
+   const auto& temp  = cache.at("temp");
+   const auto& sal   = cache.at("sal");
+   ptemp.SetUnit("degrees_C");
+   for(size_t ind = 0; ind < ptemp.N(); ind++)
    {
     if(temp.IsFill(ind) || sal.IsFill(ind))
-     temp.V(ind) = temp.Fillval();
+     ptemp.V(ind) = ptemp.Fillval();
     else
-     temp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind));
+     ptemp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), ptemp.Lon(ind), ptemp.Lat(ind));
    }
-   return temp;
+   return true;
   }
   // temp from ptemp and sal
   if(vname == "temp")
   {
-   auto temp = Read("ptemp", ip, i);
-   auto sal  = Read("sal", ip, i);
-   if(!(temp && sal)) return Data();
+   if(!(Read("ptemp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
+   cache[vname]      = cache.at("ptemp").CopyGrid();
+   auto&       temp  = cache.at(vname);
+   const auto& ptemp = cache.at("ptemp");
+   const auto& sal   = cache.at("sal");
    temp.SetUnit("degrees_C");
    for(size_t ind = 0; ind < temp.N(); ind++)
    {
-    if(temp.IsFill(ind) || sal.IsFill(ind))
+    if(ptemp.IsFill(ind) || sal.IsFill(ind))
      temp.V(ind) = temp.Fillval();
     else
-     temp.V(ind) = PTemp2Temp(temp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind));
+     temp.V(ind) = PTemp2Temp(ptemp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind));
    }
-   return temp;
+   return true;
   }
   // pdens from temp and sal
   if(vname == "pdens")
   {
    bool tempispot = HaveVar("ptemp");
-   auto temp      = Read(tempispot ? "ptemp" : "temp", ip, i);
-   auto sal       = Read("sal", ip, i);
-   if(!(temp && sal)) return Data();
-   temp.SetUnit("kg m-3");
-   for(size_t ind = 0; ind < temp.N(); ind++)
+   if(!(Read(tempispot ? "ptemp" : "temp", cache, ip, i) && Read("sal", cache, ip, i))) return false;
+   cache[vname]      = cache.at("sal").CopyGrid();
+   auto&       pdens = cache.at(vname);
+   const auto& temp  = cache.at(tempispot ? "ptemp" : "temp");
+   const auto& sal   = cache.at("sal");
+   pdens.SetUnit("kg m-3");
+   for(size_t ind = 0; ind < pdens.N(); ind++)
    {
     if(temp.IsFill(ind) || sal.IsFill(ind))
-     temp.V(ind) = temp.Fillval();
+     pdens.V(ind) = pdens.Fillval();
     else
-     temp.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind))
-                             : Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), temp.Lon(ind), temp.Lat(ind));
+     pdens.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), pdens.Lon(ind), pdens.Lat(ind))
+                              : Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), pdens.Lon(ind), pdens.Lat(ind));
    }
-   return temp;
+   return true;
   }
   // U and U2 from u and v
   if(vname == "U" || vname == "U2")
   {
    bool square = vname == "U2";
-   auto u      = Read("u", ip, i);
-   auto v      = Read("v", ip, i);
-   if(!(u && v)) return Data();
-   if(square && u.Unit().Exist()) u.SetUnit("(" + u.Unit() + ")2");
-   for(size_t ind = 0; ind < u.N(); ind++)
+   if(!(Read("u", cache, ip, i) && Read("v", cache, ip, i))) return false;
+   cache[vname]  = cache.at("u").CopyGrid();
+   auto&       U = cache.at(vname);
+   const auto& u = cache.at("u");
+   const auto& v = cache.at("v");
+   if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
+   for(size_t ind = 0; ind < U.N(); ind++)
    {
     if(u.IsFill(ind) || v.IsFill(ind))
-     u.V(ind) = u.Fillval();
+     U.V(ind) = U.Fillval();
     else
-     u.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
+     U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
    }
-   return u;
+   return true;
   }
-  return Data();
+  return false;
  }
 
  // Direct read
+ bool nodepth = false;
+ Data data;
  auto head = nc[id]->Header();
  for(const auto& v: head.Variables())
   if(v.Name() == name)
   {
    if(v.Dimensions().size() == 3) nodepth = true;
-   if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
-   if(v.Type().Id() == NC_INT) return ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
-   if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
-   if(v.Type().Id() == NC_DOUBLE) return ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
+   if(v.Type().Id() == NC_SHORT) data = ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
+   if(v.Type().Id() == NC_INT) data = ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
+   if(v.Type().Id() == NC_FLOAT) data = ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
+   if(v.Type().Id() == NC_DOUBLE) data = ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
+   if(data)
+   {
+    cache[vname] = std::move(data);
+    return true;
+   }
   }
 
- return Data();
+ return false;
 }
 
 template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const