odm/include/ncfuncs.h

#pragma once
#include "DataAdapters/ncfilealt.h"
#include "basedata.h"
#include "mdatetime.h"
#include "ncsimple.h"
#include "zarr.h"
#include <map>
#include <set>
#include <tuple>

using michlib::MDateTime;
using michlib::MString;
using michlib::NCFileA;

class NCFuncs
{
 public:
 struct CoordNames
 {
  MString lonname, latname, depthname, timename;
  size_t  nx, ny, nz, nt;
 };

 static MString                             StName2Name(const MString& stname);
 static MString                             Name2StName(const MString& name);
 static MString                             Name2LongName(const MString& name);
 static std::tuple<MDateTime, time_t, bool> Refdate(const MString& refdate);
 static void                                GetVars(const NCFileA& nc, std::set<MString>& vars);
 static CoordNames                          GetCNames(const NCFileA& nc);
 static CoordNames                          GetDNames(const NCFileA& nc);
 static bool                                HaveVar(const NCFileA& nc, const MString& vname);

 template<class NcZarrFunctions> static void       GetVars(const NcZarrRead<NcZarrFunctions>& nc, std::set<MString>& vars);
 template<class NcZarrFunctions> static CoordNames GetCNames(const NcZarrRead<NcZarrFunctions>& nc);
 template<class NcZarrFunctions> static CoordNames GetDNames(const NcZarrRead<NcZarrFunctions>& nc);
 template<class NcZarrFunctions> static bool       HaveVar(const NcZarrRead<NcZarrFunctions>& nc, const MString& vname);

 template<class HV> static VarPresence CheckVar(const MString& vname, HV hv)
 {
  if(!hv(vname))
  {
   auto varexist = VarPresence::NONE;
   if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = VarPresence::DERIVED;
   if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = VarPresence::DERIVED;
   if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = VarPresence::DERIVED;
   if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = VarPresence::DERIVED;
   if((vname == "ugeo" || vname == "vgeo") && hv("ssh")) varexist = VarPresence::DERIVED;
   return varexist;
  }
  return VarPresence::INTERNAL;
 }

