diff --git a/include/basedata.h b/include/basedata.h
index 3bfac64..2a4e445 100644
--- a/include/basedata.h
+++ b/include/basedata.h
@@ -49,9 +49,34 @@ class BaseData
  explicit operator bool() const { return N() != 0; }
 };
 
+template<class T>
+concept HaveU = requires(T t) {
+ {
+  t.template U(0)
+ } -> std::convertible_to<real>;
+};
+
+template<class T>
+concept HaveV = requires(T t) {
+ {
+  t.template V(0)
+ } -> std::convertible_to<real>;
+};
+
 template<class Data> class Averager: public Data
 {
- std::vector<size_t> count;
+ std::vector<size_t>   count;
+ static constexpr bool haveu = HaveU<Data>;
+ static constexpr bool havev = HaveV<Data>;
+
+ static bool DataOk(const Data* d, size_t i)
+ {
+  bool uok = true;
+  bool vok = true;
+  if constexpr(haveu) uok = d->U(i) != Data::Fillval();
+  if constexpr(havev) uok = d->V(i) != Data::Fillval();
+  return uok && vok;
+ }
 
  public:
  void Add(const Data& d)
@@ -61,23 +86,25 @@ template<class Data> class Averager: public Data
   {
    *static_cast<Data*>(this) = d;
    count.resize(Data::N());
-   for(size_t i = 0; i < Data::N(); i++) count[i] = (Data::V(i) == Data::Fillval()) ? 0 : 1;
+   for(size_t i = 0; i < Data::N(); i++) count[i] = DataOk(this, i) ? 1 : 0;
    return;
   }
 
   if(Data::N() != d.N()) return;
   for(size_t i = 0; i < Data::N(); i++)
-   if(d(i) != Data::Fillval())
+   if(DataOk(&d, i))
    {
-    if(Data::V(i) == Data::Fillval())
+    if(DataOk(this, i))
     {
-     Data::V(i) = d(i);
-     count[i]   = 1;
+     if constexpr(haveu) Data::U(i) += d.U(i);
+     if constexpr(havev) Data::V(i) += d.V(i);
+     count[i]++;
     }
     else
     {
-     Data::V(i) += d(i);
-     count[i]++;
+     if constexpr(haveu) Data::U(i) = d.U(i);
+     if constexpr(havev) Data::V(i) = d.V(i);
+     count[i] = 1;
     }
    }
  }
@@ -86,7 +113,16 @@ template<class Data> class Averager: public Data
  {
   if(!*this) return *this;
   for(size_t i = 0; i < Data::N(); i++)
-   if(count[i] != 0) Data::V(i) /= count[i];
+  {
+   if constexpr(haveu)
+   {
+    if(count[i] != 0) Data::U(i) /= count[i];
+   }
+   if constexpr(havev)
+   {
+    if(count[i] != 0) Data::V(i) /= count[i];
+   }
+  }
   return *this;
  }
 };
diff --git a/include/layereddata.h b/include/layereddata.h
index 0943721..99458ae 100644
--- a/include/layereddata.h
+++ b/include/layereddata.h
@@ -4,6 +4,7 @@
 #include "gsw.h"
 #include "mdatetime.h"
 #include "simple2ddata.h"
+#include "uvdata.h"
 #include <algorithm>
 #include <memory>
 #include <set>
@@ -148,6 +149,8 @@ class LayeredData
 
  Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
 
+ UVData ReadUV(const BaseParameters* ip, size_t i) const;
+
  bool isOk() const { return nc.size() > 0; }
 
  explicit operator bool() const { return nc.size() > 0; }
@@ -242,6 +245,8 @@ class LayeredData
   if(stname == "sea_surface_elevation") return "ssh";
   if(stname == "eastward_sea_water_velocity") return "u";
   if(stname == "northward_sea_water_velocity") return "v";
+  if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u";
+  if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v";
   if(stname == "upward_sea_water_velocity") return "w";
   return "";
  }
