michael_ul
/
odm
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Added geostrophic module

master
Michael Uleysky 2 weeks ago
parent
bba5f9c054
commit
6c29188436
4 changed files with 199 additions and 5 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 9
      actions/actioninfo.h
  2. 7
      actions/actiontsc.h
  3. 15
      include/CF.h
  4. 173
      src/CF.cpp

9
actions/actioninfo.h
Unescape View File

@ -24,10 +24,15 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data)
ad = std::move(ret.Value());
}
{ // Modules
// Modules
{
auto ret = CF::TEOS10Module(args, pars, ad);
if(!ret) return "Error adding TEOS-10 derived variables";
}
{
auto ret = CF::GeostrophicModule(args, pars, ad);
if(!ret) return "Error adding geostrophic velocity";
}
const MString titCol = F(C::BWHITE);
const MString timeCol = F(C::YELLOW);
@ -72,7 +77,7 @@ MString ActionInfo_DoAction(const CLArgs& args, D& data)
// Metainfo
message(ad.Is2D(name) ? v2dCol : v3dCol, name, C(), info->lname.Exist() ? (" - " + info->lname) : "",
info->targetunit.Exist() ? (", measured in " + unitCol + info->targetunit + C()) : "");
if(info->comment.Exist()) message(comCol," ", info->comment, C());
if(info->comment.Exist()) message(comCol, " ", info->comment, C());
if(info->stname.Exist()) message(" Standard name: " + info->stname);
// Vertical grid description

7
actions/actiontsc.h
Unescape View File

@ -29,10 +29,15 @@ MString ActionTSC_DoAction(const CLArgs& args, D& ds)
ad = std::move(ret.Value());
}
{ // Modules
// Modules
{
auto ret = CF::TEOS10Module(args, pars, ad);
if(!ret) return "Error adding TEOS-10 derived variables";
}
{
auto ret = CF::GeostrophicModule(args, pars, ad);
if(!ret) return "Error adding geostrophic velocity";
}
MString name = args.contains("out") ? args.at("out") : "out.nc";
int compress = args.contains("compress") ? args.at("compress").ToInt() : 3;

15
include/CF.h
Unescape View File

@ -28,6 +28,19 @@ class CF
virtual ~VInfoTEOS10() override = default;
};
struct VInfoGeostrophic: public Adapter::VInfoBase
{
enum Operation
{
U,
V
};
MString sshvar;
Operation op;
virtual ~VInfoGeostrophic() override = default;
};
// Get time variable name
static MString GetTVarName(const NCZarr& nc);
@ -94,6 +107,8 @@ class CF
}
static RetVal<std::shared_ptr<Data3D>> TEOS103DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it);
static RetVal<std::shared_ptr<Data2D>> Geostrophic2DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, size_t it);
static Error TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad);
static Error GeostrophicModule(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad);
};

173
src/CF.cpp
Unescape View File