 template<class D> static bool TransformationRead(const D* data, const MString& vname, std::map<MString, typename D::Data>& cache, const BaseParameters* ip, size_t i)
 {
  if constexpr(requires(const D* d, const BaseParameters* p) {
                {
                 d->Depth(p)
                } -> std::convertible_to<real>;
               })
  {
   real depth = data->Depth(ip);
   // ptemp from temp and sal
   if(vname == "ptemp")
   {
    if(!(data->Read("temp", cache, ip, i) && data->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))
      ptemp.V(ind) = ptemp.Fillval();
     else
      ptemp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), depth, ptemp.Lon(ind), ptemp.Lat(ind));
    }
    return true;
   }
   // temp from ptemp and sal
   if(vname == "temp")
   {
    if(!(data->Read("ptemp", cache, ip, i) && data->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(ptemp.IsFill(ind) || sal.IsFill(ind))
      temp.V(ind) = temp.Fillval();
     else
      temp.V(ind) = PTemp2Temp(ptemp.V(ind), sal.V(ind), depth, temp.Lon(ind), temp.Lat(ind));
    }
    return true;
   }
   // pdens from temp and sal
   if(vname == "pdens")
   {
    bool tempispot = data->CheckVar("ptemp") == VarPresence::INTERNAL;
    if(!(data->Read(tempispot ? "ptemp" : "temp", cache, ip, i) && data->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))
      pdens.V(ind) = pdens.Fillval();
     else
      pdens.V(ind) =
          tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), depth, pdens.Lon(ind), pdens.Lat(ind)) : Temp2PDens(temp.V(ind), sal.V(ind), depth, pdens.Lon(ind), pdens.Lat(ind));
    }
    return true;
   }
  }

  // U and U2 from u and v
  if(vname == "U" || vname == "U2")
  {
   bool square = vname == "U2";
   if(!(data->Read("u", cache, ip, i) && data->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();
    else
     U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
   }
   return true;
  }
  // ugeo
  if(vname == "ugeo")
  {
   if(!data->Read("ssh", cache, ip, i)) return false;
   const auto& ssh = cache.at("ssh");
   const real  mul = (ssh.Unit() == "cm") ? 1.0 : ((ssh.Unit() == "m") ? 100.0 : -1.0);
   if(mul < 0.0) return false; // Unknown units

   cache[vname] = cache.at("ssh").CopyGrid();
   auto& ugeo   = cache.at(vname);
   ugeo.SetUnit("cm/s");

   for(size_t iy = 0; iy < ssh.Ny(); iy++)
   {
    const real phi   = michlib::M_PI * ssh.Iy2Lat(iy) / 180.0;
    const real sphi  = michlib::Sin(phi);
    const real s2phi = michlib::Sin(2.0 * phi);
    const real g     = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi);                    // g in cm/s
    const real f     = 2.0 * 7.29211 * sphi;                                                                    // Coriolis parameter multiplied by 100'000
    const real db = (iy == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Iy2Lat(iy) - ssh.Iy2Lat(iy - 1)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
    const real du =
        (iy == ssh.Ny() - 1) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Iy2Lat(iy + 1) - ssh.Iy2Lat(iy)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter

    for(size_t ix = 0; ix < ssh.Nx(); ix++)
    {
     const bool bok = iy != 0 && !ssh.IsFill(ix, iy - 1) && !ssh.IsFill(ix, iy);
     const bool uok = iy != ssh.Ny() - 1 && !ssh.IsFill(ix, iy + 1) && !ssh.IsFill(ix, iy);
     if(!(bok || uok))
      ugeo(ix, iy) = ugeo.Fillval();
     else if(!bok)
      ugeo(ix, iy) = -g * mul * (ssh(ix, iy + 1) - ssh(ix, iy)) / (f * du);
     else if(!uok)
      ugeo(ix, iy) = -g * mul * (ssh(ix, iy) - ssh(ix, iy - 1)) / (f * db);
     else
     {
      const real bder = -g * mul * (ssh(ix, iy) - ssh(ix, iy - 1)) / (f * db);
      const real uder = -g * mul * (ssh(ix, iy + 1) - ssh(ix, iy)) / (f * du);
      ugeo(ix, iy)    = 0.5 * (bder + uder);
     }
    }
   }
   return true;
  }
  // vgeo
  if(vname == "vgeo")
  {
   if(!data->Read("ssh", cache, ip, i)) return false;
   const auto& ssh = cache.at("ssh");
   const real  mul = (ssh.Unit() == "cm") ? 1.0 : ((ssh.Unit() == "m") ? 100.0 : -1.0);
   if(mul < 0.0) return false; // Unknown units

   cache[vname] = cache.at("ssh").CopyGrid();
   auto& vgeo   = cache.at(vname);
   vgeo.SetUnit("cm/s");

   for(size_t iy = 0; iy < ssh.Ny(); iy++)
   {
    const real phi   = michlib::M_PI * ssh.Iy2Lat(iy) / 180.0;
    const real sphi  = michlib::Sin(phi);
    const real cphi  = michlib::Cos(phi);
    const real s2phi = michlib::Sin(2.0 * phi);
    const real g     = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s
    const real f     = 2.0 * 7.29211 * sphi;                                                 // Coriolis parameter multiplied by 100'000

    for(size_t ix = 0; ix < ssh.Nx(); ix++)
    {
     const real dl =
         (ix == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Ix2Lon(ix) - ssh.Ix2Lon(ix - 1)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
     const real dr = (ix == ssh.Nx() - 1)
                         ? 0.0
                         : (michlib::M_PI * 6371.0 * (ssh.Ix2Lon(ix + 1) - ssh.Ix2Lon(ix)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter

     const bool lok = ix != 0 && !ssh.IsFill(ix - 1, iy) && !ssh.IsFill(ix, iy);
     const bool rok = ix != ssh.Nx() - 1 && !ssh.IsFill(ix + 1, iy) && !ssh.IsFill(ix, iy);
     if(!(lok || rok))
      vgeo(ix, iy) = vgeo.Fillval();
     else if(!lok)
      vgeo(ix, iy) = g * mul * (ssh(ix + 1, iy) - ssh(ix, iy)) / (f * dr);
     else if(!rok)
      vgeo(ix, iy) = g * mul * (ssh(ix, iy) - ssh(ix - 1, iy)) / (f * dl);
     else
     {
      const real lder = g * mul * (ssh(ix, iy) - ssh(ix - 1, iy)) / (f * dl);
      const real rder = g * mul * (ssh(ix + 1, iy) - ssh(ix, iy)) / (f * dr);
      vgeo(ix, iy)    = 0.5 * (lder + rder);
     }
    }
   }
   return true;
  }
  return false;
 }
};