diff --git a/include/odm.h b/include/odm.h
index 21a95e1..f7c7a41 100644
--- a/include/odm.h
+++ b/include/odm.h
@@ -40,12 +40,9 @@ concept ReadSupported = requires(T t, const MString& vname, size_t i) {
 };
 
 template<class T>
-concept ActionTscSupported = requires(T t, michlib_internal::ParameterListEx& pars, const CLArgs& args, const BaseParameters* ip, size_t i) {
+concept ReadUVPSupported = requires(T t, const BaseParameters* ip, size_t i) {
  {
-  t.template Parameters(pars, args).first
- } -> std::convertible_to<const BaseParameters*>;
- {
-  t.template Read(ip, i)(0)
+  t.template ReadUV(ip, i).U(0, 0) + t.template ReadUV(ip, i).V(0, 0)
  } -> std::convertible_to<real>;
 };
 
@@ -163,4 +160,5 @@ class Data: public DataVariants
 
  MString ActionInfo(const CLArgs& args);
  MString ActionTsc(const CLArgs& args);
+ MString ActionUV(const CLArgs& args);
 };
diff --git a/include/uvdata.h b/include/uvdata.h
new file mode 100644
index 0000000..4b16e39
--- /dev/null
+++ b/include/uvdata.h
@@ -0,0 +1,225 @@
+#pragma once
+#include "geohelpers.h"
+#include <vector>
+
+using michlib::M_PI;
+using michlib::real;
+using michlib::Sqrt;
+
+enum class Metric
+{
+ EUCLID,
+ SPHERIC
+};
+
+class UVData
+{
+ static constexpr real fillval = 1.0e10;
+ real                  x0 = 0.0, y0 = 0.0;
+ size_t                nx = 0, ny = 0;
+ real                  xstep = 0.0, ystep = 0.0;
+
+ std::vector<real> u, v;
+ Metric            metric;
+
+ real D(real lon1, real lat1, real lon2, real lat2) const
+ {
+  switch(metric)
+  {
+  case(Metric::EUCLID): return michlib::Hypot(lon2 - lon1, lat2 - lat1);
+  // Spheric distance in km
+  case(Metric::SPHERIC): return michlib::GCD(lon1 * M_PI / 180.0, lat1 * M_PI / 180.0, lon2 * M_PI / 180.0, lat2 * M_PI / 180.0) * 6371.0;
+  }
+  return 0.0;
+ }
+
+ // For derivatives
+ real Ud(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : U(ix, iy); }
+ real Vd(size_t ix, size_t iy) const { return IsCoast(ix, iy) ? 0.0 : V(ix, iy); }
+
+ real dUdX(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return fillval;
+  if(ix == 0) return (Ud(1, iy) - Ud(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
+  if(ix == nx - 1) return (Ud(nx - 1, iy) - Ud(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
+  return 0.5 * ((Ud(ix + 1, iy) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
+                (Ud(ix, iy) - Ud(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dVdX(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return fillval;
+  if(ix == 0) return (Vd(1, iy) - Vd(0, iy)) / D(Lon(0, iy), Lat(0, iy), Lon(1, iy), Lat(1, iy));
+  if(ix == nx - 1) return (Vd(nx - 1, iy) - Vd(nx - 2, iy)) / D(Lon(nx - 2, iy), Lat(nx - 2, iy), Lon(nx - 1, iy), Lat(nx - 1, iy));
+  return 0.5 * ((Vd(ix + 1, iy) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix + 1, iy), Lat(ix + 1, iy)) +
+                (Vd(ix, iy) - Vd(ix - 1, iy)) / D(Lon(ix - 1, iy), Lat(ix - 1, iy), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dUdY(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return fillval;
+  if(iy == 0) return (Ud(ix, 1) - Ud(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
+  if(iy == ny - 1) return (Ud(ix, ny - 1) - Ud(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
+  return 0.5 * ((Ud(ix, iy + 1) - Ud(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
+                (Ud(ix, iy) - Ud(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
+ }
+ real dVdY(size_t ix, size_t iy) const
+ {
+  if(IsCoast(ix, iy)) return fillval;
+  if(iy == 0) return (Vd(ix, 1) - Vd(ix, 0)) / D(Lon(ix, 0), Lat(ix, 0), Lon(ix, 1), Lat(ix, 1));
+  if(iy == ny - 1) return (Vd(ix, ny - 1) - Vd(ix, ny - 2)) / D(Lon(ix, ny - 2), Lat(ix, ny - 2), Lon(ix, ny - 1), Lat(ix, ny - 1));
+  return 0.5 * ((Vd(ix, iy + 1) - Vd(ix, iy)) / D(Lon(ix, iy), Lat(ix, iy), Lon(ix, iy + 1), Lat(ix, iy + 1)) +
+                (Vd(ix, iy) - Vd(ix, iy - 1)) / D(Lon(ix, iy - 1), Lat(ix, iy - 1), Lon(ix, iy), Lat(ix, iy)));
+ }
+
+ // For stationary points
+ real LonR(real x, [[maybe_unused]] real y) const { return x0 + x * xstep; }
+ real LatR([[maybe_unused]] real x, real y) const { return y0 + y * ystep; }
+
+ public:
+ enum STPOINT
+ {
+  SADDLE = 0,
+  SACICFOCUS,
+  SKNOT,
+  UACICFOCUS,
+  UKNOT,
+  SCICFOCUS,
+  UCICFOCUS,
+  NOPOINT
+ };
+
+ struct StPoint
+ {
+  real    x = 0.0, y = 0.0;
+  STPOINT type = NOPOINT;
+ };
+
+ UVData() = default;
+ UVData(size_t nx_, size_t ny_, real x0_, real y0_, real xs_, real ys_, Metric m_ = Metric::SPHERIC):
+     x0(x0_), y0(y0_), nx(nx_), ny(ny_), xstep(xs_), ystep(ys_), u(nx_ * ny_), v(nx_ * ny_), metric(m_)
+ {
+ }
+
+ const real& U(size_t i) const { return u[i]; }
+ const real& V(size_t i) const { return v[i]; }
+ const real& U(size_t ix, size_t iy) const { return U(iy * nx + ix); }
+ const real& V(size_t ix, size_t iy) const { return V(iy * nx + ix); }
+
+ real& U(size_t i) { return u[i]; }
+ real& V(size_t i) { return v[i]; }
+ real& U(size_t ix, size_t iy) { return U(iy * nx + ix); }
+ real& V(size_t ix, size_t iy) { return V(iy * nx + ix); }
+
+ size_t N() const { return u.size(); }
+ size_t Nx() const { return nx; }
+ size_t Ny() const { return ny; }
+
+ real Lon(size_t ix, [[maybe_unused]] size_t iy) const { return x0 + ix * xstep; }
+ real Lat([[maybe_unused]] size_t ix, size_t iy) const { return y0 + iy * ystep; }
+
+ real Lon(size_t i) const { return Lon(i % nx, i / nx); }
+ real Lat(size_t i) const { return Lat(i % nx, i / nx); }
+
+ explicit operator bool() const { return N() != 0; }
+
+ bool IsCoast(size_t i) const { return U(i) == fillval || V(i) == fillval; }
+ bool IsCoast(size_t ix, size_t iy) const { return U(ix, iy) == fillval || V(ix, iy) == fillval; }
+
+ static real Fillval() { return fillval; }
+
+ real Div(size_t i) const { return Div(i % nx, i / nx); }
+ real Rot(size_t i) const { return Rot(i % nx, i / nx); }
+ real Div(size_t ix, size_t iy) const
+ {
+  if(!*this) return 0.0;
+  real ux = dUdX(ix, iy);
+  real vy = dVdY(ix, iy);
+  return (ux == fillval || vy == fillval) ? fillval : (ux + vy);
+ }
+ real Rot(size_t ix, size_t iy) const
+ {
+  if(!*this) return 0.0;
+  real vx = dVdX(ix, iy);
+  real uy = dUdY(ix, iy);
+  return (vx == fillval || uy == fillval) ? fillval : (vx - uy);
+ }
+ // Okubo-Weiss parameter
+ real OW(size_t i) const { return OW(i % nx, i / nx); }
+ real OW(size_t ix, size_t iy) const
+ {
+  if(!*this) return 0.0;
+  real ux = dUdX(ix, iy);
+  real uy = dUdY(ix, iy);
+  real vx = dVdX(ix, iy);
+  real vy = dVdY(ix, iy);
+
+  if(ux == fillval || uy == fillval || vx == fillval || vy == fillval) return fillval;
+
+  real sn = ux - vy;
+  real ss = vx + uy;
+  real w  = vx - uy;
+  return sn * sn + ss * ss - w * w;
+ }
+
+ auto StablePoints(size_t ix, size_t iy) const
+ {
+  std::vector<struct StPoint> points;
+  if(!*this) return points;
+  if(ix >= Nx() - 1 || iy >= Ny() - 1) return points;
+  if(IsCoast(ix, iy) || IsCoast(ix + 1, iy) || IsCoast(ix, iy + 1) || IsCoast(ix + 1, iy + 1)) return points;
+
+  real au = U(ix, iy) - U(ix, iy + 1) - U(ix + 1, iy) + U(ix + 1, iy + 1);
+  real bu = U(ix, iy + 1) - U(ix, iy);
+  real cu = U(ix + 1, iy) - U(ix, iy);
+  real du = U(ix, iy);
+
+  real av = V(ix, iy) - V(ix, iy + 1) - V(ix + 1, iy) + V(ix + 1, iy + 1);
+  real bv = V(ix, iy + 1) - V(ix, iy);
+  real cv = V(ix + 1, iy) - V(ix, iy);
+  real dv = V(ix, iy);
+
+  real D   = au * bv - av * bu;
+  real a   = (au * av * cu - au * au * cv) / D;
+  real b   = (au * bv * cu - au * bu * cv + au * av * du - au * au * dv) / D;
+  real c   = (au * bv * du - au * bu * dv) / D;
+  real des = b * b - 4.0 * a * c;
+  if(des < 0) return points;
+  des = Sqrt(des);
+
+  real p  = av * cu - au * cv;
+  real q  = av * bu - au * bv;
+  real r  = av * du - au * dv;
+  real x1 = 0.5 * (-b + des) / a;
+  real x2 = 0.5 * (-b - des) / a;
+  real y1 = -(p * x1 + r) / q;
+  real y2 = -(p * x2 + r) / q;
+
+  auto PointType = [au = au, av = av, bu = bu, bv = bv, cu = cu, cv = cv](real x, real y) -> auto
+  {
+   real a   = 1.0;
+   real b   = -y * au - x * av - bv - cu;
+   real c   = y * (au * bv - av * bu) + x * (av * cu - au * cv) + bv * cu - bu * cv;
+   real des = b * b - 4.0 * a * c;
+   if(des > 0.0)
+   {
+    if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0 && c > 0.0)
+     return UKNOT;
+    else if(std::max(-b + Sqrt(des), -b - Sqrt(des)) > 0.0)
+     return SADDLE;
+    else
+     return SKNOT;
+   }
+   else
+   {
+    bool acyclcrit = (av * y + cv > 0.0);
+    if(b < 0.0)
+     return acyclcrit ? SACICFOCUS : SCICFOCUS;
+    else
+     return acyclcrit ? UACICFOCUS : UCICFOCUS;
+   }
+  };
+
+  if(x1 >= 0.0 && x1 < 1.0 && y1 >= 0.0 && y1 < 1.0) points.emplace_back(LonR(x1 + ix, y1 + iy), LatR(x1 + ix, y1 + iy), PointType(x1, y1));
+  if(x2 >= 0.0 && x2 < 1.0 && y2 >= 0.0 && y2 < 1.0) points.emplace_back(LonR(x2 + ix, y2 + iy), LatR(x2 + ix, y2 + iy), PointType(x2, y2));
+  return points;
+ }
+};
diff --git a/src/actiontsc.cpp b/src/actiontsc.cpp
index afad0b6..14a7abc 100644
--- a/src/actiontsc.cpp
+++ b/src/actiontsc.cpp
@@ -96,6 +96,9 @@ MString Data::ActionTsc(const CLArgs& args)
  if(!name.Exist()) name = "out.bin";
 
  fw.Create(name, 3);
+ fw.SetColumnName(1, "Longitude");
+ fw.SetColumnName(2, "Latitude");
+ fw.SetColumnName(3, vname);
  fw.SetParameters(pars);
  for(size_t i = 0; i < data.N(); i++)
  {
diff --git a/src/actionuv.cpp b/src/actionuv.cpp
new file mode 100644
index 0000000..8c3fb0f
--- /dev/null
+++ b/src/actionuv.cpp
@@ -0,0 +1,169 @@
+#define MICHLIB_NOSOURCE
+#include "GPL.h"
+#include "odm.h"
+
+MString Data::ActionUV(const CLArgs& args)
+{
+ if(args.contains("time") && (args.contains("timeb") || args.contains("timee"))) return "Time must be set via time parameter or timeb and timee parameter but not via both";
+ if(!(args.contains("time") || (args.contains("timeb") && args.contains("timee")))) return "Time must be set via time parameter or timeb and timee parameter";
+
+ michlib_internal::ParameterListEx pars;
+ pars.UsePrefix("");
+
+ pars.SetParameter("source", args.at("source"));
+ TIndex tindexes;
+
+ if(args.contains("time"))
+ {
+  tindexes = GetTIndexes(args.at("time")); // Regular expression or single date
+  if(tindexes.size() == 0) return "There are no times matching the regular expression";
+  if(tindexes.size() == 1)
+   pars.SetParameter("time", Time(tindexes[0]).ToString());
+  else
+   pars.SetParameter("timeregex", args.at("time"));
+ }
+ else
+ {
+  MDateTime tb(args.at("timeb")), te(args.at("timee"));
+  tindexes = GetTIndexes(tb, te);
+  if(tindexes.size() == 0) return "There are no times between " + tb.ToString() + " and " + te.ToString();
+  pars.SetParameter("timeb", tb.ToString());
+  pars.SetParameter("timee", te.ToString());
+ }
+
+ std::unique_ptr<const BaseParameters> sourcepars;
+ {
+  auto ppar = std::visit(
+      [&pars = pars, &args = args](const auto& arg) -> std::pair<const BaseParameters*, MString>
+      {
+       using T = std::decay_t<decltype(arg)>;
+       if constexpr(ParametersSupported<T>)
+        return arg.Parameters(pars, args);
+       else
+        return {nullptr, ""};
+      },
+      *this);
+  if(ppar.second.Exist()) return ppar.second;
+  sourcepars.reset(ppar.first);
+ }
+
+ auto   p = sourcepars.get();
+ size_t ind;
+ auto   read = [p = p, &ind = ind](const auto& arg) -> UVData
+ {
+  using T = std::decay_t<decltype(arg)>;
+  if constexpr(ReadUVPSupported<T>)
+   return arg.ReadUV(p, ind);
+  else
+   return UVData();
+ };
+
+ UVData data;
+ if(tindexes.size() == 1)
+ {
+  ind = tindexes[0];
+  michlib::message("Time: " + Time(ind).ToTString());
+  data = std::visit(read, *this);
+ }
+ else
+ {
+  Averager<UVData> out;
+  bool             ok = true;
+  for(size_t i = 0; i < tindexes.size(); i++)
+  {
+   if(!ok) break;
+   ind = tindexes[i];
+   michlib::message("Time: " + Time(ind).ToTString());
+   auto dat = std::visit(read, *this);
+   if(dat)
+    out.Add(dat);
+   else
+    ok = false;
+  }
+  if(ok) data = out.Div();
+ }
+ if(!data) return "Can't read data";
+
+ // Main file
+ MString name = args.contains("out") ? args.at("out") : "";
+ if(name.Exist())
+ {
+  BFileW fw;
+  fw.Create(name, 7);
+  fw.SetColumnName(1, "Longitude");
+  fw.SetColumnName(2, "Latitude");
+  fw.SetColumnName(3, "u, cm/s");
+  fw.SetColumnName(4, "v, cm/s");
+  fw.SetColumnName(5, "Div, (cm/s)/km");
+  fw.SetColumnName(6, "Rot, (cm/s)/km");
+  fw.SetColumnName(7, "Okubo-Weiss parameter, (cm^2/s^2)/km^2");
+  fw.SetParameters(pars);
+  for(size_t i = 0; i < data.N(); i++)
+  {
+   fw.Write(data.Lon(i));
+   fw.Write(data.Lat(i));
+   fw.Write(data.U(i) == data.Fillval() ? NAN : data.U(i));
+   fw.Write(data.V(i) == data.Fillval() ? NAN : data.V(i));
+   fw.Write(data.Div(i) == data.Fillval() ? NAN : data.Div(i));
+   fw.Write(data.Rot(i) == data.Fillval() ? NAN : data.Rot(i));
+   fw.Write(data.OW(i) == data.Fillval() ? NAN : data.OW(i));
+  }
+  fw.Finalize();
+  fw.Close();
+ }
+
+ // Filtered vectors file
+ name = args.contains("velout") ? args.at("velout") : "";
+ if(name.Exist())
+ {
+  size_t shiftx = args.contains("shiftx") ? args.at("shiftx").ToInteger<size_t>() : 0;
+  size_t shifty = args.contains("shifty") ? args.at("shifty").ToInteger<size_t>() : 0;
+  size_t skipx  = args.contains("skipx") ? args.at("skipx").ToInteger<size_t>() : 1;
+  size_t skipy  = args.contains("skipy") ? args.at("skipy").ToInteger<size_t>() : 1;
+
+  BFileW vel;
+  vel.Create(name, 4);
+  vel.SetColumnName(1, "Longitude");
+  vel.SetColumnName(2, "Latitude");
+  vel.SetColumnName(3, "u, cm/s");
+  vel.SetColumnName(4, "v, cm/s");
+
+  for(size_t ix = shiftx; ix < data.Nx(); ix += skipx)
+   for(size_t iy = shifty; iy < data.Ny(); iy += skipy)
+   {
+    vel.Write(data.Lon(ix, iy));
+    vel.Write(data.Lat(ix, iy));
+    vel.Write(data.U(ix, iy) == data.Fillval() ? NAN : data.U(ix, iy));
+    vel.Write(data.V(ix, iy) == data.Fillval() ? NAN : data.V(ix, iy));
+   }
+  vel.Finalize();
+  vel.Close();
+ }
+
+ // Stationary points
+ name = args.contains("stpout") ? args.at("stpout") : "";
+ if(name.Exist())
+ {
+  BFileW stp;
+  stp.Create(name, 3);
+  stp.SetColumnName(1, "Longitude");
+  stp.SetColumnName(2, "Latitude");
+  stp.SetColumnName(3, "Stability (0 - saddle, 1 - st. anticicl. focus, 2 - st. knot, 3 - unst. anticicl. focus, 4 - unst. knot, 5 - st. cicl. focus, 6 - unst. cicl. focus)");
+
+  for(size_t ix = 0; ix < data.Nx() - 1; ix++)
+   for(size_t iy = 0; iy < data.Ny() - 1; iy++)
+   {
+    auto points = data.StablePoints(ix, iy);
+    for(size_t i = 0; i < points.size(); i++)
+    {
+     stp.Write(points[i].x);
+     stp.Write(points[i].y);
+     stp.Write(michlib::int_cast<size_t>(points[i].type));
+    }
+   }
+  stp.Finalize();
+  stp.Close();
+ }
+
+ return "";
+}
diff --git a/src/layereddata.cpp b/src/layereddata.cpp
index 559a356..02d0e91 100644
--- a/src/layereddata.cpp
+++ b/src/layereddata.cpp
@@ -280,6 +280,28 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters*
  return Data();
 }
 
+UVData LayeredData::ReadUV(const BaseParameters* ip, size_t i) const
+{
+ if(!isOk()) return UVData();
+
+ auto u = Read("u", ip, i);
+ auto v = Read("v", ip, i);
+ if(!(u && v)) return UVData();
+
+ UVData out{u.Nx(), u.Ny(), u.Lon(0, 0), u.Lat(0, 0), lonstep, latstep};
+ for(size_t i = 0; i < out.N(); i++)
+ {
+  if(u(i) == Data::Fillval() || v(i) == Data::Fillval())
+   out.U(i) = out.V(i) = UVData::Fillval();
+  else
+  {
+   out.U(i) = u(i);
+   out.V(i) = v(i);
+  }
+ }
+ return out;
+}
+
 template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::Parameters* p) const
 {
  real     unitmul = 1.0;
diff --git a/src/odm.cpp b/src/odm.cpp
index 1dd7105..62dbfa7 100644
--- a/src/odm.cpp
+++ b/src/odm.cpp
@@ -4,7 +4,7 @@ inline void Usage(const MString& arg0)
 {
  message(arg0 + " (<key>=<value>...)");
  message("Keys are:");
- message(" action. What the program should do. May be: info, tsc. Default: tsc.");
+ message(" action. What the program should do. May be: info, tsc, uv. Default: tsc.");
  message(" Keys for action=info. Print some information about dataset.");
  message("  source. Required. May be: NEMO, HYCOM");
  message("  Keys for source=NEMO");
@@ -12,7 +12,7 @@ inline void Usage(const MString& arg0)
  message("  Keys for source=HYCOM");
  message("   dataset. Can be Forecast, Hindcast or Reanalysis. Default: Forecast");
  message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset.");
- message("  source. Required. May be: NEMO");
+ message("  source. Required. May be: NEMO, HYCOM");
  message("  var. Required. May be: temp, ptemp, pdens, sal, chl, mld, ssh or w.");
  message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
  message("  timeb, timee. Time interval. If present, time must be absent");
@@ -25,6 +25,17 @@ inline void Usage(const MString& arg0)
  message("   dataset. Can be Forecast, Hindcast or Reanalysis. Default: Forecast");
  message("   layer and/or depth. Layer or depth of HYCOM dataset. If depth is specified, layer is ignored. Both ignored, if var=mld or var=ssh. Default: layer=0");
  message("   lonb, lone, latb, late. Required. Region of interest");
+ message(" Keys for action=uv. Get velocity field and its derivatives.");
+ message("  source. Required. May be: NEMO, HYCOM");
+ message("  time. Time moment or regular expression. If present, timeb and timee must be absent");
+ message("  timeb, timee. Time interval. If present, time must be absent");
+ message("  out. Output file for components of velocity field, divergency, rotor and Okubo-Weiss parameter. If absent, this data not calculated.");
+ message("  velout. Output file for sparse velocity field. If present, parameters shiftx, shifty, skipx, skipy controls the sparsing.");
+ message("  shiftx, shifty. Shift of the sparsed grid. Used only if velout parameter is present. Default: 0");
+ message("  skipx, skipy. How many poinst skipped in the sparsed grid. Used only if velout parameter is present. Default: 1, output each point");
+ message("  stpout. Output file for stationary points. If absent, this data not calculated.");
+ message("  Keys for source=NEMO. See section action=tsc.");
+ message("  Keys for source=HYCOM. See section action=tsc.");
 }
 
 int main(int argc, char** argv)
@@ -53,6 +64,11 @@ int main(int argc, char** argv)
   ret = data.ActionInfo(args);
  else if(action == "tsc")
   ret = data.ActionTsc(args);
+ else if(action == "uv")
+ {
+  args["var"] = "U";
+  ret         = data.ActionUV(args);
+ }
  else
  {
   errmessage("Unknown action " + action);