From e89e52b269acfb8ef2240fc239762452d3198858 Mon Sep 17 00:00:00 2001 From: Michael Uleysky Date: Mon, 10 Apr 2023 15:28:20 +1000 Subject: [PATCH] Added support for conversion from in-situ temperature to potential temperature and vice versa, as well as potential density calculations. --- include/basedata.h | 2 + include/gsw.h | 25 ++++++++++++ include/layereddata.h | 3 +- src/layereddata.cpp | 95 ++++++++++++++++++++++++++++++++++++++++++- src/odm.cpp | 2 +- 5 files changed, 123 insertions(+), 4 deletions(-) create mode 100644 include/gsw.h diff --git a/include/basedata.h b/include/basedata.h index 3db88a2..3bfac64 100644 --- a/include/basedata.h +++ b/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; } diff --git a/include/gsw.h b/include/gsw.h new file mode 100644 index 0000000..c1422b9 --- /dev/null +++ b/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); } diff --git a/include/layereddata.h b/include/layereddata.h index d75eb91..2d03338 100644 --- a/include/layereddata.h +++ b/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 @@ -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"; diff --git a/src/layereddata.cpp b/src/layereddata.cpp index 3ff2097..14dd073 100644 --- a/src/layereddata.cpp +++ b/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 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(); 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(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(¶ms, i); + params.varname = "sal"; + auto sal = Read(¶ms, 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(¶ms, i); + params.varname = "sal"; + auto sal = Read(¶ms, 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(¶ms, i); + params.varname = "sal"; + auto sal = Read(¶ms, 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) { diff --git a/src/odm.cpp b/src/odm.cpp index 610156f..1dd7105 100644 --- a/src/odm.cpp +++ b/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");