@ -456,7 +456,7 @@ Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListE
RetVal<std::shared_ptr<Data3D>> CF::TEOS103DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it)
{
static const MString pref = "CF::TempSalDer3DReader";
static const MString pref = "CF::Teos103DReader";
const auto* v = dynamic_cast<const struct VInfoTEOS10*>(vinfo);
if(v == nullptr) return {pref, "Incorrect argument type"};
@ -550,9 +550,119 @@ RetVal<std::shared_ptr<Data3D>> CF::TEOS103DReader(Adapter& ad, const Adapter::V
return out;
}
RetVal<std::shared_ptr<Data2D>> CF::Geostrophic2DReader(Adapter& ad, const Adapter::VInfoBase* vinfo, std::shared_ptr<Projection> proj, size_t it)
{
static const MString pref = "CF::Geostrophic2DReader";
const auto* v = dynamic_cast<const struct VInfoGeostrophic*>(vinfo);
if(v == nullptr) return {pref, "Incorrect argument type"};
std::shared_ptr<Data2D> out(new Data2D(proj, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment}));
std::shared_ptr<Data2D> ssh;
{
auto retssh = ad.Read2D(v->sshvar, ad.CurrentTime());
if(!retssh) return {pref, "Can't read ssh variable " + v->sshvar};
ssh = retssh.Value();
}
if(ssh->Ny() != out->Ny()) return {pref, MString("For ssh variable Ny=") + ssh->Ny() + ", for output variable Ny=" + out->Ny()};
if(ssh->Nx() != out->Nx()) return {pref, MString("For ssh variable Nx=") + ssh->Nx() + ", for output variable Nx=" + out->Nx()};
const UtUnit infr(ssh->Unit()), into("cm");
const UtUnit outfr("cm/s"), outto(out->Unit());
const auto cnvin = Converter(infr, into);
const auto cnvout = Converter(outfr, outto);
if(!cnvin) return {pref, "Can't convert units from \"" + ssh->Unit() + "\" to \"cm\""};
if(!cnvout) return {pref, "Can't convert units from \"cm/s\" to \"" + ssh->Unit() + "\""};
const bool needsconvertin = infr != into;
const bool needsconvertout = outfr != outto;
auto In = [needsconvertin, &cnvin = cnvin](real r) -> real { return needsconvertin ? cnvin(r) : r; };
auto Out = [needsconvertout, &cnvout = cnvout](real r) -> real { return needsconvertout ? cnvout(r) : r; };
switch(v->op)
{
case(VInfoGeostrophic::U):
{
for(size_t iy = 0; iy < ssh->Ny(); iy++)
{
const real phi = michlib::M_PI * ssh->Lat(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->Lat(iy) - ssh->Lat(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->Lat(iy + 1) - ssh->Lat(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))
(*out)(ix, iy) = out->Fillval();
else if(!bok)
(*out)(ix, iy) = Out(-g * (In((*ssh)(ix, iy + 1)) - In((*ssh)(ix, iy))) / (f * du));
else if(!uok)
(*out)(ix, iy) = Out(-g * (In((*ssh)(ix, iy)) - In((*ssh)(ix, iy - 1))) / (f * db));
else
{
const real bder = -g * (In((*ssh)(ix, iy)) - In((*ssh)(ix, iy - 1))) / (f * db);
const real uder = -g * (In((*ssh)(ix, iy + 1)) - In((*ssh)(ix, iy))) / (f * du);
(*out)(ix, iy) = Out(0.5 * (bder + uder));
}
}
}
break;
}
case(VInfoGeostrophic::V):
{
for(size_t iy = 0; iy < ssh->Ny(); iy++)
{
const real phi = michlib::M_PI * ssh->Lat(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->Lon(ix) - ssh->Lon(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->Lon(ix + 1) - ssh->Lon(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))
(*out)(ix, iy) = out->Fillval();
else if(!lok)
(*out)(ix, iy) = Out(g * (In((*ssh)(ix + 1, iy)) - In((*ssh)(ix, iy))) / (f * dr));
else if(!rok)
(*out)(ix, iy) = Out(g * (In((*ssh)(ix, iy)) - In((*ssh)(ix - 1, iy))) / (f * dl));
else
{
const real lder = g * (In((*ssh)(ix, iy)) - In((*ssh)(ix - 1, iy))) / (f * dl);
const real rder = g * (In((*ssh)(ix + 1, iy)) - In((*ssh)(ix, iy))) / (f * dr);
(*out)(ix, iy) = Out(0.5 * (lder + rder));
}
}
}
break;
}
}
return out;
}
Error CF::TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad)
{
static const MString pref = "TempSalDerModule";
static const MString pref = "TEOS10Module";
MString sal, temp, ptemp;
@ -708,3 +818,62 @@ Error CF::TEOS10Module(const CLArgs& args, michlib_internal::ParameterListEx& pa
return {};
}
Error CF::GeostrophicModule(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad)
{
static const MString pref = "GeostrophicModule";
MString ssh;
for(const auto& v: ad.Vars())
{
if(!ad.Is2D(v.first)) continue;
const auto& info = v.second.info;
if(!ssh.Exist() && info->stname == "sea_surface_height_above_geoid") ssh = v.first;
if(!ssh.Exist() && info->stname == "sea_surface_elevation") ssh = v.first;
}
if(!ssh.Exist()) return {};
const auto& svar = ad.Vars().at(ssh);
const auto& sproj = svar.proj;
const auto& rinfo = svar.info->rinfo;
if(!sproj->IsRectangular()) return {}; // We deal only with rectangular grids for now
{ // U-component
auto vinfo = new VInfoGeostrophic();
vinfo->sshvar = ssh;
vinfo->op = VInfoGeostrophic::U;
vinfo->stname = "surface_geostrophic_eastward_sea_water_velocity";
vinfo->lname = "Geostrophic eastward velocity";
vinfo->comment = "Module Geostrophic: Geostrophic eastward velocity calculated from sea surface height (" + ssh + ")";
vinfo->targetunit = "cm/s";
vinfo->rinfo = rinfo;
vinfo->fillval = NAN;
const bool addsuc = ad.Add2DVariable("ugs", Adapter::VInfo(vinfo), sproj, Geostrophic2DReader);
if(!addsuc) return {pref, "Variable with name ugs already exist"};
}
{ // V-component
auto vinfo = new VInfoGeostrophic();
vinfo->sshvar = ssh;
vinfo->op = VInfoGeostrophic::V;
vinfo->stname = "surface_geostrophic_northward_sea_water_velocity";
vinfo->lname = "Geostrophic northward velocity";
vinfo->comment = "Module Geostrophic: Geostrophic northward velocity calculated from sea surface height (" + ssh + ")";
vinfo->targetunit = "cm/s";
vinfo->rinfo = rinfo;
vinfo->fillval = NAN;
const bool addsuc = ad.Add2DVariable("vgs", Adapter::VInfo(vinfo), sproj, Geostrophic2DReader);
if(!addsuc) return {pref, "Variable with name vgs already exist"};
}
return {};
}

Write Preview
Loading…
Cancel
Save