Added support for conversion from in-situ temperature to potential temperature and vice versa, as well as potential density calculations.

include/basedata.h

@@ -42,6 +42,8 @@ class BaseData
 
  size_t N() const { return data.size(); }
 
+ bool IsFill(size_t i) const { return V(i) == Fillval(); }
+
  static real Fillval() { return fillval; }
 
  explicit operator bool() const { return N() != 0; }

include/gsw.h

@@ -0,0 +1,25 @@
+#pragma once
+#include "gswteos-10.h"
+
+inline double Temp2PTemp(double temp, double psal, double depth, double lon, double lat) { return gsw_pt0_from_t(gsw_sa_from_sp(psal, depth, lon, lat), temp, depth); }
+
+inline double PTemp2Temp(double ptemp, double psal, double depth, double lon, double lat)
+{
+ double sa = gsw_sa_from_sp(psal, depth, lon, lat);
+ return gsw_t_from_ct(sa, gsw_ct_from_pt(sa, ptemp), depth);
+}
+
+inline double PTemp2PDens(double ptemp, double psal, double depth, double lon, double lat) { return gsw_rho_t_exact(gsw_sa_from_sp(psal, depth, lon, lat), ptemp, 0.0) - 1000.0; }
+
+inline double Temp2PDens(double temp, double psal, double depth, double lon, double lat)
+{
+ return gsw_pot_rho_t_exact(gsw_sa_from_sp(psal, depth, lon, lat), temp, depth, 0.0) - 1000.0;
+}
+
+inline double PTemp2SSpeed(double ptemp, double psal, double depth, double lon, double lat)
+{
+ double sa = gsw_sa_from_sp(psal, depth, lon, lat);
+ return gsw_sound_speed_t_exact(sa, gsw_t_from_ct(sa, gsw_ct_from_pt(sa, ptemp), depth), depth);
+}
+
+inline double Temp2SSpeed(double temp, double psal, double depth, double lon, double lat) { return gsw_sound_speed_t_exact(gsw_sa_from_sp(psal, depth, lon, lat), temp, depth); }

include/layereddata.h

@@ -1,6 +1,7 @@
 #pragma once
 #include "DataAdapters/ncfilealt.h"
 #include "ParseArgs.h"
+#include "gsw.h"
 #include "mdatetime.h"
 #include "simple2ddata.h"
 #include <algorithm>
@@ -220,7 +221,7 @@ class LayeredData
 
  static MString StName2Name(const MString& stname)
  {
-  if(stname == "sea_water_potential_temperature") return "temp";
+  if(stname == "sea_water_potential_temperature") return "ptemp";
   if(stname == "sea_water_temperature") return "temp";
   if(stname == "sea_water_salinity") return "sal";
   if(stname == "ocean_mixed_layer_thickness_defined_by_sigma_theta") return "mld";

src/layereddata.cpp

@@ -18,6 +18,10 @@ MString LayeredData::Info() const
    if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
   }
  }
+ if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
+ if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
+ if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");
+
  MString svars;
  {
   bool first = true;
@@ -136,7 +140,14 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
 
  if(!args.contains("var")) return {nullptr, "Variable not specified"};
  ppar->varname = args.at("var");
- if(!HaveVar(ppar->varname)) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"};
+ if(!HaveVar(ppar->varname))
+ {
+  bool varexist = false;
+  if(ppar->varname == "temp" && HaveVar("ptemp") && HaveVar("sal")) varexist = true;
+  if(ppar->varname == "ptemp" && HaveVar("temp") && HaveVar("sal")) varexist = true;
+  if(ppar->varname == "pdens" && (HaveVar("ptemp") || HaveVar("temp")) && HaveVar("sal")) varexist = true;
+  if(!varexist) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"};
+ }
  if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger<size_t>();
  if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"};
  real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer);
@@ -192,7 +203,87 @@ LayeredData::Data LayeredData::Read(const BaseParameters* ip, size_t i) const
 
  auto p               = dynamic_cast<const struct Parameters*>(ip);
  auto [name, id, tid] = VarNameLoc(p->varname, times[i]);
- auto head            = nc[id]->Header();
+ if(!name.Exist()) // Conversion read
+ {
+  // ptemp from temp and sal
+  if(p->varname == "ptemp")
+  {
+   struct Parameters params = *p;
+   params.varname           = "temp";
+   auto temp                = Read(&params, i);
+   params.varname           = "sal";
+   auto sal                 = Read(&params, i);
+   if(!(temp && sal)) return Data();
+   auto out = temp;
+   for(size_t ind = 0; ind < out.N(); ind++)
+   {
+    if(out.Lon(ind) != temp.Lon(ind) || out.Lon(ind) != sal.Lon(ind) || out.Lat(ind) != temp.Lat(ind) || out.Lat(ind) != sal.Lat(ind))
+    {
+     michlib::errmessage("Internal error");
+     exit(1);
+    }
+    if(temp.IsFill(ind) || sal.IsFill(ind))
+     out.V(ind) = out.Fillval();
+    else
+     out.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
+   }
+   return out;
+  }
+  // temp from ptemp and sal
+  if(p->varname == "temp")
+  {
+   struct Parameters params = *p;
+   params.varname           = "ptemp";
+   auto temp                = Read(&params, i);
+   params.varname           = "sal";
+   auto sal                 = Read(&params, i);
+   if(!(temp && sal)) return Data();
+   auto out = temp;
+   for(size_t ind = 0; ind < out.N(); ind++)
+   {
+    if(out.Lon(ind) != temp.Lon(ind) || out.Lon(ind) != sal.Lon(ind) || out.Lat(ind) != temp.Lat(ind) || out.Lat(ind) != sal.Lat(ind))
+    {
+     michlib::errmessage("Internal error");
+     exit(1);
+    }
+    if(temp.IsFill(ind) || sal.IsFill(ind))
+     out.V(ind) = out.Fillval();
+    else
+     out.V(ind) = PTemp2Temp(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
+   }
+   return out;
+  }
+  // pdens from temp and sal
+  if(p->varname == "pdens")
+  {
+   struct Parameters params    = *p;
+   bool              tempispot = HaveVar("ptemp");
+   params.varname              = tempispot ? "ptemp" : "temp";
+   auto temp                   = Read(&params, i);
+   params.varname              = "sal";
+   auto sal                    = Read(&params, i);
+   if(!(temp && sal)) return Data();
+   auto out = temp;
+   for(size_t ind = 0; ind < out.N(); ind++)
+   {
+    if(out.Lon(ind) != temp.Lon(ind) || out.Lon(ind) != sal.Lon(ind) || out.Lat(ind) != temp.Lat(ind) || out.Lat(ind) != sal.Lat(ind))
+    {
+     michlib::errmessage("Internal error");
+     exit(1);
+    }
+    if(temp.IsFill(ind) || sal.IsFill(ind))
+     out.V(ind) = out.Fillval();
+    else
+     out.V(ind) = tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind))
+                            : Temp2PDens(temp.V(ind), sal.V(ind), Depth(p->layer), out.Lon(ind), out.Lat(ind));
+   }
+   return out;
+  }
+  return Data();
+ }
+
+ // Direct read
+ auto head = nc[id]->Header();
  for(const auto& v: head.Variables())
   if(v.Name() == name)
   {

src/odm.cpp

@@ -13,7 +13,7 @@ inline void Usage(const MString& arg0)
  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("  var. Required. May be: temp, sal, chl, mld, ssh or w.");
+ 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");
  message("  out. Output file. Default: out.bin");